Delta function at the origin in polar coordinates

I have some problems understanding what the best way of dealing with the delta functions in polar coordinates (I know there are many questions on the subjects on this website but they are all not satisfactory).

In (Delta function integrated from zero), they claim that the delta function is given by $\delta^{(2)}=\frac{\delta(r)}{\pi r}$

while in (Dirac delta in polar coordinates) it is claimed that $\delta^{(2)}=\frac{\delta(r)}{2\pi r}$.

However, the confusion probably comes from the fact that when evaluating a delta function in polar coordinates, one ends up with the expression $\int_0^\infty f(x)\delta(x)$. This expression is ill-defined as far as I can tell, since using different limiting functions for the delta function can give different results, and thus none of the above expressions can be a well-defined definition of the delta function in polar coordinates.

So my question is, if I want to write down the delta function in polar coordinates, what is the best representation for working with it? In my particular case, I want to be able to start with the delta function in polar coordinates and then do coordinate transformations to obtain it in other coordinate systems, without any ambiguities.

edit: The best representation I can come up with would be to regularize the radial direction, and write the delta function as $\delta=\frac{1}{r}\delta(r-\epsilon)\delta(\theta-\theta_0)$ for some arbitrary $\theta_0$ and then let $\epsilon\rightarrow0$ in the end.

see http://mathworld.wolfram.com/DeltaFunction.html eqn 46. The result given there corresponds to your first equation, $\delta^{(2)}=\frac{\delta(r)}{\pi r}$. However it can be more complicated: $\delta^{(2)}=\frac{\delta(r)}{\pi r}$ is only for a Dirac Delta located at the origin.
See the pdf at https://www.google.com/#q=06_notes_2dfunctions page 18. This shows the result given in your first equation is for an Dirac Delta at the origin, but your final equation $\delta=\frac{1}{r}\delta(r-\epsilon)\delta(\theta-\theta_0)$ represents an Dirac Delta radially offset from the origin by $\epsilon$ and rotated through the angle $\theta_0$.
So your equation should work for you, maybe rewritten $\delta(r-r_0)=\frac{1}{r_0}\delta(r-r_0)\delta(\theta-\theta_0)$ where $\epsilon$ is replaced with $r_0$ and no limiting process is needed.