We can treat one variable delta function as
$$\delta(f(x)) = \sum_i\frac{1}{|\frac{df}{dx}|_{x=x_i}} \delta(x-x_i).$$
Then how do we treat two variable delta function, such as $\delta(f(x,y))$?
for example, how calculate $\int\int \delta(x-y)$ ?
I first thought using $\int f(x)\delta(x-y)dx = f(y)$
$$\int\int \delta(x-y)dxdy = \int\int1 * \delta(x-y)dx dy = \int dy.$$
but this is nonsense, since we can also think it as $\int\int \delta(x-y)dxdy = \int dx$
On
$$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{d}x\,\mathrm{d}y\, f(x)g(y)\delta(x,y) = \int_{-\infty}^{+\infty}\mathrm{d}y\, f(y)g(y),$$ or $$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{d}x\,\mathrm{d}y\, f(x)g(y)\delta(x,y) = \int_{-\infty}^{+\infty}\mathrm{d}x\, f(x)g(x),$$ which is the same since the integration limits are the same, and the integration variable is just a dummy variable, whose actual name doesn't matter. Now what if the limits are not the same? For example: $$\int_{a}^{b}\mathrm{d}x\,\int_{c}^{d} \mathrm{d}y\, f(x)g(y)\delta(x,y).$$ Then we can use the Heaviside step function to write: $$h(x)=f(x)(\theta(x-a)-\theta(x-b)),$$ and $$l(y)=g(y)(\theta(y-c)-\theta(y-d)),$$ and so: $$\int_{a}^{b}\mathrm{d}x\,\int_{c}^{d} \mathrm{d}y\, f(x)g(y)\delta(x,y) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{d}x\,\mathrm{d}y\, h(x)l(y)\delta(x,y),$$ and we can use the previous result. Note that I have used one dimensional Cartesian variables. When dealing with non-Cartesian coordinates (for example, spherical coordinates), the Dirac delta $\delta(x,y)$ can (depending on your notation), become a scalar in the first variable $x$ and a scalar density in the second variable $y$, and you may need to insert a factor of the square root of the determinant of the metric to make everything work. If you need more details let me know - coincidentally a couple of people had asked me about this recently, so I have a lot of it LaTeXed, just not handy at the moment.
The simplest way to start is to think about how the delta function is defined in 1D, $$\int_{-\infty}^{\infty} f(x)\delta(x-a)dx = f(a)$$ which you correctly identified above. This is saying that some function $f(x)$ multiplied by the a shifted delta function $\delta(x-a)$ and integrated everywhere returns the same function evaluated at the value given by the shift, $f(a)$. We want the same sort of definition to apply when we integrate everywhere in 2D, so we have to have a 2-variable function $f(x,y)$ and it has to have a 2-variable delta function $\delta(x,y)$. (Notice that $\delta(x-y)\neq \delta(x,y)$ in general.) Applying it to the 2D case in this way we should give $$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)\delta(x-a,y-b) \ dx dy = f(a,b)$$
which should be true whether you integrate either $x$ or $y$ first as shown below.
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)\delta(x-a,y-b) \ dx dy = \int_{-\infty}^{\infty} f(a,y)\delta(y-b) \ dy = f(a,b)$$
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)\delta(x-a,y-b) \ dx dy = \int_{-\infty}^{\infty} f(x,b)\delta(x-a) \ dx = f(a,b)$$
These are easiest to solve when $f(x,y)$ is a separable function like $f(x,y)=g(x) h(y)$, but some harder functions are still quite doable.