Projection of point to plane (using least squares) I want to project 3D points to a plane. First, I was thinking of using the projection obtained by least squares (link)
$$\mathbf{\hat{z}} =  \mathbf {X} \left(\mathbf {X} ^{\textsf {T}}\mathbf {X} \right)^{-1}\mathbf {X} ^{\textsf {T}} \mathbf{{z}}  $$
based on the plane equation
$$ ax + by + d = z.$$
At that point I was wondering how to transform the $x$- and $y$-coordinate into the respective subspace?
Furthermore, I found this procedure and I'm wondering where the connection is.

The formula you quoted is, in a sense, the wrong tool for the job, but it can be used. I would instead add a multiple of the normal vector $(a, b, -1)$ to whatever you’re projecting, so that it lies on the plane. But, let’s do it with the LLS formula, as requested.
The first thing to note is that the formula is for projecting onto a linear subspace of $\Bbb{R}^n$, and the plane equation above is only a subspace if $d = 0$. This is not a problem; we just have to translate the problem. We take a particular point $\mathbf{v}_0 = (x_0, y_0, z_0)$ in the plane (it doesn’t matter which one; an easy choice would be $(0, 0, d)$), and then our projection of an arbitrary point $\mathbf{v}$ will be given by:
$$\mathbf{v_0} + \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top(\mathbf{v} - \mathbf{v}_0).$$
The $\mathbf{X}$ simply needs to be a matrix whose columns form a basis for the subspace we are projecting onto, i.e. $\{(x, y, z) : ax + by = z\}$ (we need the columns to be linearly independent, not just spanning, so $(\mathbf{X}^\top \mathbf{X})^{-1}$ exists, otherwise we’d need to compute the Moore-Penrose pseudoinverse). We just need to find two linearly independent vectors that satisfy $ax + by = z$, and there are two relatively natural choices:
$$(x, y, z) = (1, 0, a), (0, 1, b).$$
Of course, there are other possibilities too! But, this is the choice we’re using, so our $\mathbf{X}$ becomes
$$\mathbf{X} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ a & b \end{pmatrix}.$$
Just for fun, we then get
$$\mathbf{X}^\top\mathbf{X} = \pmatrix{1 + a^2 & ab \\ ab & 1 + b^2} \implies (\mathbf{X}^\top \mathbf{X})^{-1} = \frac{1}{1 + a^2 + b^2}\pmatrix{1 + b^2 & -ab \\ -ab & 1 + a^2}.$$
Thus,
\begin{align*}
\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top &= \frac{1}{1 + a^2 + b^2} \pmatrix{1 & 0 \\ 0 & 1 \\ a & b}\pmatrix{1 + b^2 & -ab \\ -ab & 1 + a^2} \pmatrix{1 & 0 & a \\ 0 & 1 & b} \\
&= \frac{1}{1 + a^2 + b^2} \pmatrix{1 & 0 \\ 0 & 1 \\ a & b}\pmatrix{1 + b^2 & -ab & a \\ -ab & 1 + a^2 & b} \\
&= \frac{1}{1 + a^2 + b^2} \pmatrix{1 + b^2 & -ab & a \\ -ab & 1 + a^2 & b \\ a & b & a^2 + b^2}
\end{align*}
Thus, the projection formula onto the original plane becomes:
$$\pmatrix{x \\ y \\ z} \mapsto \pmatrix{0 \\ 0 \\ d} + \frac{1}{1 + a^2 + b^2} \pmatrix{1 + b^2 & -ab & a \\ -ab & 1 + a^2 & b \\ a & b & a^2 + b^2}\pmatrix{x \\ y \\ z - d}.$$
