Dirac Delta in spherical coordinates, Gaussian representation.

Question

I am trying to understand something seemingly simple about the Dirac Delta distribution. I understand that for a delta, $\delta(x)$, we can write it as $$ \delta_\epsilon(x) = \lim_{\epsilon \rightarrow 0} \frac{1}{\sqrt{2\pi}\epsilon}\exp(-x^2/2\epsilon^2) $$ which I am comfortable with. Where I get tripped up as is the following for spherical coordinates: $$ \delta(\vec{r}) = \frac{\delta(r-r_0)}{r^2}\frac{\delta(\phi - \phi_0)}{1}\frac{\delta(\theta - \theta_0)}{\sin\theta}, $$ such that $\int_V \delta(\vec{r})d^3\vec{r} = 1$. Where I am confused is how to apply the Gaussian representation to the angular parts such that $$ \int_0^{2\pi} \delta_\epsilon(\phi - \phi_0)d\phi = 1, \\ \int_0^\pi \delta_\epsilon(\theta - \theta_0) d\theta = \int_0^\pi\delta_\epsilon(\cos\theta - \cos\theta_0)\sin\theta d\theta = 1, $$ where I have used $\delta(f(x)) = \delta(x) / |df/dx|$ above. The Gaussian over $r$ makes sense to me because I can integrate from 0 to $\infty$ and get the correct result, but I am lost on how to transform the angular parts. I've looked into wrapped distributions, but I am unsure of how to invoke them here. Could someone perhaps point me in the right direction? I am looking to write a code that utilizes this formulation, but I am quite confused on how to move forward. Thanks for any help!