What is the measure of the set
$$ A=\bigcup_{\substack{1 \le z_x \le N_x \\ 1 \le z_y \le N_y}} A(z_x,z_y) $$ where $A(z_x,z_y)=\{ (x,y) \in [0,1]^2 : \lvert axz_x-byz_y \rvert \le c \}$ for some $c>0$ and where $a,b>1$?
$A(z_x,z_y)$ is a set of points between two parallel lines we can refer to as a strip.
The two figures below show examples of shapes of the region $A$.
Below is the shape of a region of small $c$.

Below is the shape of a region of large $c$.

What I tried: I used the following subadditivity (union) bound \begin{align} \mu(A)=\mu\left(\bigcup_{z_x,z_y} A(z_x,z_y)\right) \le \sum_{z_x,z_y} \mu(A(z_x,z_y)) \end{align}
then I was able to bound $\mu(A(z_x,z_y))$ (by finding its area) and found that $$\mu(A) \le 2c \min\left(\frac{N_x(1+\ln(N_y))}{b},\frac{N_y(1+\ln(N_x))}{a}\right).$$
However, I believe that a tighter bound can be found of the form $\mu(A) \le 2c \min\left(\frac{N_x}{b},\frac{N_y}{a}\right)$, because the subadditivity bound overcounts intersections. I was thinking about using the inclusion-exclusion formula but it seems very messy and difficult – but I can be wrong. Has anyone encountered something similar before? Are there better bounds on unions of events?
Thank you, in advance, for any help or ideas you can give me. Please, see an interesting approach by @epimorphic
If I understood your question correctly, the answer is in the negative: There is no bound of the form $Kc \min\left(\frac{N_x}{b}, \frac{N_y}{a}\right)$, where $K>0$ is a constant, that works for all values of parameters $a,b,c,N_x,N_y$.
We first note that the width of strip $A_{z_x, z_y}$ is $\left. c \middle/ \sqrt{a^2 z_x^2 + b^2 z_y^2} \right.$. Now for fixed $a$ and $b$, if $N_x = 1$, then for any $N_y$ we can always choose $c$ to be small enough (i.e. make the strips thin enough) that the strips are essentially disjoint rectangles (i.e. the intersections and the difference between strips and their approximating rectangles are insignificant in proportion to the total area). Since these approximating rectangles can each be chosen to have length at least one, we see that for any $N_y$ there always exists a $c$ such that $$\frac{1}{c}\mu(A) \gtrsim \sum_{z_y=1}^{N_y} \frac{1}{\sqrt{a^2 + b^2 z_y^2}}.$$ The right side is unbounded as $N_y \to \infty$, so we can find a pair $(N_y, c)$ who, together with the $A$ it defines, satisfies $$\mu(A) > \frac{Kc}{b} = Kc \min\left(\frac{N_x}{b}, \frac{N_y}{a}\right).$$
However, we do have a similar bound \begin{equation} \tag{*} \label{lazy-bound-74908} \mu(A) < \frac{\pi c}{\sqrt{2}} \sqrt{\frac{N_x^2}{b^2} + \frac{N_y^2}{a^2}} < \frac{\pi c}{\sqrt{2}} \left(\frac{N_x}{b} + \frac{N_y}{a}\right). \end{equation}
Proof. Each strip $A_{z_x, z_y}$ is contained in a rectangle of length $\sqrt{2}$ and width $\left. c \middle/ \sqrt{a^2 z_x^2 + b^2 z_y^2} \right.$, so $$ \mu(A) < 2^{1/2} c \sum_{\substack{1 \leq z_x \leq N_x \\ 1 \leq z_y \leq N_y}} \frac{1}{\sqrt{a^2 z_x^2 + b^2 z_y^2}}. $$ The sum can be bounded by a corresponding integral: $$ \sum_{\substack{1 \leq z_x \leq N_x \\ 1 \leq z_y \leq N_y}} \frac{1}{\sqrt{a^2 z_x^2 + b^2 z_y^2}} < \int\limits_{[0, N_x)\times[0, N_y)} \frac{dx\,dy}{\sqrt{a^2 x^2 + b^2 y^2}}. $$ The level sets of the integrand are quadrants of ellipses centered at the origin. Let $$Q := \{(x,y) \in \mathbb R_{\geq 0}^2 : a^2 x^2 + b^2 y^2 < a^2 N_x^2 + b^2 N_y^2\} \supset [0, N_x) \times [0, N_y).$$ The infinitesimal area of the elliptic level set strips, parametrized by their $x$-intercepts, is $dV(x) = \left(\frac{d}{dx} \left(\frac{\pi}{4} \cdot x \cdot \frac{a}{b}x\right)\right)dx = \frac{\pi a x}{2b}dx$.$^{**}$ Thus $$\begin{align*} \int\limits_{[0, N_x)\times[0, N_y)} \frac{dx\,dy}{\sqrt{a^2 x^2 + b^2 y^2}} &< \int_Q \frac{dx\,dy}{\sqrt{a^2 x^2 + b^2 y^2}} \\ &= \int_0^{\frac{1}{a}\sqrt{a^2 N_x^2 + b^2 N_y^2}} \frac{1}{\sqrt{a^2 x^2}} \cdot \frac{\pi a x}{2b} dx \\ &= \frac{\pi}{2ab} \sqrt{a^2 N_x^2 + b^2 N_y^2}. \end{align*}$$ Putting it all together gives $\eqref{lazy-bound-74908}$.
$^{**}$ As requested, I'll provide more details for this step. Denote the map $(x,y) \mapsto \frac{1}{\sqrt{a^2 x^2 + b^2 y^2}}$ on $\mathbb R_{\geq 0}^2 \setminus \{(0,0)\}$ by $f$, and for $\tilde x > 0$, let $Q(\tilde x) := \{(x,y) \in \mathbb R_{\geq 0}^2 : a^2 x^2 + b^2 y^2 < a^2 \tilde x^2\}$. Then $ f^{-1}\bigl(f(\tilde x + \epsilon,0), f(\tilde x,0)\bigr] = Q(\tilde x + \epsilon) \setminus Q(\tilde x). $ So intuitively (apologies for the incoming abuse of notation), the area of the infinitesimal level set strip $ f^{-1}\bigl(f(\tilde x + d\tilde x,0), f(\tilde x,0)\bigr] $ should be \begin{align*} \frac{\mu\left\{f^{-1}\bigl(f(\tilde x + d\tilde x,0), f(\tilde x,0)\bigr]\right\}}{d\tilde x} d\tilde x &= \frac{\mu(Q(\tilde x + d\tilde x)) - \mu(Q(\tilde x))}{d\tilde x}d\tilde x \\ &= \frac{d(\frac{\pi a}{4b} \tilde x^2)}{d\tilde x}d\tilde x \\ &= \frac{\pi a \tilde x}{2b} d\tilde x, \end{align*} where we used the fact that the area of an ellipse is $\pi$ times the product of the semi-axes (in particular, $\mu(Q(\tilde x)) = \frac{1}{4}(\pi \cdot \tilde x \cdot \frac{a}{b}\tilde x)$).
One way to make this rigorous is through the change of variables formula. We define an "elliptic-polar" coordinate system $(\tilde x(x,y), \theta(x,y))$, where $\tilde x(x,y) = \frac{1}{a}\sqrt{a^2 x^2 + b^2 y^2}$ is defined so that $(\tilde x, 0)$ lies on the same level set as $(x,y)$ (thus $\tilde x$ identifies level sets), and where $\theta$ is just the standard polar angle. By the polar parametrization of the ellipse, the reverse transformations are $$ x(\tilde x, \theta) = \frac{\frac{a}{b} \tilde x^2 \cos \theta}{\sqrt{(\frac{a}{b}\tilde x \cos \theta)^2 + (\tilde x \sin \theta)^2}} = \frac{a \tilde x \cos \theta}{\sqrt{(a \cos \theta)^2 + (b \sin \theta)^2}} $$ and $$ y(\tilde x, \theta) = \frac{a \tilde x \sin \theta}{\sqrt{(a \cos \theta)^2 + (b \sin \theta)^2}}, $$ which yields (according to Mathematica) the Jacobian determinant $$ \lvert \det(D(x,y))(\tilde x, \theta) \rvert = \frac{a^2 \tilde x}{(a \cos \theta)^2 + (b \sin \theta)^2}. $$ Thus \begin{align*} \int_Q f(x,y) \, dx \, dy &= \int_0^{\frac{1}{a}\sqrt{a^2 N_x^2 + b^2 N_y^2}} \left( \int_0^{\pi/2} \tilde f(\tilde x, \theta) \cdot \frac{a^2 \tilde x}{(a \cos \theta)^2 + (b \sin \theta)^2} d\theta \right) d\tilde x \\ &= \int_0^{\frac{1}{a}\sqrt{a^2 N_x^2 + b^2 N_y^2}} f(\tilde x, 0) \cdot a^2 \tilde x \cdot \left( \int_0^{\pi/2} \frac{d\theta}{(a \cos \theta)^2 + (b \sin \theta)^2} \right) d\tilde x \end{align*} and we have $$ \int_0^{\pi/2} \frac{d\theta}{(a \cos \theta)^2 + (b \sin \theta)^2} = \frac{\pi}{2ab} $$ (again according to Mathematica; we should be able to manually verify it using the tangent half-angle substitution) as desired.
And finally, some ruminations on intersections – be warned though that I have little to offer in the way of rigorous results. As I see it, there are two types of intersections. The first is the nesting of parallel strips: We have $$A(z_x, z_y) \subseteq A\left(\frac{z_x}{\gcd(z_x, z_y)}, \frac{z_y}{\gcd(z_x, z_y)}\right)$$ for all $z_x$ and $z_y$, leading to the sharper (but more complicated) bound $$ \mu(A) \leq \sum_{\substack{1 \leq z_x \leq N_x \\ 1 \leq z_y \leq N_y \\ \gcd(z_x, z_y) = 1}} \mu(A(z_x, z_y)) < 2^{1/2} c \sum_{\substack{1 \leq z_x \leq N_x \\ 1 \leq z_y \leq N_y \\ \gcd(z_x, z_y) = 1}} \frac{1}{\sqrt{a^2 z_x^2 + b^2 z_y^2}}. $$ Imagine now that we could take $N_x = N_y = \infty$. We would be able to write \begin{align*} \sum_{1\leq z_x, z_y} \frac{1}{\sqrt{a^2 z_x^2 + b^2 z_y^2}} &= \sum_{d=1}^\infty \sum_{\substack{1 \leq z_x, z_y \\ \gcd(z_x, z_y) = d}} \frac{1}{\sqrt{a^2 z_x^2 + b^2 z_y^2}} \\ &= \left(\sum_{d=1}^\infty \frac{1}{d}\right) \sum_{\substack{1 \leq k, l \\ \gcd(k, l) = 1}} \frac{1}{\sqrt{a^2 k^2 + b^2 l^2}}. \end{align*} Note that $\sum_{d=1}^N \frac{1}{d}$ is of order $O(\log N)$. This suggests that a similar procedure for finite $N_x$ and $N_y$ might be able to reduce the order of the estimate \eqref{lazy-bound-74908}. Obviously we can't factor cleanly the way we were able to above, only $$ \sum_{\substack{1 \leq z_x \leq N_x \\ 1 \leq z_y \leq N_y}} \frac{1}{\sqrt{a^2 z_x^2 + b^2 z_y^2}} = \sum_{d=1}^\infty \left(\frac{1}{d}\sum_{\substack{1 \leq k \leq N_x/d \\ 1 \leq l \leq N_y/d \\ \gcd(k, l) = 1}} \frac{1}{\sqrt{a^2 k^2 + b^2 l^2}}\right), $$ so it seems that we would have to carefully estimate the expression within the parentheses for this approach to work.
After discarding the redundant parallel strips, we would be left with the other type of intersection: those between nonparallel strips. These... I have no idea how to deal with.
Analytic number theory and combinatorial geometry, neither of which I have any background in, may be able to shed further light on these two types of intersections.