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$? The title sums it up. It's simple to prove, but I'm wondering if there is a geometric interpretation?
 A: 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. 
A: 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.
A: 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$)
A: This is a really lovely question, precisely because when you first look at it you think "Oh, that's trivial" -- and then you realize that no, wait a minute, there's something weird here after all.  (See also the last sentence of James Cook's response and the comments on it.)
To get a grip on this, let's consider a very simple case:  Suppose $f$ and $g$ are both constant functions; in fact let's be really specific and say $f(x)=3$ and $g(x)=2$.  (If we can understand what's happening for the case of constant functions, we will probably be able to extend that interpretation to the general case by thinking of Riemann integrals as lots of little infinitesimal "towers" on the $xy$-plane.)
Furthermore, let's say $[a,b]=[0,5]$ and $[c,d]=[0,7]$.
Now, we can interpret $\int_a^b f(x)  \, \mathrm{d}x$ as a 5 x 3 rectangle (area 15), $\int_c^d g(y)  \, \mathrm{d}y$ as a 7 x 2 rectangle (area 14), and $\int_a^b \int_c^d f(x)g(y) \, \mathrm{d}y \, \mathrm{d}x$ as a rectangular solid with a 5 x 7 base and a height of 6 (volume = 210).  So your question is:
Can we somehow "see" the 5 x 7 x 6 rectangular solid as the product of the 5 x 3 rectangle and the 7 x 2 rectangle?
Having gone through all of this, I think the answer (surprise!) is No, at least not naturally.  And here's why:
You may have noticed that something weird is going on with our dimensions -- or our units, if you prefer to think of them that way.  Suppose $x,y,f(x)$ and $g(y)$ all have dimensions of length.  In that case $\int_a^b f(x)  \, \mathrm{d}x$ and $\int_c^d g(y)  \, \mathrm{d}y$ each have dimensions of area (length$^2$), and $\int_a^b \int_c^d f(x)g(y) \, \mathrm{d}y \, \mathrm{d}x$ ought to have dimensions of length$^4$... but for some reason we are visualizing it as a volume.
And that, right there, points us at what is wrong with our setup, and the reason why we can't geometrically interpret the double-integral as a product of two single-integrals:  The function $f(x)g(y)$ has dimensions of area, but when we interpret it as the height of a surface above the xy-plane we are plotting it as if it were a length.  If we could visualize four dimensions, we would see $\int_a^b \int_c^d f(x)g(y) \, \mathrm{d}y \, \mathrm{d}x$ like this:  To each point in the 5 x 7 rectangle over which we are integrating, there would be attached not a 6-unit tall vertical line segment, but rather a 2 x 3 rectangular area (both dimensions orthogonal to the $xy$-plane).  So the whole thing would look like a 5 x 7 x 2 x 3 four-volume, which could easily be recognized as the product of the two 5 x 3 and 7 x 2 rectangles.  But because we can't visualize four dimensions, we replace the 2 x 3 rectangular area with a 6-unit tall line segment.
In doing so we "smoosh" a 4-dimensional hypersolid down to a 3-dimensional solid.  Collapsing two dimensions down into one mixes together two independent dimensions so that we can no longer clearly see the relationship between the two rectangular factors and their four-dimensional product.
Now go read David K's answer, which I think is a good complement to mine.
A: For the geometric interpretation of the formula in the title of the question assume $$f(x)\geq0\quad(a\leq x\leq b),\qquad g(y)\geq0\quad(c\leq y\leq d)\ .$$
Then the integrals $\int_a^b f(x)\ dx$ and $\int_c^d g(y)$ can be viewed as areas of the two shapes $$A:=\{(x,u)\ |\ a\leq x\leq b, \ 0\leq u\leq f(x)\},\quad B:=\{(y,v)\ |\ c\leq x\leq d, \ 0\leq v\leq g(y)\}\ .$$ The four dimensional volume of the cartesian product $X:=A\times B$ is equal to
$${\rm vol}_4(X)={\rm vol}_2(A)\cdot{\rm vol}_2(B)\ ,$$
much as the area of a rectangle is width times height. On the other hand $X$ can be viewed as attaching to each point $(x,y)\in[a,b]\times[c,d]$ a rectangle of width $f(x)$ and height $g(y)$ pointing out into four-dimensional space. The volume of $X$ is then computed as
$${\rm vol}_4(X)=\int_{[a,b]\times[c,d]} f(x)\>g(y)\ {\rm d}(x,y)\ ,$$
much as the area of $A$ is computed as $\int_a^b f(x)\ dx$.
When $f$ and $g$ can take negative values as well the idea of "signed" areas or volumes is not very helpful. Its more convenient to argue in an algebraical way, making use of the "totally Cartesian product situation" here. It is then seen that, apart from convergence questions, there is just the distributive law  at work.
Consider a partition of $[a,b]$ into subintervals $P_j$ of length $\lambda(P_j)$, and similarly a partition of $[c,d]$ into subintervals $Q_k$ of length $\lambda(Q_k)$. Then the rectangles $R_{j,k}:=P_j\times Q_k$ have area $\mu(R_{j,k})=\lambda(P_j)\lambda(Q_k)$, and together they form a partition of $[a,b]\times[c,d]$.
Choose sample points $\xi_j\in P_j$ and $\eta_k\in Q_k$. Putting $fg=:h$ we then have the exact equality
$$\sum_{j,\>k}h(\xi_j,\eta_k)\mu(R_{j,k})=\sum_{j,\>k}f(\xi_j)g(\eta_k)\>\lambda(P_j)\lambda(Q_k)=\sum_jf(\xi_j)\lambda(P_j)\ \sum_k g(\eta_k)\lambda(Q_k)\ .$$
The sums appearing on the left and right hand sides of this equation can be seen as Riemann sums for the integrals in the title of the question, and as the two sides  are always equal their limits have to be equal, too.
A: Any geometric interpretation of this quantity must be quite nuanced.  James S. Cook has suggested this is a "signed volume".  Let's take that idea to a logical conclusion.
If $D = \ell \times m$ is the 2d domain of integration, composed of a rectangle of two intervals $\ell$ and $m$, and we take $V$ as a "signed volume" such that
$$V = \int_D f(x) g(y) \, dx \, dy = \left[\int_\ell f(x) \,d x \right] \times \left[ \int_m g(y) \, dy \right]$$
then what are the "units" of $f, g$?
In other words, if we are to interpret $V$ as a volume, then we associate it with the dimensions of length cubed.  Then the most symmetric way to break down the 1d integrals is to associate them with length to the three-halves, not with areas (length squared).

One should remember that the interpretation of a 1d integral as a signed area is getting the punchline backwards:  the integral describes this signed area quantity, yes, but that does not mean an integral is this signed area.  Integrals are just continuum versions of weighted sums, and as such, the picture of the integral here is quite clear:  the domain can be cut into horizontal strips of constant $y = y_0$, and the weight function on those strips is just $f(x) g(y_0)$.  At that value of $y_0$, integration over $x$ collapses the strip to a point at $y=y_0$ with weight $g(y_0) \int_\ell f(x) \, dx$.  The resulting line with its associated weight function can be integrated over again to get the final answer.
A: Here is my heretic attempt to give a geometrical explanation. First, have a look at the following diagrams, where we have an illustration of
$$
\int_a^b  \int_c^d x^2 cos(y) dx dy = \int_a^b x^2 dx \int_c^d cos(y) dy 
$$
for a, b = 1, 2 and c, d = 0, 1.4. In the diagram below I draw a rectangle with sides $ \int_a^b x^2 dx $ and $\int_c^d cos(y) dy$. So, the integrals are in this drawing a length instead of an area! This is possible because areas are real numbers and as such they can be seen as numbers (lengths) on the real line. For a complete heresy you can think of $f(x)$ and $dx$ units as being $\sqrt{length}$.
With this in mind, each segment in bottom side, say, is $a_i = f(x_i) \frac{1}{n}$, where we are attempting to approximate the integral in $\frac{1}{n}$ steps. Likewise along the vertical sides of the figure. There, a single square (in blue) is $a_i b_j$ and we have, by distributive law as mentioned by someone:
$$
\sum_{i=1}^{n} \sum_{j=1}^{m} a_i b_j =_1 \sum_{i=1}^{n} a_i (\sum_{j=1}^{m} b_j) =_2 (\sum_{i=1}^{n} a_i) (\sum_{j=1}^{m} b_j)   
$$
In $=_1$ we used that $a_i$ does not vary in $\sum_{j=1}^{n}$, and in $=_2$ we used the summation $\sum_{j=1}^{n} b_j$ is fixed for all $i$. The sum of the blue rectangles is the same as the product of the sides. To complete the argument, $a_i = (x_i)^2 \frac{1}{n}$ and $b_j = cos(x_j)^2 \frac{1}{m}$ are the summation terms in simple partitions of the corresponding integrals.

Now the final, interesting part. Let's increase the number of segments in each partition. Here a picture:

As we increase the number of partitions, the sides become closer to the real value, i.e. the length (heresy!) should be around the value of the corresponding integral and the geometric argument above still applies.
Is this geometric enough? ;-)
Edit: In this geometrical view, a negative "area" (length) can be seen as the product of two directed $\sqrt{length}$ measurements. I think one should be able to picture some intervals going from right to left when the corresponding partial sum (or term) is negative by using appropriate examples, but the image eludes me for now.
A: Let us try this interpretation, which is really rephrasing what others have already pointed out (particularly caya's answer). 
The geometric interpretation fof,
$$\int_{A}^{B} dx \int_{C}^{D} dy = (B-A)(D-C)$$
is clear. However, we have $f(x)\,g(y)$ in the integrand, which is not necessarily semi-positive or semi-negative definite functions over the integration domain. Now, suppose that $f(x)$ and $g(y)$ are "nice enough" functions so that we can break up the respective intervals $[a,b]$ and $[c,d]$ into smaller intervals $[a_i,b_i]$ and $[c_i,d_i]$ where the functions are semi-positive or semi-negative definite. 
Let's focus on a pair of such intervals, say, $[a_\ell,b_\ell]$ and $[c_j,d_j]$, where, for example, $f(x)$ is positive definite (except at $a_\ell$ and $b_\ell$) and $g(y)$ is negative-definite (except at $c_j$ and $d_j$). Then we can think of these integrands as Jacobians, i.e.,
$$ \int_{a_\ell}^{b_\ell} f(x)\,dx \,=\, \int_{A_\ell}^{B_\ell} d\tilde{x}$$
and, in this example where $g(y)$ is negative definite over $(c_j,d_j)$,
$$ \int_{c_j}^{d_j} g(y)\,dy \,=\, -\int_{C_j}^{D_j} d\tilde{y} $$
A: Not a geometric proof, but a hint: product of sums, distributivity.
$$(\sum_ia_i).(\sum_jb_j)=\sum_i\sum_j(a_i.b_j)$$
A: I am not good in plotting, but consider function $z=F(x,y)=f(x)*g(y)$ and plot it in 3 dimensional space. The volume this surface "covers" is exactly the integral we consider ( we allow negative volume too). The fact that $F(x,y)=f(x)*f(y)$ means that any section of this graph by $x=const$ or $y=costnt$ plane will get your functions $f(x),g(y)$ up to multiplier.
