Merge pull request #805 from desihub/ADM-cols-mtl-hp

Add an option to only read a subset of columns from the HEALPixel-split MTLs

doc/changes.rst

@@ -6,7 +6,15 @@ desitarget Change Log
 ------------------
 
 * Reproducible output from :func:`io.read_targets_in_hp()` [`PR #806`_].
-
+* Option to read subset of columns from HEALPix-split MTLs [`PR #805`_].
+    * Specifically in :func:`io.read_mtl_in_hp()`.
+    * Addresses `issue #802`_.
+    * Partial fix for `issue #804`_.
+    * Also address some deprecation warnings.
+
+.. _`issue #802`: https://github.com/desihub/desitarget/issues/802
+.. _`issue #804`: https://github.com/desihub/desitarget/issues/804
+.. _`PR #805`: https://github.com/desihub/desitarget/pull/805
 .. _`PR #806`: https://github.com/desihub/desitarget/pull/806
 
 2.6.1 (2023-08-23)

py/desitarget/QA.py

@@ -981,7 +981,7 @@ def qagaia(cat, objtype, qadir='.', fileprefix="gaia", nobjscut=1000, seed=None)
     plt.xlabel(r'$PM_{RA}\,(mas\,yr^{-1})$')
     plt.ylabel(r'$PM_{DEC}\,(mas\,yr^{-1})$')
 
-    cmap = plt.cm.get_cmap('RdYlBu')
+    cmap = matplotlib.colormaps.get_cmap('RdYlBu')
 
     ralim = (-25, 25)
     declim = (-25, 25)
@@ -1244,7 +1244,7 @@ def mock_qanz(cat, objtype, qadir='.', area=1.0, dndz=None, nobjscut=1000,
     plt.clf()
     plt.xlabel(zlabel)
 
-    cmap = plt.cm.get_cmap('RdYlBu')
+    cmap = matplotlib.colormaps.get_cmap('RdYlBu')
 
     if 'LRG' in objtype:
         band = 'z'
@@ -1444,7 +1444,7 @@ def qso_colorbox(ax, plottype='grz', verts=None):
     zW1lim = (-1.5, 3.5)
     W1W2lim = (-1.5, 1.5)
 
-    cmap = plt.cm.get_cmap('RdYlBu')
+    cmap = matplotlib.colormaps.get_cmap('RdYlBu')
 
     objcolor = {'ALL': 'black', objtype: 'black'}
     type2color = {**_type2color, **objcolor}

py/desitarget/io.py

@@ -142,8 +142,8 @@ def add_photsys(indata):
         # ADM add PHOTSYS to the data model.
         # ADM the fitsio check is a hack for the v0.9 to v1.0 transition
         # ADM (v1.0 now converts all byte strings to unicode strings).
-        from distutils.version import LooseVersion
-        if LooseVersion(fitsio.__version__) >= LooseVersion('1'):
+        from packaging.version import Version
+        if Version(fitsio.__version__) >= Version('1'):
             pdt = [('PHOTSYS', '<U1')]
         else:
             pdt = [('PHOTSYS', '|S1')]
@@ -2041,12 +2041,12 @@ def check_fitsio_version(version='0.9.8'):
     ImportError
         If the fitsio version is insufficiently recent.
     """
-    from distutils.version import LooseVersion
+    from packaging.version import Version
     #
     # LooseVersion doesn't handle rc1 as we want, so also check for 0.9.8xxx.
     #
     if (
-        LooseVersion(fitsio.__version__) < LooseVersion(version) and
+        Version(fitsio.__version__) < Version(version) and
         not fitsio.__version__.startswith(version)
     ):
         raise ImportError(('ERROR: fitsio >{0}rc1 required ' +
@@ -2794,7 +2794,7 @@ def write_mtl_tile_file(filename, data):
 
 
 def read_mtl_ledger(filename, unique=True, isodate=None,
-                    initial=False, leq=False):
+                    initial=False, leq=False, columns=None):
     """Wrapper to read individual MTL ledger files.
 
     Parameters
@@ -2820,6 +2820,8 @@ def read_mtl_ledger(filename, unique=True, isodate=None,
         If ``True``, restrict the ledger to entries BEFORE or EQUAL TO
         `isodate` instead of the default behavior of strictly before
         `isodate`. Only relevant if `isodate` is passed.
+    columns : :class:`list`, optional
+        Only return these target columns.
 
     Returns
     -------
@@ -2897,9 +2899,17 @@ def read_mtl_ledger(filename, unique=True, isodate=None,
         if not initial:
             mtl = np.flip(mtl)
         _, ii = np.unique(mtl["TARGETID"], return_index=True)
-        return mtl[ii]
-    else:
-        return mtl
+        mtl = mtl[ii]
+
+    # ADM limit to certain columns, if only a subset were requested.
+    if columns is not None:
+        dt = [col for col in mtl.dtype.descr if col[0] in columns]
+        mtlredux = np.zeros(len(mtl), dtype=dt)
+        for col in mtlredux.dtype.names:
+            mtlredux[col] = mtl[col]
+        return mtlredux
+
+    return mtl
 
 
 def read_target_files(filename, columns=None, rows=None, header=False,
@@ -3132,7 +3142,7 @@ def find_mtl_file_format_from_header(hpdirname, returnoc=False,
 
 
 def read_mtl_in_hp(hpdirname, nside, pixlist, unique=True, isodate=None,
-                   returnfn=False, initial=False, leq=False):
+                   returnfn=False, initial=False, leq=False, columns=None):
     """Read Merged Target List ledgers in a set of HEALPixels.
 
     Parameters
@@ -3166,6 +3176,10 @@ def read_mtl_in_hp(hpdirname, nside, pixlist, unique=True, isodate=None,
         If ``True``, restrict the ledger to entries BEFORE or EQUAL TO
         `isodate` instead of the default behavior of strictly before
         `isodate`. Only relevant if `isodate` is passed.
+    columns : :class:`list`, optional
+        Only return these target columns. `RA` and `DEC` will always be
+        included in the output, whether or not they're passed, as they
+        are critical for determining if a location is in a HEALPixel.
 
     Returns
     -------
@@ -3184,6 +3198,12 @@ def read_mtl_in_hp(hpdirname, nside, pixlist, unique=True, isodate=None,
     if isinstance(pixlist, int):
         pixlist = [pixlist]
 
+    # ADM if columns was passed without "RA" and "DEC", add them.
+    if columns is not None:
+        for radec in ["RA", "DEC"]:
+            if radec not in columns:
+                columns.append(radec)
+
     # ADM if a directory was passed, do fancy HEALPixel parsing...
     outfns = {}
     fileform = find_mtl_file_format_from_header(hpdirname)
@@ -3198,8 +3218,9 @@ def read_mtl_in_hp(hpdirname, nside, pixlist, unique=True, isodate=None,
         for pix in filepixlist:
             fn = fileform.format(pix)
             try:
-                targs = read_mtl_ledger(
-                    fn, unique=unique, isodate=isodate, initial=initial, leq=leq)
+                targs = read_mtl_ledger(fn, unique=unique, isodate=isodate,
+                                        initial=initial, leq=leq,
+                                        columns=columns)
                 mtls.append(targs)
                 outfns[pix] = fn
             except FileNotFoundError:
@@ -3209,7 +3230,7 @@ def read_mtl_in_hp(hpdirname, nside, pixlist, unique=True, isodate=None,
         if len(mtls) == 0:
             fns = iglob(fileform.format("*"))
             fn = next(fns)
-            mtl = read_mtl_ledger(fn)
+            mtl = read_mtl_ledger(fn, columns=columns)
             outly = np.zeros(0, dtype=mtl.dtype)
             if returnfn:
                 return outly, outfns
@@ -3218,8 +3239,8 @@ def read_mtl_in_hp(hpdirname, nside, pixlist, unique=True, isodate=None,
         mtl = np.concatenate(mtls)
     # ADM ...if a directory wasn't passed, just read in the targets.
     else:
-        mtl = read_mtl_ledger(
-            hpdirname, unique=unique, isodate=isodate, initial=initial, leq=leq)
+        mtl = read_mtl_ledger(hpdirname, unique=unique, isodate=isodate,
+                              initial=initial, leq=leq, columns=columns)
 
     # ADM restrict the targets to the actual requested HEALPixels...
     ii = is_in_hp(mtl, nside, pixlist)

py/desitarget/randoms.py

@@ -517,8 +517,8 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir,
                             imsigma[ivar == 0] = 0
                         # ADM aperture photometry at requested radius (aprad).
                         apxy = np.vstack((x, y)).T
-                        aper = photutils.CircularAperture(apxy, aprad)
-                        p = photutils.aperture_photometry(img.data, aper, error=imsigma)
+                        aper = photutils.aperture.CircularAperture(apxy, aprad)
+                        p = photutils.aperture.aperture_photometry(img.data, aper, error=imsigma)
                         # ADM store the results.
                         qdict[qout+'_'+filt] = np.array(p.field('aperture_sum'))
                         err = p.field('aperture_sum_err')

py/desitarget/skyfibers.py

@@ -472,8 +472,8 @@ def sky_fibers_for_brick(survey, brickname, nskies=144, bands=['g', 'r', 'z'],
             imsigma[coiv == 0] = 0
         apxy = np.vstack((skyfibers.x, skyfibers.y)).T
         for irad, rad in enumerate(apertures):
-            aper = photutils.CircularAperture(apxy, rad)
-            p = photutils.aperture_photometry(coimg, aper, error=imsigma)
+            aper = photutils.aperture.CircularAperture(apxy, rad)
+            p = photutils.aperture.aperture_photometry(coimg, aper, error=imsigma)
             apflux[:, irad] = p.field('aperture_sum')
             err = p.field('aperture_sum_err')
             # ADM where the error is 0, that actually means infinite error