  Given a set $P$ of points and a set $U$ of axis-parallel unit squares in the
Euclidean plane, a minimum ply cover of $P$ with $U$ is a subset of $U$ that
covers $P$ and minimizes the number of squares that share a common
intersection, called the minimum ply cover number of $P$ with $U$. Biedl et al.
[Comput. Geom., 94:101712, 2020] showed that determining the minimum ply cover
number for a set of points by a set of axis-parallel unit squares is NP-hard,
and gave a polynomial-time 2-approximation algorithm for instances in which the
minimum ply cover number is constant. The question of whether there exists a
polynomial-time approximation algorithm remained open when the minimum ply
cover number is $\omega(1)$. We settle this open question and present a
polynomial-time $(8+\varepsilon)$-approximation algorithm for the general
problem, for every fixed $\varepsilon>0$.
