In my opinion, you never get a True in this way. The delta function is the kernel of the identity map by an integral

f[x] == Integrate[ f[y] delta[x.y]], {x-inf,inf}},{y-inf,inf}, {z,-inf, inf]} 

I am not wrting delta(x-y) but delta(x,y) because one cannnot add points in curvilinear coordinates.

Now, transform the Integral to spherical coordinates.

f[r, theta, phi ] ==
 Integrate[  r'^2  Sin[theta'] * f[ r',theta', phi' ] * delta[{r ,theta, phi},{ r' theta', phi'} ], {r', 0, inf}, {theta',0, Pi},  {phi', 0, 2 Pi]}

Next replace the the point delta by the product of coordinate deltas divided by the volume element

Integrate[  r'^2  Sin[theta']  * f[ r', theta', phi' ] * delta[r-r']  delta[ theta-theta'] *delta[phi-phi']\(r'2 Sin[theta'], {r', 0, inf}, {theta',0, P},  {phi', 0, 2 Pi}] 

The volume density in the denominator cancels the volume density under the integral and the two-point delta is replaced by the product of ordinary deltas of coordinate differences and the integral is now a threefold coordinate integral with volume density 1.