I was reading the answer from here, which is showing $$\nabla = \partial_r e_r + \frac{1}{r} \partial_\theta e_\theta.$$ The answer actually made the choice $$e_r = (\cos\theta, \sin\theta), \quad e_\theta = (-\sin\theta, \cos\theta)$$ I would like to clarify this step more naturally.
More generally, let $g(x,y) = (u,v)$ be the change of variable map. Given a real valued function $f(x,y)$, we know in the new variable $\tilde f (u,v)= f(g^{-1}(u,v))$.
Now given a vector field $$v(x,y) = v_1\frac{\partial}{\partial x} + v_2\frac{\partial}{\partial y}$$ we would like to calculate $\tilde v(u,v)$ in the new variable, I think this is given by $dg(v)$ following the change of variable formula $$dg\left(\frac{\partial}{\partial x}\right) = \frac{\partial u}{\partial x} \frac{\partial}{\partial u}+\frac{\partial v}{\partial x} \frac{\partial}{\partial v}$$ $$dg\left(\frac{\partial}{\partial y}\right) = \frac{\partial u}{\partial y} \frac{\partial}{\partial u}+\frac{\partial v}{\partial y} \frac{\partial}{\partial v}$$ Let us denote the coordinates of $v$ as $[v] = (v_1, v_2)$, we get $$[\tilde v(u,v)] = (Jg)(g^{-1}(u,v))[v(g^{-1}(u,v))].$$
I tried the following calculation. Given $f(x,y)$, $h(r,\theta) = (x,y)$ the change of variable map, and $\tilde f = f \circ h$.
Given the vector field: gradient of $\tilde f$ $$\tilde v(r, \theta) = \frac{\partial \tilde f}{\partial r}\frac{\partial }{\partial r}+ \frac{1}{r^2}\frac{\partial \tilde f}{\partial \theta}\frac{\partial }{\partial \theta}$$ Then we see that $$dh(\tilde v(r, \theta)) = \frac{\partial \tilde f}{\partial r}dh\left(\frac{\partial }{\partial r}\right)+ \frac{1}{r^2}\frac{\partial \tilde f}{\partial \theta}dh\left(\frac{\partial }{\partial \theta}\right)\\ = \frac{\partial \tilde f}{\partial r}\left(\cos\theta \frac{\partial }{\partial x}+ \sin\theta \frac{\partial }{\partial y}\right)+ \frac{1}{r^2}\frac{\partial \tilde f}{\partial \theta}\left(-{r}\sin\theta \frac{\partial }{\partial x}+ {r}\cos\theta \frac{\partial }{\partial y}\right)\\ = \left(\cos\theta\frac{\partial \tilde f}{\partial r} -\frac{1}{r}\sin\theta \frac{\partial \tilde f}{\partial \theta}\right) \frac{\partial }{\partial x}+ \left(\sin\theta\frac{\partial \tilde f}{\partial r} +\frac{1}{r}\cos\theta \frac{\partial \tilde f}{\partial \theta}\right) \frac{\partial }{\partial y}\\ = \frac{\partial f }{\partial x}\frac{\partial }{\partial x}+ \frac{\partial f }{\partial y}\frac{\partial }{\partial y}. $$
Now $e_r = dh\left(\frac{\partial }{\partial r}\right)$ and $e_\theta = \frac{1}{r} dh\left(\frac{\partial }{\partial \theta}\right)$