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.
I'll take a different approach here. Suppose $\vec{\varepsilon}$ is fixed vector, or constant vector field with respect to the cartesian frame. Then as we write $\vec{\varepsilon}$ in the spherical frame the coefficients manifest a point-dependence. In particular, there exist functions $\varepsilon_{\rho},\varepsilon_{\phi},\varepsilon_{\theta}$ such that
$$\vec{\varepsilon} = \varepsilon_{\rho}\widehat{e}_{\rho}+\varepsilon_{\phi}\widehat{e}_{\phi}+\varepsilon_{\theta}\widehat{e}_{\theta} $$
we can calculate these coefficient functions by $\varepsilon_{\rho}=\vec{\varepsilon} \cdot \widehat{e}_{\rho}$, $\varepsilon_{\phi}=\vec{\varepsilon} \cdot \widehat{e}_{\phi}$ and $\varepsilon_{\theta} =\vec{\varepsilon} \cdot \widehat{e}_{\theta}$. This follows from the orthonormality of the spherical frame. (you can check, and I think you already realize $\widehat{e}_{\rho} \cdot \widehat{e}_{\rho}=1, \widehat{e}_{\rho} \cdot \widehat{e}_{\phi}=0$ etc...). However, you should also realize that $\{ \widehat{e}_{\rho}, \widehat{e}_{\phi}, \widehat{e}_{\theta} \}$ forms a right-handed-triple in the sense that the cross-products of $\widehat{e}_{\rho}, \widehat{e}_{\phi}, \widehat{e}_{\theta}$ share the same patterns as that of the standard cartesian frame:
$$ \widehat{e}_{\rho} \times \widehat{e}_{\phi} = \widehat{e}_{\theta}, \ \
\widehat{e}_{\phi} \times \widehat{e}_{\theta} = \widehat{e}_{\rho}, \ \
\widehat{e}_{\theta} \times \widehat{e}_{\rho} = \widehat{e}_{\phi} $$
Now, as you point out, $\vec{x} = \rho \widehat{e}_{\rho}$ thus,
$$ \vec{F} = \vec{\varepsilon} \times \vec{x} = (\varepsilon_{\rho}\widehat{e}_{\rho}+\varepsilon_{\phi}\widehat{e}_{\phi}+\varepsilon_{\theta}\widehat{e}_{\theta} ) \times \rho \widehat{e}_{\rho} = -\rho \varepsilon_{\phi}\widehat{e}_{\theta}
+ \rho \varepsilon_{\theta}\widehat{e}_{\phi}$$
now, we just need to calculate those dot-products and we're done.
That said, I'd rather use some notation like $\vec{A}$ instead of $\vec{ \varepsilon}$ since $e$ and $\varepsilon$ look so similar. In $\vec{A}$-notation,
$$ \vec{F} = \vec{A} \times \vec{x} = (A_{\rho}\widehat{e}_{\rho}+A_{\phi}\widehat{e}_{\phi}+A_{\theta}\widehat{e}_{\theta} ) \times \rho \widehat{e}_{\rho} = -\rho A_{\phi}\widehat{e}_{\theta}
+ \rho A_{\theta}\widehat{e}_{\phi}$$
Best Answer
A general way to proceed would be to expand the Green Function in spherical harmonics. To that end, we have
$$G(\vec r|\vec r')=\frac{1}{4\pi|\vec r-\vec r'|}=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{1}{2l+1}\left(\frac{r_<^l}{r>^{l+1}}\right)Y^*_{lm}(\theta',\phi')Y_{lm}(\theta,\phi)$$
where $r_<$ ($r_>$) is the smaller (larger) of $r$ and $r'$ and the spherical harmonic functions $Y_{lm}$ are given by
$$Y_{lm}(\theta,\phi)=\sqrt{ \frac{(2l+1)}{4\pi} \frac{(l-m)!}{(l+m)!} }P_l^m(\cos \theta)e^{im\phi}$$
Here, $P_l^m$ are the associated Legendre Polynomials which satisfy the orthogonality relationship
$$\int_{-1}^1 P_l^m(x)P_{l'}^{m'}(x)dx= \frac{2}{2l+1}\frac{(l+|m|)!}{(l-|m|)!}\delta_{ll'}$$
Note that the proposed function $f$ is independent of the azimuthal angle $\phi$, which renders the integration over $\phi$ trivial. It will reduce the problem to an integration over $r$ and $\theta$, which involves a single series over $l$ with only the Legendre polynomials $P_l^0(\cos \theta)=P_l(\cos \theta)$. Can you take it from here?