Consider the following two integrals:
\begin{align} I_1&:=\iint_{0\leq x,\, y\leq 1,\,x=y}\,dx\, dy,\\ I_2&:=\int_0^1\int_0^1 \delta(x-y)\,dx \,dy. \end{align}
I believe $I_1=0$ because it is the measure "surface area" of the line segment $y=x$ inside a two dimensional rectangle.
From the properties of the delta function, we have
$$I_2=\int_0^1\int_0^1 \delta(x-y)\, dx\, dy=\int_0^1 1 \,dy=1$$
as $$\int_0^1 \delta(x-y)\, dx =1$$ for any $y\in ]0,1[.$
Question: I understand very well that the integrands are different in $I_1$ and $I_2$. In $I_1$, the integrand is the characteristic function $\chi_{0\leq x,y\leq 1,x=y}(x,y)$. But these integrands are different on a negligible set... I am confused!
There's no need to make this about double integrals. Your question would apply just as well to the two integrals $$\int_{[0,1]}\chi_{\{0\}}(x)\,dx=0$$ and $$\int_{[0,1]}\delta(x)\,dx=1,$$ which are different even though $\delta$ and $\chi_{\{0\}}$ are the same everywhere except at $0$. The explanation for this is simply that $\delta$ is not actually a function and its "integral" is not actually an integral in the usual sense. It's a theorem that if two integrable functions agree except on a set of measure zero, then their integrals are the same. But $\delta$ is not a function at all, so that theorem doesn't apply. Instead, $\int_{[0,1]}\delta(x)\,dx$, though written as an "integral", is actually simply defined as the evaluation of the distribution $\delta$ on the constant function $1$. Since $\delta$ is by definition the distribution that takes a function and outputs that function evaluated at $0$, this gives the answer of $1$.
Intuitively, $\delta$ is a "function" that gets infinitely large at the point $0$, such that even though it's just one single point, the "area" under it somehow manages to be positive. But it's not actually literally the function $\mathbb{R}\to[-\infty,\infty]$ that is $\infty$ at $0$ and $0$ elsewhere; that is an actual genuine integrable function whose integral is $0$.