Explain $\iint \mathrm dx\,\mathrm dy = \iint r \,\mathrm \,d\alpha\,\mathrm dr$

Generally speaking, a double integral always has an area differential, so that you're integrating $\iint dA$. Another way to view the question, then, is why $dA = dx dy$ in Cartesian coordinates and $dA = r dr d\theta$ in polar coordinates.

An area element in Cartesian is a rectangle, as Qiaochu describes in his answer. The area of the rectangle is the small change in $x$ times the small change in $y$, or $\Delta x \Delta y$.

An area element in polar, however, is a piece of a circle sector. There's a nice picture below taken from here. (The area element is the shaded part.)

If the angle is measured in radians, we know that the area of a sector of angle $\theta$ of a circle of radius $r$ is $\frac{1}{2}r^2 \theta$. So the area of the shaded piece in the picture is $$\Delta A = \frac{1}{2}(r + \Delta r)^2 \Delta \theta - \frac{1}{2}r^2 \Delta \theta = r \, \Delta r \, \Delta \theta + \frac{1}{2}(\Delta r)^2 \, \Delta \theta.$$ The quadratic factor of $\Delta r$ makes the second term negligible compared to the first term for small enough $\Delta r$ and $\Delta \theta$. Thus in the limit we get $dA = r \, dr \, d\theta$.

As others have said, you can also use the multivariate change of variables formula involving the Jacobian of the transformation directly. I like the geometric argument when first introducing the polar element in a calculus course, though.

It's a special case of the multivariate change of variables formula. Intuitively you can think of it as follows: starting from a point $(x, y)$, you make an infinitesimal change in $x$ and then an infinitesimal change in $y$ to get to $(x + \delta x, y + \delta y)$. Those changes trace out a little square whose area is $\delta x \delta y$, and you're summing over a bunch of these little squares.

So what happens when you do the same thing in polar coordinates? Starting from $(r, \theta)$ you move to $(r + \delta r, \theta + \delta \theta)$. Now moving $\delta r$ is just like moving $\delta x$ in an appropriately rotated set of axes. But when you move $\delta \theta$, the actual distance you move by is multiplied by $r$ (draw a diagram to see this), and it's in a direction orthogonal to the direction you move when changing $r$. You're actually moving by $r \, \delta \theta$. The result is not quite a square, but for small enough values of $\delta r, \delta \theta$ it approaches a square with side lengths $\delta r$ and $r \, \delta \theta$.

Let us approach this problem like some physicists would, that is, let us consider the differential $\mathrm{d}$ as an operator which obeys the rules of derivation, only with a sign.

In the present case, $x=r\cos\alpha$ and $y=r\sin\alpha$ hence $$ \mathrm{d}x=\cos\alpha\mathrm{d}r-r\sin\alpha\mathrm{d}\alpha, \quad \mathrm{d}y=\sin\alpha\mathrm{d}r+r\cos\alpha\mathrm{d}\alpha. $$ Now, there are some magic rules which allow to multiply $\mathrm{d}r$ and $\mathrm{d}\alpha$ elements. These are $$ \mathrm{d}r\mathrm{d}r=0,\quad \mathrm{d}\alpha\mathrm{d}r=-\mathrm{d}r\mathrm{d}\alpha,\quad \mathrm{d}\alpha\mathrm{d}\alpha=0. $$ Hence, $$ \mathrm{d}x\mathrm{d}y=(\cos\alpha\mathrm{d}r-r\sin\alpha\mathrm{d}\alpha)(\sin\alpha\mathrm{d}r+r\cos\alpha\mathrm{d}\alpha). $$ The $\mathrm{d}r\mathrm{d}r$ and $\mathrm{d}\alpha\mathrm{d}\alpha$ terms disappear and the terms which interest us are the $\mathrm{d}r\mathrm{d}\alpha$ and $\mathrm{d}\alpha\mathrm{d}r$ ones. One gets $$ \mathrm{d}x\mathrm{d}y=(\cos\alpha\cdot r\cos\alpha-r\sin\alpha\cdot (-1)\sin\alpha)\mathrm{d}r\mathrm{d}\alpha, $$ that is, $$ \color{green}{\mathrm{d}x\mathrm{d}y=r\mathrm{d}r\mathrm{d}\alpha}, $$ a formula which yields $$ \color{red}{\iint f(x,y)\mathrm{d}x\mathrm{d}y=\iint f(r\cos\alpha,r\sin\alpha)r\mathrm{d}r\mathrm{d}\alpha}. $$ One can also transform integrals the other way round, all there is to do is to compute a formula for $\mathrm{d}r\mathrm{d}\alpha$ as a multiple of $\mathrm{d}x\mathrm{d}y$. Our formula for $\mathrm{d}x\mathrm{d}y$ in terms of $\mathrm{d}r\mathrm{d}\alpha$ yields $$ \mathrm{d}r\mathrm{d}\alpha=\frac1r\mathrm{d}x\mathrm{d}y=\frac1{\sqrt{x^2+y^2}}\mathrm{d}x\mathrm{d}y, $$ hence $$ \color{blue}{\iint f(r,\alpha)\mathrm{d}r\mathrm{d}\alpha=\iint f(x,y)\frac1{\sqrt{x^2+y^2}}\mathrm{d}x\mathrm{d}y}. $$ Once again, this only describes the computational side of the story, but this recipe is supported by a well established theory of differential forms which we omitted.