I have here a surface whose curvature properties I want to study, represented in cylindrical coordinates: $$f(r,\theta) = r^2\cos4\theta$$ The problem, however, is that the parametrization is not differentiable (indeed, not defined) at exactly the point I want to look at - the origin! I have a plot here made by Mathematica:
Trying to use the usual conversion to rectangular coordinates doesn't work either - at best, I can get $f(x,y) = (x^2 + y^2)\cos\cos^{-1}\left(\frac{x}{\sqrt{x^2+y^2}}\right)$, which is also not differentiable at the origin. Now, I'm not sure the surface itself is actually $C^\infty$, but I'm pretty sure it's at least $C^2$. How can we find a parametrization that's differentiable at the origin?

Note that by setting $z = x + iy = re^{i\theta}$, your expression is equal to
$$ f(x,y) = \operatorname{Re} \left( \frac{z^4}{|z|^2} \right) = \frac{\operatorname{Re} \left( (x + iy)^4 \right)}{x^2 + y^2} = \frac{x^4 - 6x^2y^2 + y^4}{x^2+y^2} = \frac{8x^4}{x^2+y^2} + y^2 - 7x^2. $$
Now, a priori $f$ is not defined at $(0,0)$. However, the polar expression (and your plot) of $f$ makes it clear that $\lim_{(x,y) \to (0,0)} f(x,y) = 0$ so you can define $f(0,0) = 0$ and then $f$ becomes continuous. How about higher regularity? It is somewhat easier in my opinion to continue to work in polar coordinates so let me state two useful observations:
$$ (\nabla h)(r,\theta) = \frac{\partial f}{\partial r} \hat{e}_r + \frac{1}{r} \frac{\partial f}{\partial \theta} \hat{e}_{\theta} = kr^{k-1}p(\cos(\theta),\sin(\theta))\hat{e}_r + r^{k-1} \left( \frac{\partial p}{\partial v}(\cos(\theta),\sin(\theta))\cos(\theta) - \frac{\partial p}{\partial u}(\cos(\theta),\sin(\theta))\sin(\theta) \right)\hat{e}_{\theta} = r^{k-1} \left( kp(\cos(\theta)\hat{e}_x + \sin(\theta) \hat{e}_y) + \left( \frac{\partial p}{\partial v}\cos(\theta) - \frac{\partial p}{\partial u}\sin(\theta) \right)(-\sin(\theta) \hat{e}_x + \cos(\theta) \hat{e}_y) \right) = r^{k-1} \left( \left( kp\cos(\theta) + \frac{\partial p}{\partial u}\sin^2(\theta) - \frac{\partial p}{\partial v}\sin(\theta)\cos(\theta) \right) \hat{e}_x + \\ \left( kp\sin(\theta) + \frac{\partial p}{\partial v}\cos^2(\theta) - \frac{\partial p}{\partial u}\sin(\theta)\cos(\theta) \right)\hat{e}_y \right). $$
Since $k > 1$ and the expressions involving $\theta$ are bounded, we have $\lim_{(x,y) \to (0,0)} (\nabla h)(x,y) = (0,0)$ and the first observation implies that $h$ is continuously differentiable and $(\nabla h)(0,0) = (0,0)$.
In your case, applying it to $h = f$ we see that $f$ is continously differentiable at the origin and $(\nabla f)(0,0) = 0$. However, the second derivatives of $f$ are not continuous as we have
$$ \frac{\partial^2 f}{\partial x^2} = \frac{2(x^6 + 3x^4y^2 + 27x^2y^4 - 7y^6)}{(x^2+y^2)^3}. $$
The second derivative $\frac{\partial ^2 f}{\partial x^2}$ is homogeneous of degree zero. In polar coordinates, it has the form $\frac{\partial ^2 f}{\partial x^2}(r,\theta) = q(\theta)$ and so it depends only on the angle $\theta$. Since it is non-constant, it cannot be extended continuously to the origin.