Let $F$ be a bounded subset of $\mathbb R^d$ with $d \geq 2$. We define $F_r := rF \cap \mathbb Z^d$ for any real number $r>0$ and assume that the limit
$$
\mathcal V(F) := \lim_{r \to + \infty} \frac{|F_r|}{r^d}
$$
exists (for convex subsets, this is the Lebesgue volume of $F$). We rewrite the proof of Theorem 459 of Hardy–Wright so that it yields the following more general result.

*If $\mathcal V(F)$ is well-defined, then we have*
$$
\lim_{r \rightarrow + \infty} \frac{\left| \left\{ x \in F_r, \operatorname{gcd}(x_1, \cdots, x_d) = 1 \right\} \right|}{r^d} = \frac{\mathcal V(F)}{\zeta(d)}.
$$

**Proof.**
We can and will assume that $0 \notin F$, which will not change any of the limits.
We also fix $N$ such that $F \subset [-N,N]^d$.

For every rational $r>0$, let $f(r) = \left| \left\{ x \in F_r, \operatorname{gcd}(x_1, \cdots, x_d) = 1 \right\} \right|$. As $0 \notin F$, $|F_r| = f(r)=0$ when $r<1/N$ and $f(r) \leq |F_r| \leq (2rN+1)^d \leq (3rN)^d$ for all $r \geq 1/N$, so $|F_r| \leq (3rN)^d$ in all cases. For any point $x$ of $F_r$, there is a unique integer $k \in \mathbb N$ such that the gcd of the coordinates of $x$ is $k$, and then $x/k$ contributes to $f(r/k)$. Consequently (the right-hand side being in fact a finite sum)
$$
|F_r| = \sum_{k=1}^{+ \infty} f(r/k).
$$
By Möbius inversion, we then get
$$
f(r) = \sum_{k=1}^{+ \infty} \mu(k) |F_{r/k}|.
$$

The sum of $\mu(k)/k^d$ converges absolutely towards $1/\zeta(d)$ as $d \geq 2$, so
$$
\frac{f(r)}{r^d} - \frac{\mathcal V(F)}{\zeta(d)} = \sum_{k=1}^{+ \infty} \frac{\mu(k)}{k^d} \left( \frac{|F_{r/k}|}{(r/k)^d} - \mathcal V(F) \right).
$$
Let $\varepsilon>0$. By definition of $\mathcal V(F)$, we fix $n_0$ such that if $r/k \geq n_0$, $\left| \frac{|F_{r/k}|}{(r/k)^d} - \mathcal V(F) \right| \leq \varepsilon$, which gives the inequality
$$
\sum_{k=1}^{\lfloor r/n_0 \rfloor} \frac{1}{k^d} \left| \frac{|F_{r/k}|}{(r/k)^d} - \mathcal V(F) \right| \leq \zeta(d) \varepsilon.
$$

On the other hand, the bounds on $|F_{r/k}|$ give
$$
\sum_{k > \lfloor r/n_0 \rfloor} \frac{1}{k^d} \left| \frac{|F_{r/k}|}{(r/k)^d} - \mathcal V(F) \right| \leq \left((3N)^d + \mathcal V(F)\right) \times \frac{ (\lfloor r/n_0 \rfloor)^{1-d} }{d-1}.
$$
Hence, for $r$ large enough, the absolute value of $\left|\frac{f(r)}{r^d} - \frac{\mathcal V(F)}{\zeta(d)}\right|$ is smaller than $2 \zeta(d) \varepsilon$, which proves the desired convergence. $\blacksquare$

*Acknowledgements.* This post greatly benefitted from exchanges with Samuel Le Fourn.