Projection of a 3D spherical distribution function in to a 2D cartesian plane Consider a 3D spherical Gaussian distribution function that depends on radius only,
$$f(r) = \frac{1}{N} e^{-(\frac{r-R_\mu}{\sigma})^2}$$
where $R_\mu$ is the radial offset of the distribution and $\sigma$ is the usual variance parameter. $N$ is the normalization constant. To help visualize this function, see the plot below which shows $f(r)$ plotted in one quadrant of the x-y plane, using $R_\mu = 10$ and $\sigma^2 =0.5$.

What is the projection of the distribution function $f(r)$ on to the x-y plane, $f(x,y)$?
Note that it is obviously not as simple as what I have plotted above which is simply $f(x,y,0)$, the slice through the x-y plane. I'm asking for the projection of this distribution function in to the x-y plane, $f(x,y)$.
A simpler example
To try and solve this problem I have tried to evaluate the simpler case of a uniform distribution on the surface of a sphere with radius $R_\mu$, i.e. the limit of $f(r)$ above as $\sigma$ goes to zero. I believe this sample distribution would be defined as
$$
f_s(r) = \begin{cases}
\frac{1}{N}, & \text{if }r = R_\mu \\
0, & \text{if }r \neq R_\mu
\end{cases}
$$
The integral over the whole surface must equal 1.
$$
\int_\Omega f_s(r) d\mathbf{a} = \int_0^{\pi} \int_0^{2\pi} f_s(r) r^2 Sin(\phi) d\theta d\phi = 1
$$
To find the projection of this much simpler function on to the xy plane, I express it in terms of an infinitesimal area element. In spherical coordinates the area element is $ d\mathbf{a} = r^2 Sin(\phi)$. But $r^2 = R_\mu$, and $Sin(\phi)$ can be re-written in terms of the x-y coordinates since $Sin (\phi ) = \frac{z}{R_\mu}$. e.g. from the equation for a sphere
$$x^2 + y^2 + z^2 = R_\mu ^2$$
$$\therefore \frac{z^2}{R_\mu^2} = 1- \frac{x^2+y^2}{R_\mu^2}$$
$$\therefore Sin(\phi) = \sqrt{1- \frac{x^2+y^2}{R_\mu^2}}$$
Hence, I find the projection of this distribution in to the x-y plane is:
$$ f_s(x,y) = \begin{cases}
\frac{1}{N} R_\mu^2 \sqrt{1- \frac{x^2+y^2}{R_\mu^2}}, & \text{if }x^2 + y^2 < R_\mu^2 \\
0, & \text{if }x^2 + y^2 > R_\mu^2
\end{cases} $$
Is this correct? If so, how do I do the more complex case above for a spherical Gaussian in 3D projected on to the xy plane? Do I need to express $f(r)$ in Cartesian coords and integrate the elemental volume over $z$? If I use this logic can I write it in 3D as
$$f(r) = \frac{1}{N} e^{-(\frac{r-R_\mu}{\sigma})^2}$$
$$ dV = r^2 Sin(\phi) dr d\theta d\phi$$
$$ f(x,y) = \int (x^2 + y^2 +z^2) \frac{z}{\sqrt{x^2 + y^2 +z^2}} e^{-(\frac{r-R_\mu}{\sigma})^2} dz  $$
I don't think I've got this last step correct because it doesn't look right when I evaluate it numerically and plot it. I suspect it needs to visually converge on the previous result as $\sigma$ -> 0. 
Any help you can give is appreciated. This is for a research project where I plan to use these distribution functions in another expression. I expect my question is a very specific case of someone's more general question.
EDIT: Extra graph inspired by Joriki's answer below
See Joriki's answer below. Its easy to see why it makes sense by generating 1000 random unit vectors and plotting their (x-y) components on a scatter plot. You get something like this. Its clear that my original result has the wrong behaviour.

 A: No, that's not correct. You can see that it can't be correct because it's zero on the boundary and maximal in the centre, whereas it should be minimal at the centre and maximal at the boundary where the slope of the surface becomes infinite and thus concentrates a lot of mass.
There are two errors in fact. One is that you seem to have mixed different conventions for the spherical coordinates – if $z$ is proportional to the sine of the polar angle, the volume element is proportional to the cosine, and vice versa; you've got the sine for both.
The other problem is that you've simply dropped $\mathrm d\phi\mathrm d\theta$ whereas you need to transform it to $\mathrm dx\mathrm dy$. Going with the $z\propto\cos\phi$ convention, we have
$$
\begin{align}
x&=r\sin\phi\cos\theta\;,\\
y&=r\sin\phi\sin\theta
\end{align}
$$
and thus
$$
\frac{\partial(x,y)}{\partial(\phi,\theta)}=(r\cos\phi\cos\theta\,r\sin\phi\cos\theta+r\sin\phi\sin\theta\,r\cos\phi\sin\theta)=r^2\sin\phi\cos\phi\;,
$$
so the $\sin\phi$ in $\sin\phi\,\mathrm d\phi\,\mathrm d\theta$ is divided out and you get a factor $1/\cos\phi$ instead, which diverges as you approach the boundary, as it should.
A: Answer for the 3D case based on @Joriki's response. The 3D distribution function of interest is
$$f(r)=\frac{1}{N} e^{-(\frac{r-R_\mu}{\sigma})^2}$$
The integral over all space in spherical coordinates must be 1.
$$1= \iiint_V f(r) r^2 Sin(\phi) d\phi d\theta dr$$
The Jacobian for the transformation back to cartesian is $\frac{\partial (r,\theta,\phi)}{\partial (x,y,z)} = \frac{1}{r^2 Sin(\phi)}$, hence
$$\iiint_V f(r) r^2 Sin(\phi) d\phi d\theta dr = \iiint_V f(r) r^2 Sin(\phi) \frac{\partial (r,\theta,\phi)}{\partial (x,y,z)}  dz dy dx$$
$$= \iiint_V \frac{1}{N} e^{-(\frac{\sqrt{x^2+y^2+z^2}-R_\mu}{\sigma})^2}  dz dy dx$$
Therefore we can pull out the integral over dz as this gives us $f(x,y)$.
$$\therefore f(x,y) = \frac{1}{N}\int_{-\infty}^{\infty} e^{-(\frac{\sqrt{x^2+y^2+z^2}-R_\mu}{\sigma})^2}  dz $$
We can see this is correct by plotting this function for several converging values of $\sigma$. We expect the function to converge on to $\frac{1}{Cos(\phi)}$. For the plot below I used $\sigma$ values of 2,1,0.75,0.5. $R_\mu = 10$

A: Nice, although it seems that in the 'simpler case' one could derive the $\frac{1}{\cos(\phi)}$ factor straight from the circle-to-line (rather than sphere-to-plane) projection: local curvature is ignored under the infinitesimal approximation, so the cylindrical case is equivalent to spherical, in the context.
