The Hilbert space counting which gets $2D-2$ is correct. When we think about parameter counting in the way you have in this question, we are implicitly assuming that the equations are sufficiently "generic" so that intersections work the way they do in linear algebra. This is not always the case, particularly when we consider equations with singularities or which are defined on spaces with boundaries. For an extreme example, in a normed $N$-dimensional real vector space, the equation $|\vec v|^2 = 0$ is a single equation, but it reduces the $N$-dimensional space down to a single (apparently $0$-dimensional) point.
When we write the equation $\text{Tr} \rho^2 = 1$, a slightly more complicated version of the same thing is happening. If you diagonalize the density operator, you will get a set of $D$ real eigenvalues $\lambda_i$. These must each be non-negative for a density operator. Additionally, we know that $\rho$ has unit trace, meaning that $\sum_i \lambda_i = 1$. This means that each $\lambda_i \in [0,1]$. Under these conditions, $\lambda_i \ge \lambda_i^2$ with equality only for $\lambda_i \in \{0,1\}$. Thus $\text{Tr} \rho^2 = \sum_i \lambda_i^2 \le \text{Tr} \rho = 1$, and the two are equal only if all the $\lambda_i$ are either $0$ or $1$, which means that exactly one is $1$ and the others are $0$. Note that this is just saying that $\rho$ is a projection onto the single eigenvector with eigenvalue $1$, meaning that $\rho = | \psi \rangle \langle \psi |$ for some $|\psi\rangle$.
Let us also note that, in general, the boundary of the set of mixed states is not the set of pure states. Being on the boundary means that (just) one inequality becomes an equality, which means we only need one $\lambda_i = 0$. Being a pure state is a much stronger condition. I think this may be part of your confusion as the boundary of the set of mixed states does have dimension $D^2 - 2$. The Bloch sphere is an unhelpful example in this case, because since the Hilbert space is only $2$ dimensional, one eigenvalue going to $0$ is equivalent to being a pure state, but for larger $D$ that is not true.
Note that this still looks like you only impose $D$ real equations, namely one per eigenvalue, meaning the naive dimension counting still appears to be wrong. Why is that? The answer is tied to the fact that our end result has a degeneracy; specifically we have $D-1$ eigenvectors of $\rho$ with eigenvalue $0$. Thus the system, described in this way, has a fictitious $U(D-1)$ symmetry rotating those vectors. If you apply such a rotation, the density operator does not change, but our naive counting would not realize that. We would think that we should subtract the dimensionality of $U(D-1)$, namely $(D-1)^2$. But this $U(D-1)$ does not act freely; a transformation which only changes the phase of a given $|\psi_i \rangle$ leaves $|\psi_i \rangle \langle \psi_i |$ invariant, meaning that any basis for the Hilbert space is stabilized by $U(1)^{D}$, so we actually have only $(D-1)^2 - D$ real equation redundancies. Now we can finally get the right counting:
$$(D^2 - 1) - D - ((D-1)^2 - D) = 2D-2.$$
The "trick" that makes this counting work, whereas the more naive counting fails, is that these equations can all be imposed on the affine space of trace-class unit trace operators without requiring positive semi-definiteness and still get the right set of pure states. Positive semidefiniteness is the set of inequalities that bounds the space of mixed states, and if we don't impose it we don't have the issues arising before where the solutions to the equations are on the boundary of the space.
