If $X$ and therefore $X^T$ are indeed square and invertible then you are correct: the linear system $y = X\beta$ can be solved by inversion to get $\beta = X^{-1}y$ (and the residuals are all zero). But in linear regression we usually have $n \gg p$ so rather than exactly solving that system, we instead find the element of the column space of $X$ that is as close as possible to $y$, which likely is not in that $p$-dimensional subspace. This element is $\hat y = X\hat\beta$, and $\hat\beta$ is found via $\hat\beta = X^{+}y$ where $X^+$ is the pseudoinverse rather than the actual inverse.
The equation $(X^TX)^{-1} \stackrel ?= X^{-1}X^{-T}$ does not hold when $X$ is not square because $X$ is not invertible so that is not meaningful in general (as @whuber says in his comment).
For your update, we are using the fact that
$$
(ABC)^{-1}=C^{-1}B^{-1}A^{-1}
$$
for invertible $A$, $B$, and $C$, but we're still never pulling $X^TX$ apart, we're just doing
$$
((XD)^T(XD))^{-1}(XD)^T = \\
(DX^TXD)^{-1}DX \\
= D^{-1}(X^TX)^{-1}D^{-1}DX^T \\
= D^{-1}(X^TX)^{-1}X^T
$$
so once we start using inverses $X^TX$ is only ever considered as a single square matrix, the fact that it's coming from $X$ isn't relevant as far as the inversion is concerned.
