Evaluating (Uniform) Expectations over Non-simple Region 
Background. Let $V = (X,Y)$ be a random vector in 2-dimensions uniformly distributed over two disjoint regions $R_X \cup R_Y$ defined as follows:
$$
\begin{align}
R_X &= ([0,1] \times [0,1]) \setminus \left(\bigcup A\right) \\
R_Y &= [1,2] \times [1,2],
\end{align}
$$
where
$$
A = \{[.2,.4] \times [.2,.4], [.2,.4] \times [.6,.8], [.6,.8] \times [.2,.4], [.6,.8] \times [.6,.8]\}.
$$
A plot of the regions (green is the density/area associated with each region) is:


I'm interested in how to find $E(XY)$. So, letting $\lambda$ be the Lebesgue measure the associated pdf is
$$
f_{XY}(x,y) = 
\begin{cases}
\frac{1}{\lambda\left(R_X \cup R_Y\right)} = \frac{1}{\lambda([0,1] \times [0,1]) - \lambda\left(\bigcup A\right) + \lambda([1,2] \times [1,2])} \approx \frac{25}{46}, &(x,y) \in R_X\cup R_Y \\
0, &\text{otherwise}
\end{cases}
$$
Using the traditional definition
$$E(XY) = \int_{R_X\cup R_Y} xyf_{XY}(x,y)d\lambda = \frac{25}{46}\int_{R_X\cup R_Y} xyd\lambda = \frac{25}{46}\left(\int_{R_X} xydxdy + \int_{R_Y} xydxdy\right).$$
Integrating over $R_Y$ is straightforward. But for the "non-simple" region $R_X$ would we calculate it as
$$\int_{R_X} xydxdy = \int_0^1\int_0^1 xydxdy - \sum_{a \in A} \int_a xydxdy~\text{?} \tag{1}$$
What if there were countably many boxes to remove from $R_X$ instead of 4 finite ones? Does the formula for (1) generalize (I'm assuming here we could use something like the MCT/DCT/etc. to evaluate the sum).

EDIT (after @whuber answer): Quick follow-up points:

*

*It's apparent now that $E(X)$ could be found using this mixture approach too, i.e.

$$\sum_i \omega_i p_i = \frac{21}{25}$$
and
$$E(X) = \frac{25}{21}\left(\frac{1}{2} - \frac{9}{2500} - \frac{21}{2500} - \frac{21}{2500} - \frac{49}{2500}\right) = \frac{23}{42}.$$


*For the case where instead of 4 finite regions removed, we have countably many, then we can generalize as

$$E(XY) = \frac{25}{46}\left(\iint_{[0,1]\times [0,1]} xydxdy + \iint_{[1,2]\times [1,2]} xydxdy - \sum_{j}^{\infty}g(x)\mathscr{I}_j(x)\right) = \frac{25}{46}\left(\frac{1}{4} + \frac{9}{4} - \sum_{j}^{\infty}g(x)\mathscr{I}_j(x)\right),$$
for $j > i$ and assuming that $\sum_{j}^{\infty}g(x)\mathscr{I}_{j}(x) < \infty$ (i.e. converges).
 A: Let's solve the more general problem and then apply it to the specific setting as an illustration.
Suppose $g:\mathcal{X}\to\mathbb{R}$ is a measurable function.  Let $\{\mathcal{A}_i\}\subset \mathscr{P}(\mathcal{X})$ be a finite or countable collection of subsets of $\mathcal X,$ each with finite positive measure $p_i = \int_{\mathcal{A}_i}\mathrm{d}x.$  Associated with each $\mathcal A_i$ is its indicator function $\mathscr{I}_i.$  Then for any sequence of numbers $(\omega_i)$ -- essentially by construction of the integral -- we have that
$$\int_\mathcal{X} \sum_i g(x)\omega_i\mathscr{I}_i(x)\,\mathrm{d}x = \sum_i \omega_i \int_{\mathcal{A}_i} g(x)\,\mathrm{d}x =  \sum_i \omega_i p_i \int_{\mathcal{A}_i} g(x)\,\frac{\mathrm{d}x}{|p_i|}$$
When the $\omega_i$ are positive and the sum of $\omega_i p_i$ is unity, the right hand side is the expectation of $g(X)$ where the distribution of $X$ is a mixture of the uniform distributions on the $\mathcal{A}_i$ with mixture weights $\omega_i p_i.$  I will continue to use this imagery and this language even when some of the $\omega_i$ are negative.  Think of this as a "generalized mixture" if you like.
Provided the left hand integrand is never negative and is positive on some set of positive area, we may normalize it to produce a genuine distribution.  Evidently its density function is

$$f(x;(\mathcal{A}_i), (\omega_i)) = \frac{1}{\sum_i \omega_i p_i} \sum_i \omega_i g(x) \mathscr{I}_i(x).\tag{*}$$

The first formula therefore gives the expectation $E[g(Z)]$ when $Z$ has this mixture distribution.  The right hand side of the formula shows that this expectation is a linear combination of the expectations of the mixture components.

To apply this observation to the example in the question, let $\mathcal X =\mathbb{R}^2$ (with its usual Borel measure) and $g(x,y) = xy.$  Let's begin by getting all the calculations out of the way.  They amount to integrating $g$ over various rectangles $[a,b]\times[c,d].$  It is elementary to compute that

*

*The area of a rectangle $[a,b]\times[c,d]$ is $p(a,b,c,d)=(b-a)(d-c).$


*$$G(a,b,c,d)=\iint_{[a,b]\times[c,d]}xy\,\mathrm{d}x\mathrm{d}y = \frac{1}{4}(b^2-a^2)(d^2-c^2) = p(a,b,c,d)(a+b)(c+d)/4.$$
The problem can be expressed in terms of six rectangles: the two big ones (which therefore receive weights $\omega_i=1$) from which the four little ones have been removed (by applying weights $\omega_i=-1$).  Here is a table of their properties, computed using $(1)$ and $(2)$ above.
$$
\begin{array}[llrrrl]
    & i &  [a,b]& [c,d] & p & G & \omega \\ \hline
    1 &  [0,1]& [0,1] & 1 & 1/4 & 1 \\
    2 &  [1,2]& [1,2] & 1 & 9/4 & 1 \\
    3 &  [1/5,2/5] & [1/5,2/5] & 1/25 & 9/100 & -1 \\
    4 &  [3/5,4/5] & [1/5,2/5] & 1/25 & 21/100 & -1\\
    5 &  [1/5,2/5] & [3/5,4/5] & 1/25 & 21/100 & -1\\
    6 &  [3/5,4/5] & [3/5,4/5] & 1/25 & 49/100 & -1
\end{array}
$$
The denominator in $(*)$ is
$$\sum_{i=1}^6 \omega_i p_i = 1 + 1 - \frac{1}{25} - \cdots - \frac{1}{25} = \frac{46}{25}.$$
Writing $Z=(X,Y),$ the answer to the question is
$$E[XY] = E[g(Z)] = \frac{25}{46}\left(\frac{1}{4} + \frac{9}{4} - \frac{9}{2500} - \frac{21}{2500} - \frac{21}{2500} - \frac{49}{2500}\right)=\frac{123}{92}.$$
