Define $g(r,\theta,\phi) := f(r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta)$.
Our assumptions say
$$\frac{\partial f}{\partial x}r\sin\theta\sin\phi = \frac{\partial f}{\partial y}r\sin\theta\sin\phi$$
$$\frac{\partial f}{\partial y}r\cos\theta = \frac{\partial f}{\partial z}r\sin\theta\sin\phi$$
We have
$$\frac{\partial g}{\partial \phi} = -\frac{\partial f}{\partial x}r\sin\theta\sin\phi + \frac{\partial f}{\partial y}r\sin\theta\sin\phi = 0$$
so $g$ in fact does not depend on $\phi$ and we can fix $\phi = \frac\pi2$ so
$$g(r,\theta) = f(0, r\sin\theta, r\cos\theta).$$
We have
$$\frac{\partial g}{\partial \theta} = \frac{\partial f}{\partial y}r\cos\theta - \frac{\partial f}{\partial z}r\sin\theta = 0$$
so $g$ doesn't depend on $\theta$ either. Therefore $g(r,\theta,\phi) = g(r)$ and is differentiable.
Taking $(x,y,z)$ from a dense subset of $\Bbb{R}^3$ we can define $$(r,\theta,\phi) = \left(\sqrt{x^2+y^2+z^2}, \arctan\frac{\sqrt{x^2+y^2}}z, \arctan\frac{y}x\right)$$
and notice
$$f\left(x,y,z\right) = f(r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta) = g(r,\theta,\phi) = g(r)= g(\|(x,y,z)\|).$$
By continuity the result extends to all of $\Bbb{R}^3$.