Why is it that $\int_a^b \int_c^d f(x)g(y)\,dy\,dx=\int_a^b f(x)\,dx \int_c^d g(y)\,dy$?

Since each $x$-slice has signed-area $f(x)\int_c^d g(y)dy$ we can just say $\int_c^d g(y)dy=I$ and so the area of the $x$-slice of the volume in question is just $I\cdot f(x)$. To find the volume we simply add up the volume of the infinitesimal slices $I \cdot f(x) dx$ from $x=a$ to $x=b$ thus $$\iint_{R} f(x)g(y) \ dA = \int_a^b I \cdot f(x) \, dx = I\int_a^b f(x) \, dx = \int_c^d g(y) \, dy \int_a^b f(x) \, dx. $$ I'd like to give a more geometric answer, I suspect there is something more to say here. Why, geometrically, is the volume the product of these particular areas? Indeed, now that I say that, it is strange.

Consider a geometric interpretation of the much simpler expression, $f(a)g(c)$. The usual interpretation is that we have line segments of lengths $f(a)$ and $g(c)$, respectively, and we construct a rectangle using copies of those two segments as adjacent sides of the rectangle.

The rectangle could be in any orientation, but when representing a set of products $f(x)g(y)$ for various values of $x$ and $y$, we can arbitrarily say we'll orient the rectangles so that all the $f(x)$ sides are parallel to those of the other rectangles and all the $g(y)$ sides are parallel.

Now let's consider the geometric interpretation of the integral $\int_a^b f(x) \,dx$. Often we'll think of this function graphed on the $x,y$ plane, that is, we'll set $y = f(x)$, but as we're going to be combining $f(x)$ with a function of $y$ later, the $x,y$-plane interpretation is a dead end, so let's instead set $z = f(x)$ and graph the function in the $x,z$-plane. Then the integral can be interpreted by considering all the disjoint planar regions contained between the curve $z = f(x)$ and the $x$-axis (the curve $x = 0$): add together all such regions that are "above" the $x$-axis, and from this result subtract all such regions that are "below" the $x$-axis.

A similar geometric interpretation can be made of $\int_c^d g(y) \,dy$, except of course now the regions to add or subtract are "above" or "below" the $y$-axis rather than the $x$-axis.

But we have to avoid falling into the trap of visualizing a graph of $z = g(y)$ in the $y,z$-plane, because that would make all the $g(y)$-length segments parallel to all the $f(x)$-length segments in $x,y,z$-space, whereas when visualizing $f(x)g(y)$ we wanted the $g(y)$-length segment to be perpendicular to the $f(x)$-length segment. So we really want to graph $g(y)$ in a direction that is not only perpendicular to $y$, and not in the $x,y$-plane; in fact, we want it perpendicular to $x$ and to $z$ as well.

In other words, we are forced into a fourth dimension. Let's say we're now in $x,y,z,w$-hyperspace, so we can visualize $\int_c^d g(y) \,dy$ by plotting $w = g(y)$ in the $w,y$-plane.

Then for each $x \in [a,b]$ and for each $y \in [c,d]$, we have a segment of length $f(x)$ parallel to the $z$-axis, with one end at $(x,0,0,0)$ and the other at $(x,0,f(x),0)$; we also have a segment of length $g(y)$ parallel to the $w$-axis, with one end at $(0,y,0,0)$ and the other at $(0,y,0,g(y))$. We project both these segments parallel to the $x,y$-plane onto segments that each have one end at the point $(x,y,0,0)$. These projected segments are the sides of a rectangle parallel to the $w,z$-plane, and the area of that rectangle is $f(x)g(y)$.

Now we take the four-dimensional object composed of all such rectangles for all values of $(x,y) \in [a,b]\times[c,d]$, and we interpret the integral $\int_a^b \int_c^d f(x) g(y) \,dy\, dx$ as the net four-dimensional volume obtained by adding the hypervolumes of all parts of this four-dimensional object that project onto the first or third quadrants of the $w,z$-plane (that is, $wz > 0$) and subtracting the hypervolumes of all the other parts. In other words, we consider all points with $wz > 0$ to be "above" the $x,y$-plane, just as when interpreting $\int_a^b f(x) \,dx$ in the $x,z$-plane we consider all points with $z > 0$ to be "above" the $x$-axis.

A more formal definition of the four-dimensional object as a subset of $\mathbb{R}$ has already been offered, so I won't repeat that definition, but the interpretation of such an object in the case where $f$ and $g$ are not both non-negative functions should now be clear.

Consider the subset of $\Bbb R^4$ defined by $X = \{ (x,y,u,v) \mid x \in [a;b], y \in [c;d], u \in [0 ; f(x)], v \in [0 ; g (y)]\}$.
This looks like a family of rectangles $[0 ; f(x)]\times[0;g(y)]$ indexed by two coordinates $x,y$.

To get the measure of this $4$-dimensional shape, you can do some standard integration and find $\int_x \int_y \int_u \int_v 1_X(x,y,u,v)\; dv \,du\,dy\,dx = \int_a^b \int_c^d \int_0^{f(x)} \int_0^{g(y)} 1\; dv \,du\,dy\,dx = \int_a^b \int_c^d f(x)g(y) \;dy \,dx$.

You can also notice that $X$ is the cartesian product of $Y = \{(x,u) \mid x \in [a;b], u \in [0 ; f(x)]\}$ with $Z = \{(y,v) \mid y \in [c;d], v \in [0 ; g(y)]\} $

So the measure of $X$ is the product of the measures of those two sets, which you can easily see to be $(\int_a^b f(x)dx) (\int_c^d g(y) dy)$ ($Y$ and $Z$ are quite literally the areas under the curves $f$ and $g$)