First, $\mathbf{F} = x\mathbf{\hat i} + y\mathbf{\hat j} + z\mathbf{\hat k}$ converted to spherical coordinates is just $\mathbf{F} = \rho \boldsymbol{\hat\rho} $. This is because $\mathbf{F}$ is a radially outward-pointing vector field, and so points in the direction of $\boldsymbol{\hat\rho}$, and the vector associated with $(x,y,z)$ has magnitude $|\mathbf{F}(x,y,z)| = \sqrt{x^2+y^2+z^2} = \rho$, the distance from the origin to $(x,y,z)$.
You also asked about where
$$\begin{bmatrix}\boldsymbol{\hat\rho} \\ \boldsymbol{\hat\theta} \\ \boldsymbol{\hat\phi} \end{bmatrix} = \begin{bmatrix} \sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\ \cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\ -\sin\phi & \cos\phi & 0 \end{bmatrix} \begin{bmatrix} \mathbf{\hat x} \\ \mathbf{\hat y} \\ \mathbf{\hat z} \end{bmatrix}$$
comes from. Let's look at the simpler 2D case first. For a point $(x,y)$ it helps to imagine that you're on a circle centered at the origin. In this case, the two fundamental directions you can move are perpendicular to the circle or along the circle. For the perpendicular direction we use the outward-pointing radial unit vector $\mathbf{\hat{r}}$. For the other direction, moving along the circle means (instantaneously) that you're moving tangent to it, and we take the unit vector in this case to be $\boldsymbol{\hat\theta}$, pointing counterclockwise. For example, suppose you're at the point $(1/\sqrt{2},1/\sqrt{2})$. Then, in the graph below, $\mathbf{\hat{r}}$ is in red and $\boldsymbol{\hat\theta}$ is in yellow.
Note that this means that, unlike the unit vectors in Cartesian coordinates, $\mathbf{\hat{r}}$ and $\boldsymbol{\hat{\theta}}$ aren't constant; they change depending on the value of $(x,y)$.
Now, what about a formula for $\mathbf{\hat{r}}$? If we move perpendicular to the circle we're keeping $\theta$ fixed in the polar coordinate representation $(r \cos \theta, r \sin \theta)$. The vector $\mathbf{\hat{r}}$ is the unit vector in the direction of this motion. If we interpret $r$ as time, taking the derivative with respect to $r$ will give us the velocity vector, which we know points in the direction of motion. Thus we want the unit vector in the direction of $\frac{d}{dr} (r \cos \theta, r \sin \theta) = (\cos \theta, \sin \theta)$. This is already a unit vector, so $\mathbf{\hat{r}} = \cos \theta \mathbf{\hat{x}} + \sin \theta \mathbf{\hat{y}} $. Similarly, moving counterclockwise along the circle entails keeping $r$ fixed in the polar coordinate representation $(r \cos \theta, r \sin \theta)$. Thus to find $\boldsymbol{\hat\theta}$ we take $\frac{d}{d\theta} (r \cos \theta, r \sin \theta) = (-r \sin \theta, r \cos \theta)$. This is not necessarily a unit vector, and so we need to normalize it. Doing so yields $\boldsymbol{\hat\theta}= -\sin \theta \mathbf{\hat{x}} + \cos \theta \mathbf{\hat{y}} $. In matrix form, this is $\begin{bmatrix} \mathbf{\hat{r}} \\ \boldsymbol{\hat\theta}\end{bmatrix} = \begin{bmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{bmatrix} \begin{bmatrix} \mathbf{\hat{x}} \\ \mathbf{\hat{y}} \end{bmatrix}$.
Moving up to spherical coordinates, for a given point $(x,y,z)$, imagine that you're on the surface of a sphere. The three fundamental directions are perpendicular to the sphere, along a line of longitude, or along a line of latitude. The first corresponds to $\boldsymbol{\hat\rho}$, the second to $\boldsymbol{\hat\theta}$, and the third to $\boldsymbol{\hat{\phi}}$. (This is using the convention in the Wikipedia page, which has $\theta$ and $\phi$ reversed from what you have.) Thus to find $\boldsymbol{\hat\rho}$, $\boldsymbol{\hat\theta}$, and $\boldsymbol{\hat{\phi}}$, we take the derivative of the spherical coordinate representation $(\rho \sin \theta \cos \phi, \rho \sin \theta \sin \phi, \rho \cos \theta)$ with respect to $\rho$, $\theta$, and $\phi$, respectively, and then normalize each one. That's where the matrix
$$\begin{bmatrix}\boldsymbol{\hat\rho} \\ \boldsymbol{\hat\theta} \\ \boldsymbol{\hat\phi} \end{bmatrix} = \begin{bmatrix} \sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\ \cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\ -\sin\phi & \cos\phi & 0 \end{bmatrix} \begin{bmatrix} \mathbf{\hat x} \\ \mathbf{\hat y} \\ \mathbf{\hat z} \end{bmatrix}$$
comes from.
As Henning Makholm points out, one way to view what we're doing here is that we're rotating the $\mathbf{\hat{x}}, \mathbf{\hat{y}}, \mathbf{\hat{z}}$ vectors. The transformation matrix can thus be considered a change-of-basis matrix. This means you could also (and more generally) convert $\mathbf{F} = x\mathbf{\hat i} + y\mathbf{\hat j} + z\mathbf{\hat k}$ to spherical coordinates via
\begin{align}
\mathbf{F}
&=
\begin{bmatrix} F_{\rho}\\ F_{\theta}\\ F_{\phi}\end{bmatrix}
\\
&=
\begin{bmatrix} \sin \theta \cos \phi & \sin \theta \sin \phi & \cos \theta \\ \cos \theta \cos \phi & \cos \theta \sin \phi & -
\sin \theta \\ - \sin \phi & \cos \phi & 0 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix}
\\
&=
\begin{bmatrix} \sin \theta \cos \phi & \sin \theta \sin \phi & \cos \theta \\ \cos \theta \cos \phi & \cos \theta \sin \phi & - \sin \theta \\ - \sin \phi & \cos \phi & 0 \end{bmatrix} \begin{bmatrix} \rho \sin \theta \cos \phi \\ \rho \sin \theta \sin \phi \\ \rho \cos \theta \end{bmatrix}
\\
&=
\begin{bmatrix} \rho \sin^2 \theta \cos^2 \phi + \rho \sin^2 \theta \sin^2 \phi + \rho \cos^2 \theta \\ \rho \sin \theta \cos \theta \cos^2 \phi + \rho \sin \theta \cos \theta \sin^2 \phi - \rho \sin \theta \cos \theta \\ - \rho \sin \theta \sin \phi \cos \phi + \rho \sin \theta \sin \phi \cos \phi + 0 \end{bmatrix}
\\
&= \begin{bmatrix} \rho \\ 0 \\ 0\end{bmatrix}.
\end{align}
So we get $\mathbf{F} = \rho \boldsymbol{\hat\rho} + 0 \boldsymbol{\hat\theta} + 0 \boldsymbol{\hat{\phi}} = \rho \boldsymbol{\hat{\rho}}$, just as we argued at the beginning.
Since you (the OP) haven't accepted an answer, I'm posting this, but consider this as a supplement to amd's answer, since his/her contribution made me understood this problem, about which I was recurrently thinking for two days. Also the comment of Jahan Claes helped, which was the same that one of my teachers told me. I would also like to note, that the reasoning given here is an intuitive, physicist's reasoning, perhaps not really fit for a mathematics site. However, I think that the question also refers to acquiring the volume element with this kind of reasoning, in an intuitive way.
Consider a point in 3D space described by coordinates $(r,\theta,\phi)$ such as in Fig. 1. 
Figure 1.
Now consider an infinitesimal increase $dr$ in $r$. Whatever the value of $theta$ and $phi$ are, the length of this increase in the geometrical space will be the same. See Fig. 2. 
Figure 2.
Now, let's forget about the increase in $r$ for a moment, and consider an infinitesimal increase in $\theta$, and let's just foucus on that. Whatever the value of $\phi$ is, it won't change the arc length, but it strongly depends on the value of $r$, namely as $rd\theta$. You can see this by intersecting the 3D space with a plane perpendicular to the $xy$ plane, and going through the line which connets the orgin $(0,0,0)$ and the point $(r,\theta,\phi)$. For illustration see Fig. 3.
Figure 3.
Let's forget about the increse in $\theta$ too, and let's just focus on what happens, if there is and infinitesimal increase $d\phi$ in $\phi$. To see what happens, I suggest looking on the setup from "above", as shown in Fig. 4. (the $z$ axis is pointing outside of the screen). You can see now, that if $\theta$ is smaller ($\theta_1$), the arclength due to the increase in $\phi$ will also be smaller. However, if $\theta$ is larger, ($\theta_2$), than from above you can see that the arclength will also be larger. If you determine the exact relation, you will see that the arc length is proportional to the projection of the position vector on the $xy$ plane, and $d\phi$. That is the arclength is $|Proj_{\theta}(\vec{r})|d\phi=r\sin\theta d\phi$. 
Figure 4.
Now, to put it all together, refer to Fig. 5. The volume element comes from the infinitesimal increase in $r$, $\phi$, and $\theta$. This is of course a trapesoidlike thing with pieces of spherical surfaces on its top and bottom. However, since it is infinitely small, it may be approximated by a cuboid (the precise justification of this approximation I can't do for the moment), and its volume may be written as the length of its sides multiplied by eachother, that is its volume is \begin{equation}dV=drrd\theta r\sin\theta d\phi=r^2\sin\theta drd\theta d\phi \end{equation} 
Figure 5.
Two sidenotes:
- I really like amd's idea, the same one that one of my teachers told me. That is the term $r^2\sin\theta$ which we have acquired intuitively, is the so called "volume distortion factor", which shows you that if you were to parametrize the geometrical space with three numbers denoted by $r$, $\theta$ and $\phi$, than the parameter space would have three axes denoted by the same letters. If in that paremeter space you would take an infinitesimal cuboid with sides $dr$, $d\theta$ and $d\phi$, with one of its corners residing at $(r,\theta,\phi)$, than the volume of this infinitesimal increment after transforming it back to the real, geometrical space would appear as $r^2\sin\theta drd\theta d\phi$. In genral this volume distortion factor is given by the Jacobian. For the case of transforming spherical coordinates into Descartes coordinates, the value of the Jacobian is indeed $r^2\sin\theta$.
- One may argue, that in my approximation of the strange trapesoid-like shape as a cuboid, I could have taken the length of the upper arches in computing its volume, that is its volume could have been $dV=dr*(r+dr)d\theta*(r+dr)sin\theta d\phi$. Or one could also argue why I have used $\theta$ in the calculations, not $\theta+d\theta$. That is the volume could also have been $dV=dr*(r+dr)d\theta*(r+dr)sin(\theta+d\theta) d\phi$. However the limiting process done while integrating makes those terms disappear, which contain the squares of infinitesimally small elements, if they appear in a sum that contains other infinitesimally small elements on a lower power. That is e.g. $const*dr+const*dr^2\approx const*dr$, $const*r+const*dr\approx const*r$, and so on. In the second case of computing the volume element, where I have used the infinitesimally increased $\theta$, one would also have to use a trigonometric identity on the sine, and use the approximations $\sin d\theta\approx d\theta$ and $\cos d\theta\approx 1$ to see this.
Best Answer
Your intuition did not fail you, the formula in your question was incorrect indeed. The correct formula for a 3D metric in polar coordinates is:
$${\rm d}s=\sqrt{{\rm d}r^2 + (r\,{\rm d}\theta)^2 + (r \sin(\theta) \,{\rm d}\phi)^2}$$
In case of an arc on the surface of a sphere, $\,{\rm d}r=0\,$ and you have:
$${\rm d}s=r\sqrt{{\rm d}\theta^2 + (sin(\theta) \,{\rm d}\phi)^2}$$
Both still are the same Pythagorian theorem, if you look at them closely.