I'm going to give a rough outline that you can use to remember this, rather than a derivation. In an orthogonal coordinate system, the scale factors determine all the formulae for differential operators: these come from the definition
$$ \frac{\partial \mathbf{r}}{\partial q_i} = h_i \mathbf{e}_i, $$
where $\mathbf{e}_i$ is a unit vector, and $h_i$ is a function taken as positive. We set it up so that the $\mathbf{e}_i$ at a point are orthogonal (this is actually what it means to say the coordinate system is orthogonal),
$$\mathbf{e}_i \cdot \mathbf{e}_i = \delta_{ij}.$$
Then the line element looks like
$$ d\mathbf{r} = \sum_i \frac{\partial \mathbf{r}}{\partial q_i} \, dq_i = \sum_i \mathbf{e}_i h_i \, dq_i. $$
Then we can define the gradient by using the definition $df = \operatorname{grad}{f} \cdot d\mathbf{r}$ (think in Cartesians: this makes sense by the chain rule formula):
$$ df = \sum_i \frac{\partial f}{\partial q_i} \, dq_i = \sum_i \frac{1}{h_i}\frac{\partial f}{\partial q_i} h_i \, dq_i = \sum_{i,j} \frac{1}{h_i}\frac{\partial f}{\partial q_i} \mathbf{e}_i \cdot \mathbf{e}_j h_j \, dq_j = \left( \sum_i \frac{\mathbf{e}_i}{h_i} \frac{\partial f}{\partial q_i} \right) \cdot d\mathbf{r}, $$
using the orthogonality of the $\mathbf{e}_i$, so
$$ \operatorname{grad}f = \sum_i \frac{\mathbf{e}_i}{h_i} \frac{\partial f}{\partial q_i}. $$
Right, now the divergence. Let's think about what we want. My way of remembering what the divergence does is to think of it as a sort of formal adjoint to the gradient operator, because
$$ \operatorname{div}(f\mathbf{u}) = \operatorname{grad}{f} \cdot \mathbf{u} + f \operatorname{div}{\mathbf{u}}, $$
so the divergence theorem gives
$$ \int_V \operatorname{grad}{f} \cdot \mathbf{u} \, dV = - \int_V f \operatorname{div}{\mathbf{u}} \, dV + \text{(boundary terms)}. $$
The volume element is a cuboid spanned by three line elements, so $dV = h_1 h_2 h_3 \, dq_1 \, dq_2 \, dq_3 $.
Let's just look at one term: the others will work in the same way by the symmetry in the operators and linearity: suppose $u = u_1 \mathbf{e}_i$, and then
$$ \int_V \operatorname{grad}{f} \cdot \mathbf{u} \, dV = \iiint \frac{1}{h_1}\frac{\partial f}{\partial q_i} \mathbf{e}_i \cdot \mathbf{e}_1 u_1 \, h_1 h_2 h_3 \, dq_1 \, dq_2 \, dq_3 \\
= \iiint \frac{\partial f}{\partial q_1} u_1 \, h_2 h_3 \, dq_1 \, dq_2 \, dq_3
$$
Now, apply integration by parts to the $q_1$ integral:
$$ \iiint \frac{\partial f}{\partial q_1} u_1 \, h_2 h_3 \, dq_1 \, dq_2 \, dq_3 = \text{(boundary terms)} - \iiint f \frac{\partial}{\partial q_1} (h_2 h_3 u_1) \, dq_1 \, dq_2 \, dq_3 $$
Okay, but now I want to get back to a $dV$, so multiply and divide by $h_1 h_2 h_3$ to obtain
$$ \iiint \frac{\partial f}{\partial q_1} u_1 \, h_2 h_3 \, dq_1 \, dq_2 \, dq_3 = \text{(boundary terms)} - \iiint f \left( \frac{1}{h_1 h_2 h_3} \frac{\partial}{\partial q_1} (h_2 h_3 u_1) \right) \, dV \\
= \text{(boundary terms)} - \iiint f \operatorname{div} \mathbf{u} \, dV, $$
so we conclude that the divergence is given by the expression
$$ \operatorname{div} \mathbf{u} = \frac{1}{h_1 h_2 h_3} \left( \frac{\partial}{\partial q_1} (h_2 h_3 u_1) + \frac{\partial}{\partial q_2} (h_3 h_1 u_2) + \frac{\partial}{\partial q_3} (h_1 h_2 u_3) \right) $$
Right, fine. So what is the divergence in cylindrical coordinates, like you actually asked? Well, we just have to find the scale factors. $\mathbf{r} = (r\cos{\theta},r\sin{\theta},z)$, so
$$ h_r \mathbf{e}_r = 1(\cos{\theta},\sin{\theta},0)\\
h_{\theta} \mathbf{e}_{\theta} = r(-\sin{\theta},\cos{\theta},0)\\
h_z \mathbf{e}_z = 1(0,0,1), $$
so the only nontrivial $h$ is $h_{\theta}=r$. Thus
$$ \operatorname{div} \mathbf{u} = \frac{1}{r} \left( \frac{\partial}{\partial r} (r u_r) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (r u_z) \right) = \frac{1}{r} \frac{\partial}{\partial r} (r u_r) + \frac{1}{r}\frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (u_z). $$
"...I feel like this is too good to be true". Your question is perfectly legit since changes of coordinates are always tricky for differential operators. But luckily there is a general and well known theory, dubbed "curvilinear coordinates".
The answer is: YES, it is correct, but for a scalar field only. Do not try to extend naively this formula to higher-rank fields. Example: the directional derivative of a vector field in cylindrical coordinates is much more complex than this. Take a look at:
https://link.springer.com/content/pdf/bbm%3A978-3-0348-8579-9%2F1.pdf
The general topic of your question is how "orthogonal curvilinear coordinates" work: this includes the standard cartesian, cylindrical and spherical coordinates.
Best Answer
Mass continuity equation. The first line of the system gives indeed the scalar conservation law $$ \partial_t h + \mathrm{div} (h \boldsymbol{v}) = 0 . $$ Using the vector calculus identity $\mathrm{div} (h \boldsymbol{v}) = h \,\mathrm{div}\, \boldsymbol{v} + \mathbf{grad}\, h\cdot \boldsymbol{v}$ and the definition of the material time derivative $D_t$, we have also $D_t h = -h \,\mathrm{div}\, \boldsymbol{v}$.
Momentum equation. The remaining lines give an equation on $h \boldsymbol v$. In the absence of fluid rotation and with a flat bottom surface, the conservation of momentum reads $D_t \boldsymbol{v} = -\mathbf{grad}(g h)$. After a bit of tensor calculus, one should arrive at $$ \partial_t (h \boldsymbol{v}) + \mathbf{div}(h \boldsymbol{v}\otimes \boldsymbol{v}) = -\mathbf{grad}(\tfrac12 gh^2) , $$ where $\otimes$ is the tensor product. Using the identity $\mathbf{grad}\psi = \mathbf{div}(\psi \boldsymbol{I})$ for the scalar field $\psi = \tfrac12 gh^2$ where $\boldsymbol{I}$ is the identity tensor, we can rewrite the above equation as a single conservation law.
Of course, it is possible to solve this system in polar coordinates by setting $x = r\cos\theta$, $y = r\sin\theta$. With $\boldsymbol{u} = h\boldsymbol{v}$ and $\boldsymbol{T} = \frac1h \boldsymbol{u}\otimes \boldsymbol{u} + \frac12 gh^2\boldsymbol{I} = \boldsymbol{T}^\top$, the above equations of motion lead to the system \begin{aligned} \partial_t h + \partial_r u_{r} + \partial_\theta (\tfrac{1}{r} u_{\theta}) &= -\tfrac{1}{r} u_r \\ \partial_t u_r + \partial_r T_{rr} + \partial_\theta (\tfrac{1}{r} T_{r\theta}) &= -\tfrac{1}{r} (T_{rr} - T_{\theta\theta}) \\ \partial_t u_\theta + \partial_r T_{r\theta} + \partial_\theta (\tfrac{1}{r}T_{\theta\theta}) &= -\tfrac{2}{r}T_{r\theta} \end{aligned} which may be written \begin{aligned} \partial_t h + \partial_r (hv_{r}) + \tfrac{1}{r} \partial_\theta ( hv_{\theta}) &= -\tfrac{h}{r} v_r \\ \partial_t (hv_r) + \partial_r (h{v_r}^2 + \tfrac12 gh^2) + \tfrac{1}{r}\partial_\theta ( hv_rv_\theta) &= -\tfrac{h}{r} ({v_r}^2 - {v_\theta}^2) \\ \partial_t (hv_\theta) + \partial_r (hv_rv_\theta) + \tfrac{1}{r} \partial_\theta (h{v_\theta}^2 + \tfrac12 gh^2) &= -\tfrac{2 h}{r}v_r v_\theta \end{aligned} You may have a look at (1) for examples with cylindrical symmetry, i.e. where $\boldsymbol v$ is along the radial direction and where the flow is independent on $\theta$. In that case, the previous system becomes $$ \begin{aligned} \partial_t h + \partial_r u_r &= -\tfrac{1}{r} u_r \\ \partial_t u_r + \partial_r T_{rr} &= - \tfrac{1}{r h} {u_r}^2 \end{aligned} $$ with $T_{rr} = \frac1h {u_r}^2 + \frac12 gh^2$. The speeds of sound of this $2\times 2$ hyperbolic system are $v_r \pm \sqrt{g h}$.
The previous case of cylindrical waves is pretty easy to solve numerically (cf. the mentioned article (1)), e.g. by using dedicated finite-volume methods. However, the full system in polar coordinates is not as easy to solve, mostly because of the $1/r$ factors in the flux derivatives w.r.t. $\theta$. These factors make this problem space dependent in a certain way.
(1) P. Glaister, "Shallow water flow with cylindrical symmetry", J. Hydraul. Res. 29:2 (1991), 219-227. doi:10.1080/00221689109499005