Here is the proof that adapts Glasser's argument for 1-D case:
Theorem. Let $a, b, c, d$ be real numbers such that $a, d > 0$ and $ad-bc> 0$. Also, define
$$ T_1(x) = x - \sum_{k=1}^{K} \frac{p_k}{x - a_k}, \qquad T_2(y) = y - \sum_{l=1}^{L} \frac{q_l}{y - b_l}, $$
where $p_k$'s and $q_l$'s are positive and $a_k$'s and $b_l$'s are real numbers. Then for any non-negative measurable function $f(x, y)$,
$$ \int_{\mathbb{R}^2} f(a T_1(x) + by, cx + dT_2(y)) \, \mathrm{d}x\mathrm{d}y = \frac{1}{ad-bc} \int_{\mathbb{R}^2} f(x, y) \, \mathrm{d}x\mathrm{d}y. $$
Proof. If $b = 0$ or $c = 0$, the claim immediately follows from the Glasser's Master Theorem. Therefore, we may assume without loss of generality that $b, c \neq 0$.
Defining the Inverses
Let $U_k = (a_k, a_{k+1})$ for $k \in \{0, \dots, K\}$ and $V_l = (b_l, b_{l+1})$ for $l \in \{0, \dots, L\}$, using the convention $a_0 = b_0 = -\infty$ and $a_{K+1} = b_{L+1} = \infty$. Each restriction $T_1|_{U_k}$ possesses an inverse $\varphi_k$, and similarly, each $T_2|_{V_l}$ has an inverse $\psi_l$. An illustration of $T_1(x)$ is provided below:
We now apply the change of variables:
$$
\biggl\{\begin{aligned}
u &= aT_1(x) + by, \\
v &= cx + dT_2(y).
\end{aligned}
\tag{1}\label{e:1} $$
On each cell $U_k \times V_l$, and for fixed $(u, v)$, the system $\eqref{e:1}$ defines two curves:
$$ x = \varphi_k \left( \frac{u - by}{a} \right), \qquad y = \psi_l \left( \frac{v - cx}{d} \right). \tag{2}\label{e:2} $$
Since $T_1'(x) > 1$ and $T_2'(y) > 1$, the mapping $y \mapsto \varphi_k \bigl( \frac{u - by}{a} \bigr)$ is a $\frac{|b|}{a}$-Lipschitz map from $V_l$ to $U_k$, and simiarly, the mapping $x \mapsto \psi_l \bigl( \frac{v - cx}{d} \bigr)$ is a $\frac{|c|}{d}$-Lipschitz map from $U_k$ to $V_l$. Now using the conditions that $a, d > 0$ and $ad-bc>0$, it is not hard to show that $\eqref{e:2}$ has a unique solution $(x, y)$ for each given $(u, v)$. This allows defining $(x_{k,l}, y_{k,l})$ as a function of $(u, v)$ taking values in $U_k \times V_l$.
The figure below illustrates $(u, v)$ as a function of $(x, y)$. In the left plot, hue represents the argument of $u + iv$ and saturation represents its magnitude. The $xy$-plane is clearly partitioned into rectangles $U_k \times V_l$. The right plot depicts the two curves from $\eqref{e:2}$ for a specific choice of $(u, v)$.
Calculating the Jacobian
We now calculate the Jacobian determinant $\frac{\partial(x_{k,l}, y_{k,l})}{\partial(u, v)}$, suppressing the subscripts for brevity whenever it causes no ambiguity. Differentiating $\eqref{e:1}$ directly yields:
$$ (adT_1'(x)T_2'(y) - bc) \frac{\partial(x, y)}{\partial(u, v)} = 1. $$
Alternatively, by differentiating the rearranged form of $\eqref{e:1}$,
$$
\biggl\{\begin{aligned}
aT_1(x) &= u - by, \\
dT_2(y) &= v - cx,
\end{aligned} $$
we obtain:
$$ \begin{align*}
ad T_1'(x)T_2'(y) \frac{\partial(x, y)}{\partial(u, v)}
= \begin{vmatrix}
1 - b \frac{\partial y}{\partial u} &
-b \frac{\partial y}{\partial v} \\
-c\frac{\partial x}{\partial u} &
1-c\frac{\partial x}{\partial v}
\end{vmatrix}
\end{align*} $$
Simplifying this determinant gives:
$$ (adT_1'(x)T_2'(y) + bc) \frac{\partial(x, y)}{\partial(u, v)} = 1 - b \frac{\partial y}{\partial u} - c \frac{\partial x}{\partial v} $$
Combining these two identities and re-introducing subscripts, we arrive at:
$$ \frac{\partial(x_{k,l}, y_{k,l})}{\partial(u, v)} = - \frac{1}{2c} \frac{\partial y_{k,l}}{\partial u} - \frac{1}{2b} \frac{\partial x_{k,l}}{\partial v} \tag{3}\label{e:3} $$
Reduction
To utilize $\eqref{e:3}$, we propose the following:
Claim. The following summations hold:
$$ \sum_{k,l} \frac{\partial y_{k,l}}{\partial u} = -\frac{c}{ad-bc} \qquad\text{and} \qquad \sum_{k,l} \frac{\partial x_{k,l}}{\partial v} = -\frac{b}{ad-bc}. $$
Once this claim has been established, it follows that
$$\begin{align*}
\int_{\mathbb{R}^2} f(a T_1(x) + by, cx + dT_2(y)) \, \mathrm{d}x\mathrm{d}y
&= \sum_{k,l} \int_{U_k\times V_l} f(a T_1(x) + by, cx + dT_2(y)) \, \mathrm{d}x\mathrm{d}y \\
&= \sum_{k,l} \int_{\mathbb{R}^2} f(u, v) \, \frac{\partial(x_{k,l}, y_{k,l})}{\partial(u,v)} \mathrm{d}u\mathrm{d}v \\
&= \frac{1}{ad-bc} \int_{\mathbb{R}^2} f(u, v) \, \mathrm{d}u\mathrm{d}v
\end{align*}$$
as desired. So it remains to prove the claim:
Proof of Claim. Note that $\eqref{e:1}$ can be rearranged as
$$ by - u + aT_1\left( \frac{v - dT_2(y)}{c} \right) = 0. $$
Plugging the definition of $T_1$ and $T_2$ and manipulating a bit, it follows that
$$ (ad-bc)y + cu - av + \sum_{l=1}^{L} \frac{adq_l}{y - b_l} + \sum_{k=1}^{K} \frac{ac^2p_k}{v - dT_2(y) - a_kc} = 0. $$
This might look intimidating, but all we need to know is that this equation simplifies to
$$ (ad - bc)y + cu - av + \frac{P(y)}{Q(y)} = 0, $$
where $\deg P < \deg Q = K(L+1) + L$ and the coefficients of $P(y)$ and $Q(y)$ depend only on $v$. Also, assume for simplicity that $Q(y)$ is monic. Then the above equation can be further rearranged as
$$ (ad - bc)y Q(y) + (cu - av) Q(y) + P(y) = 0. $$
Writing $N = (K+1)(L+1)$ and Expanding the left-hand side,
$$ (ad - bc) y^N + (cu + \text{[function of $v$]})y^{N-1} + \mathcal{O}(y^{N-2}) = 0. $$
We already know that the solutions of this equation are $y_{k,l}$ for $k \in \{0,\ldots,K\}$ and $l \in \{0, \ldots, L\}$. Therefore, by the Vieta's formula,
$$ \sum_{k,l} y_{k,l} = -\frac{cu}{ad-bc} + \text{[function of $v$]} $$
and the first statement of the claim follows. The second claim follows by the symmetry.
