As Lubos points out, you are allowed to make any substitutions you want as long as the new function is equivalent to the old one. In this case, substituting $u(\rho)=\rho^l e^{-\rho}v(\rho)$ is allowed because you can recover $u$ if you know $v$.
The question is, though, why would you choose exactly such a form? If you're just given it outright, it does look pretty random. However, there's method behind it.
Essentially, you've got an equation that's got very particular behaviour at both of its endpoints. Take, for example the $\rho\ll1$ limit: then the equation reads
$$
\frac{d^2u}{d\rho^2}\approx\frac{\ell(\ell+1)}{\rho^2}u,
$$
whose solutions are $u=\rho^l$ and $u=\rho^{-(l+1)}$ (which we reject as unphysical). Of course, this says nothing exact about the full solution but it does inform us about its asymptotic behaviour as $\rho\rightarrow0$.
Similarly, in the $\rho\gg l$ limit you have $\frac{d^2u}{d\rho^2}\approx u$, so $u=e^{-\rho}$ (as $u=e^\rho$ is unphysical). Again, this only gives information about the asymptotics about $\rho\rightarrow\infty$.
So what have our limit games gained us? Formally, we know nothing about the exact solution, but we do know that it has to behave something like $u=\rho^le^{-\rho}\times\mathrm{some\, function}$. Plus, this new function will be better behaved than $u$ at the endpoints - constant at the origin and up to polynomial at infinity.
Thus, we make an informed guess and substitute $u=\rho^le^{-\rho}v(\rho)$, transforming the equation to the new variables, and hoping as we do so that the result will be a simpler equation. As it happens, it is both simpler and better behaved, with better-behaved singularities, and its solutions are now simple enough (the complexity being absorbed into the prefactors just discussed) that they are in fact polynomials.
In general, the exact form of the substitution will depend on the exact problem. One usually tries to do this kind of asymptotic analysis to try and distil away, slowly, the problem's complexity.
To find this Fourier transform, you start with
$$\Psi(\vec{p})=\frac{1}{(2\pi \hbar)^{3/2}}\int{e^{-i\vec{p}\cdot\vec{r} / \hbar }\psi(\vec{r})\mathrm d\vec r},$$
as you've written, but you need to be a bit more careful after that.
To begin with, you can't assume that $\vec p$ is along the $z$ axis, because $\vec p$ has been given to you as the argument of $\Psi(\vec p)$, and you can't change it. What you can do, however, is set up the coordinate system for the $\vec r$ integration so that its axes simplify your life, i.e. with $z$ along the direction of $\vec p$.
(The reason you can do this, by the way, is because the ground state $\psi(\vec r)$ is spherically symmetric. If it's not symmetric, such as e.g. anything that's not an $s$ state, then you need to use a decomposition of the plane wave into partial (spherical) waves, like the one in Jackson, 2nd ed., p. 471, eq. (10.43).)
Once you've done that, then yes, you can write $\vec p = (0,0,p)$, and you can further use spherical polar coordinates $\vec r = r(\sin(\theta) \cos(\phi), \sin(\theta)\sin(\phi), \cos(\theta))$ for your position vector, write
\begin{align}
\Psi(\vec{p})
& =
\frac{1}{(2\pi \hbar)^{3/2}}
\int_0^\infty \int_0^\pi \int_0^{2\pi}
e^{-ipr\cos(\theta) / \hbar }\psi(\vec{r})
\:r^2\sin(\theta)\:\mathrm d\phi\mathrm d\theta\mathrm dr
\\& =
\frac{1}{(2\pi \hbar)^{3/2}}\sqrt{\frac{1}{\pi a_{0}^{3}}}
\int_0^\infty \int_0^\pi \int_0^{2\pi}
e^{-ipr\cos(\theta) / \hbar }e^{-r/a_0}
\:r^2\sin(\theta)\:\mathrm d\phi\mathrm d\theta\mathrm dr,
\end{align}
and integrate away.
Alternatively (and this is really for any future visitors who chance upon here looking for an actual answer for the title question), you can cheat and look up the answer, which is found in e.g.
The Momentum Distribution in Hydrogen-Like Atoms. B. Podolsky and L. Pauling. Phys. Rev., 34, 109 (1929).
(Not you, Marc, by the way. Go and do the integral.)
Best Answer
I've not really used Griffiths, but what I typically see is that one defines/derives the $2n=e^2k/E$ and $k=\sqrt{-2mE/\hbar^2}$ terms you have from the Schrodinger equation, then squares the $n$ term & substitutes, $$ 4n^2=\frac{e^4k^2}{E^2}=\frac{e^4}{E^2}\cdot-\frac{2mE}{\hbar^2} $$ to find that $E\propto1/n^2$. You can then use this definition of $E$ to determine $k$: $$ k=\sqrt{-\frac{2mE}{\hbar^2}}=\sqrt{-\frac{2m}{\hbar^2}\cdot-\frac{me^4}{2\hbar^2n^2}}=\frac{1}{an} $$ where $a=\hbar^2/me^2$ is the Bohr radius. So I would say it is not an assumption, but a result of the derivation. A complete derivation in this fashion can be found at Michael Fowler's page on the hydrogen atom.