Derivatives with respect to Cartesian coordinates of functions in spherical coordinates

Question

This seems to me to be a common task, and the more surprised I was to not have found much on this.

Suppose I have an expression f in spherical coordinates r,theta,phi, and want to compute partial derivatives of f with respect to Cartesian coordinates x,y,z, but with the result expressed in r,theta,phi again.

I have solved it by defining the following function (in Mathematica 10):

myD[f_,n_Integer/;n>=1&&n<=3]:=
  With[{jac=Transpose@CoordinateTransformData["Spherical"->"Cartesian","InverseMappingJacobian",{r,theta,phi}]},
    jac[[n]].(D[f,#]&/@{r,theta,phi})
]

With this I can compute $\frac{\partial f}{\partial z}$ by writing myD[f,3]. Take for example f=r, you should get: $$ \frac{\partial f}{\partial z}=myD[f,3]=\frac{\partial r}{\partial z}=\frac{z}{r}=\frac{r\cos\theta}{r}=\cos\theta $$

I was wondering if this is really the way to go for such a mundane task? I suspect that there is a nicer functionality built into Mathematica which I am not aware of. Does anyone have an idea?