I calculated the circumference $S_1$ of a circle of radius $r$ using polar coordinates $(x, y)\rightarrow(\rho\cos\theta, \rho\sin\theta)$:
$$ \begin{align} S_1 &= \int\int_{\{(x, y)|x^2+y^2=r^2\}}dxdy \\ &=\int_0^{2\pi}\int_0^\infty\delta(\rho-r)|\frac{\partial(x, y)}{\partial(\rho, \theta)}|d\rho d\theta \\ &=\int_0^{2\pi}d\theta\int_0^\infty\delta(\rho-r)\rho d\rho \\ &= 2\pi r. \end{align} $$
If I replaced $t=r^2$, $(x, y)\rightarrow (\sqrt{\tau}\cos\theta, \sqrt{\tau}\sin{\theta})$, then Jacobian is
$$ |\frac{\partial(x, y)}{\partial(\tau, \theta)}| = \det\left(\begin{bmatrix}\frac{\cos\theta}{2\sqrt{\tau}} & -\sqrt{\tau}\sin\theta\\\frac{\sin\theta}{2\sqrt{\tau}} & \sqrt{\tau}\cos\theta \end{bmatrix}\right) = \frac{1}{2}, $$ So, $$ \begin{align} S_1 &= \int_{\{(x, y)|x^2+y^2=t\}}dxdy \\ &=\int_0^{2\pi}\int_0^\infty\delta(\tau-t)|\frac{\partial(x, y)}{\partial(\tau, \theta)}|d\tau d\theta \\ &=\frac{1}{2}\int_0^{2\pi}d\theta\int_0^\infty\delta(\tau-t)d\tau\\ &= \pi. \end{align} $$ Why is the second result $S_1\ne 2\pi\sqrt{t}\ ?$
When using the Dirac delta-function, you have to be careful. In particular, if you have a function $f(x)$ with a single zero at $x_0$, you might think that you can just replace $$ \int g(x) \delta(f(x)) \, dx \stackrel{?}{=} \int g(x) \delta(x - x_0) \, dx. $$ But this is not the case; the proper rule is $$ \int g(x) \delta(f(x)) \, dx = \int g(x) \frac{\delta(x - x_0)}{|f'(x_0)|} \, dx. $$ This can be proven by restricting the integration to a small region about $x = x_0$, defining an inverse function $x(f)$ in that region (note that this requires $f'(x_0) \neq 0$), and changing integration variables from $x$ to $f$. The factor of the derivative in the denominator comes from this variable change.
Thus, if we change variables from $\rho$ to $\tau$, we can view $\rho = f(\tau) = \sqrt{\tau}$ and $f'(\tau) = \frac{1}{2} \tau^{-1/2}$; so we have $$ \delta(\rho - r) = \frac{\delta(\tau - r^2)}{\frac{1}{2} \tau^{-1/2}} = 2 \sqrt{\tau} \delta(\tau - r^2). $$ When integrated over all $\tau$ and $\theta$, this will yield the result $2 \pi r$ as expected.
More generally, when doing an multi-variable integral over a single delta-function that is a function of several variables, we have the identity $$ \int_{\mathbb{R}^n} g(\mathbf{x}) \delta(f(\mathbf{x})) \,d^n \mathbf{x} = \int_{f^{-1}(0)} \frac{g(\mathbf{x})}{|\nabla f|} da, $$ where $f^{-1}(0)$ is the $(n-1)$-dimensional surface on which $f$ vanishes and $da$ is the area element on this surface. In your case, you when you transferred from $(x,y)$ to $(\tau,\theta)$ coordinates, you had implicitly $f(x,y) = x^2 + y^2 - t$; thus, there should have been a factor of $|2 \sqrt{x^2 + y^2}|$ in the denominator of your integral over $\tau$ and $\theta$.