Merge ff9ca26 into dd353c6

bin/select_secondary

@@ -1,10 +1,9 @@
 #!/usr/bin/env python
 
-import sys
 import fitsio
 import numpy as np
-import argparse
 from desitarget.secondary import select_secondary, _get_scxdir
+from desitarget.brightmask import is_in_bright_mask, get_recent_mask_dir
 from desitarget import io
 import os
 from glob import glob
@@ -29,7 +28,14 @@ ap.add_argument("--scnddir",
 ap.add_argument("--writeall",
                 action='store_true',
                 help="Default behavior is to split targets by bright/dark-time surveys. Set this to ALSO write a file of ALL targets")
+ap.add_argument("--nomasking", action='store_true',
+                help="Masking occurs by default. If this is set, do NOT use a bright star mask to mask the secondary targets")
+ap.add_argument("--maskdir",
+                help="Name of the specific directory (or file) containing the bright star mask (defaults to the most recent directory in $MASK_DIR)",
+                default=None)
+
 ns = ap.parse_args()
+do_mask = not(ns.nomasking)
 
 # ADM Sanity check that priminfodir exists.
 if not os.path.exists(ns.priminfodir):
@@ -61,6 +67,27 @@ try:
 except (ValueError, IndexError, AttributeError):
     drint = "X"
 
+# ADM if secondaries need masked, grab the mask directory and find which
+# ADM targets are masked. Also add the masking information to the header.
+mdcomp = None
+if do_mask:
+    # ADM grab the mask directory
+    maskdir = get_recent_mask_dir(ns.maskdir)
+    log.info("Masking secondaries using bright star masks in {}".format(maskdir))
+    # ADM read in the masks across the entire sky.
+    Mx = io.read_targets_in_quick(maskdir, shape='box',
+                                  radecbox=[0.0, 360.0, -90.0, 90.0])
+    log.info("Checking {} targets against {} masks".format(len(scx), len(Mx)))
+    in_mask, _ = is_in_bright_mask(scx, Mx, inonly=True)
+    # ADM the output in_mask is a list. We want the array it contains.
+    in_mask = in_mask[0]
+    log.info("{} targets are in bright star masks".format(np.sum(in_mask)))
+    # ADM a compact version of the maskdir name.
+    md = maskdir.split("/")
+    mdcomp = "/".join(md[md.index("masks"):])
+hdr["masked"] = do_mask
+hdr["maskdir"] = mdcomp
+
 # ADM write out the secondary targets, with bright-time
 # ADM and dark-time targets written separately.
 obscons = ["BRIGHT", "DARK"]
@@ -69,9 +96,14 @@ if ns.writeall:
     obscons.append(None)
 for obscon in obscons:
     # ADM the dict() here is to make hdr immutable.
-    ntargs, outfile = io.write_secondary(ns.dest, scx, primhdr=dict(hdr),
-                                         scxdir=scxdir, obscon=obscon,
-                                         drint=drint)
+    if obscon is None:
+        ntargs, outfile = io.write_secondary(ns.dest, scx, primhdr=dict(hdr),
+                                             scxdir=scxdir, obscon=obscon,
+                                             drint=drint)
+    else:
+        ntargs, outfile = io.write_secondary(ns.dest, scx[~in_mask],
+                                             primhdr=dict(hdr), scxdir=scxdir,
+                                             obscon=obscon, drint=drint)
 
     log.info('{} standalone secondary targets written to {}...t={:.1f}mins'
              .format(ntargs, outfile, (time()-time0)/60.))

bin/select_skies

@@ -159,13 +159,13 @@ else:
         skies = resolved
 
     # ADM mask the sky locations using a bright star mask.
+    mdcomp = None
     if do_mask:
         maskdir = get_recent_mask_dir(ns.maskdir)
         skies = mask_targets(skies, maskdir, nside=ns.nside, pixlist=pixlist)
-
-    # ADM a compact version of the maskdir name.
-    md = maskdir.split("/")
-    mdcomp = "/".join(md[md.index("masks"):])
+        # ADM a compact version of the maskdir name.
+        md = maskdir.split("/")
+        mdcomp = "/".join(md[md.index("masks"):])
 
     # ADM extra header keywords for the output fits file.
     extra = {k: v for k, v in zip(["masked", "maskdir"],

doc/changes.rst

@@ -5,9 +5,21 @@ desitarget Change Log
 0.50.1 (unreleased)
 -------------------
 
+* Main Survey, secondary, outside-footprint target updates [`PR #677`_]:
+    * Update priorities for some secondary programs.
+    * MTL changes to reobserve 0.7 < z < 2.1 QSOs at low priority.
+        * Requested by two secondary programs.
+    * Allow secondaries to be truncated to first N targets in a file.
+    * Add fluxes/mags to the file of "standalone" secondaries.
+    * Apply bright star mask and a mag > 16 limit to secondary targets.
+    * Mask Main Survey MWS targets on BRIGHT and GALAXY using MASKBITS.
+    * Debug reading MTL ledger headers, now that they have units.
+    * Apply non-linear correction to de-extinct Gaia magnitudes TWICE.
+    * Add ``PHOT_G_N_OBS`` to the Gaia edr3 data model.
 * Migrated unit tests to GitHub Actions [`PR #675`_].
 
 .. _`PR #675`: https://github.com/desihub/desitarget/pull/675
+.. _`PR #677`: https://github.com/desihub/desitarget/pull/677
 
 
 0.50.0 (2021-01-29)

py/desitarget/cuts.py

@@ -652,7 +652,7 @@ def isSTD(gflux=None, rflux=None, zflux=None, primary=None,
           gaia=None, astrometricexcessnoise=None, paramssolved=None,
           pmra=None, pmdec=None, parallax=None, dupsource=None,
           gaiagmag=None, gaiabmag=None, gaiarmag=None, bright=False,
-          usegaia=True, south=True):
+          usegaia=True, maskbits=None, south=True):
     """Select STD targets using color cuts and photometric quality cuts (PSF-like
     and fracflux).  See isSTD_colors() for additional info.
 
@@ -703,9 +703,8 @@ def isSTD(gflux=None, rflux=None, zflux=None, primary=None,
         mask : boolean array, True if the object has colors like a STD star.
 
     Notes:
-        - Gaia-based quantities are as in `the Gaia data model`_.
-        - Current version (01/15/21) is version 236 on `the wiki`_.
-
+    - Gaia-based quantities are as in `the Gaia data model`_.
+    - Current version (02/18/21) is version 246 on `the wiki`_.
     """
     if primary is None:
         primary = np.ones_like(gflux, dtype='?')
@@ -724,6 +723,9 @@ def isSTD(gflux=None, rflux=None, zflux=None, primary=None,
     # ADM apply type=PSF cut
     std &= _psflike(objtype)
 
+    # ADM don't target standards in Legacy Surveys mask regions.
+    std &= imaging_mask(maskbits, mwsmask=True)
+
     # ADM apply fracflux, S/N cuts and number of observations cuts.
     fracflux = [gfracflux, rfracflux, zfracflux]
     fluxivar = [gfluxivar, rfluxivar, zfluxivar]
@@ -763,7 +765,7 @@ def isMWS_main(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
                obs_rflux=None, objtype=None, gaia=None,
                gaiagmag=None, gaiabmag=None, gaiarmag=None,
                gaiaaen=None, gaiadupsource=None, paramssolved=None,
-               primary=None, south=True):
+               primary=None, south=True, maskbits=None):
     """Set bits for main ``MWS`` targets.
 
     Args:
@@ -778,7 +780,7 @@ def isMWS_main(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
             ``True`` if and only if the object is a ``MWS_MAIN_BLUE`` target.
 
     Notes:
-        - as of 11/2/18, based on version 158 on `the wiki`_.
+    - Current version (02/18/21) is version 246 on `the wiki`_.
     """
     if primary is None:
         primary = np.ones_like(gaia, dtype='?')
@@ -789,7 +791,8 @@ def isMWS_main(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
     # if south:
     mws &= notinMWS_main_mask(gaia=gaia, gfracmasked=gfracmasked, gnobs=gnobs,
                               gflux=gflux, rfracmasked=rfracmasked, rnobs=rnobs,
-                              rflux=rflux, gaiadupsource=gaiadupsource, primary=primary)
+                              rflux=rflux, gaiadupsource=gaiadupsource,
+                              primary=primary, maskbits=maskbits)
 
     # ADM pass the mws that pass cuts as primary, to restrict to the
     # ADM sources that weren't in a mask/logic cut.
@@ -805,7 +808,7 @@ def isMWS_main(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
 
 
 def notinMWS_main_mask(gaia=None, gfracmasked=None, gnobs=None, gflux=None,
-                       rfracmasked=None, rnobs=None, rflux=None,
+                       rfracmasked=None, rnobs=None, rflux=None, maskbits=None,
                        gaiadupsource=None, primary=None):
     """Standard set of masking-based cuts used by MWS target selection classes
     (see, e.g., :func:`~desitarget.cuts.isMWS_main` for parameters).
@@ -814,6 +817,9 @@ def notinMWS_main_mask(gaia=None, gfracmasked=None, gnobs=None, gflux=None,
         primary = np.ones_like(gaia, dtype='?')
     mws = primary.copy()
 
+    # ADM don't target MWS-like targets in Legacy Surveys mask regions.
+    mws &= imaging_mask(maskbits, mwsmask=True)
+
     # ADM apply the mask/logic selection for all MWS-MAIN targets
     # ADM main targets match to a Gaia source
     mws &= gaia
@@ -941,7 +947,7 @@ def isMWS_nearby(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
 def isMWS_bhb(primary=None, objtype=None,
               gaia=None, gaiaaen=None, gaiadupsource=None, gaiagmag=None,
               gflux=None, rflux=None, zflux=None,
-              w1flux=None, w1snr=None,
+              w1flux=None, w1snr=None, maskbits=None,
               gnobs=None, rnobs=None, znobs=None,
               gfracmasked=None, rfracmasked=None, zfracmasked=None,
               parallax=None, parallaxerr=None):
@@ -961,7 +967,7 @@ def isMWS_bhb(primary=None, objtype=None,
     - Criteria supplied by Sergey Koposov
     - gflux, rflux, zflux, w1flux have been corrected for extinction
       (unlike other MWS selections, which use obs_flux).
-    - Current version (12/21/20) is version 149 on `the SV wiki`_.
+    - Current version (02/18/21) is version 246 on `the wiki`_.
     """
     if primary is None:
         primary = np.ones_like(gaia, dtype='?')
@@ -990,6 +996,9 @@ def isMWS_bhb(primary=None, objtype=None,
     gmr = gmag-rmag
     rmz = rmag-zmag
 
+    # ADM don't target MWS-like targets in Legacy Surveys mask regions.
+    mws &= imaging_mask(maskbits, mwsmask=True)
+
     # APC must be a Legacy Surveys object that matches a Gaia source
     mws &= gaia
     # APC type must be PSF
@@ -1121,8 +1130,9 @@ def isMWS_WD(primary=None, gaia=None, galb=None, astrometricexcessnoise=None,
     return mws
 
 
-def isMWSSTAR_colors(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None, primary=None, south=True):
-    """Select a reasonable range of g-r colors for MWS targets. Returns a boolean array.
+def isMWSSTAR_colors(gflux=None, rflux=None, zflux=None,
+                     w1flux=None, w2flux=None, primary=None, south=True):
+    """A reasonable range of g-r for MWS targets. Returns a boolean array.
 
     Args:
         gflux, rflux, zflux, w1flux, w2flux: array_like
@@ -2273,7 +2283,7 @@ def set_target_bits(photsys_north, photsys_south, obs_rflux,
                     w1flux=w1flux, w1snr=w1snr,
                     gnobs=gnobs, rnobs=rnobs, znobs=znobs,
                     gfracmasked=gfracmasked, rfracmasked=rfracmasked, zfracmasked=zfracmasked,
-                    parallax=parallax, parallaxerr=parallaxerr
+                    parallax=parallax, parallaxerr=parallaxerr, maskbits=maskbits
              )
 
         # ADM run the MWS target types for (potentially) both north and south.
@@ -2283,7 +2293,7 @@ def set_target_bits(photsys_north, photsys_south, obs_rflux,
                     gflux=gflux, rflux=rflux, obs_rflux=obs_rflux, objtype=objtype,
                     gnobs=gnobs, rnobs=rnobs, gfracmasked=gfracmasked,
                     rfracmasked=rfracmasked, pmra=pmra, pmdec=pmdec,
-                    parallax=parallax, parallaxerr=parallaxerr,
+                    parallax=parallax, parallaxerr=parallaxerr, maskbits=maskbits,
                     paramssolved=gaiaparamssolved, primary=primary, south=south
             )
     mws_broad_n, mws_red_n, mws_blue_n = mws_classes[0]
@@ -2311,7 +2321,7 @@ def set_target_bits(photsys_north, photsys_south, obs_rflux,
         for bright in [False, True]:
             std_classes.append(
                 isSTD(
-                    primary=primary, zflux=zflux, rflux=rflux, gflux=gflux,
+                    primary=primary, zflux=zflux, rflux=rflux, gflux=gflux, maskbits=maskbits,
                     gfracflux=gfracflux, rfracflux=rfracflux, zfracflux=zfracflux,
                     gfracmasked=gfracmasked, rfracmasked=rfracmasked, objtype=objtype,
                     zfracmasked=zfracmasked, gnobs=gnobs, rnobs=rnobs, znobs=znobs,

py/desitarget/data/targetmask.yaml

@@ -120,9 +120,9 @@ mws_mask:
 #- ADM variable, e.g. $SECONDARY_DIR/veto.fits for VETO targets.
 scnd_mask:
     - [VETO,        0, "Never observe, even if a primary target bit is set",
-       {obsconditions: DARK|GRAY|BRIGHT|POOR|TWILIGHT12|TWILIGHT18, filename: 'veto'}]
+       {obsconditions: DARK|GRAY|BRIGHT|POOR|TWILIGHT12|TWILIGHT18, filename: 'veto', downsample: 1}]
     - [DR16Q,       1, "Known quasars from the SDSS DR16Q catalog",
-       {obsconditions: DARK|GRAY, filename: 'dr16q'}]
+       {obsconditions: DARK|GRAY, filename: 'dr16q', downsample: 1}]
 
 #- Observing conditions
 #- These are a bitmask to allow target bits to specify multiple conditions
@@ -143,12 +143,16 @@ obsconditions:
 #- if a target passes more than one target bit, it is possible that one bit
 #- will consider the observations "DONE" while another bit will want "MORE".
 #- DONOTOBSERVE trumps MORE which trumps DONE
+# ADM The MORE_MIDZQSO priority is driven by secondary programs from Gontcho a Gontcho (1.4 < z < 2.1)
+# ADM and Weiner et al. (0.7 < z < 2.1) to reobserve confirmed quasars where possible. The priority
+# ADM of 100 should only be higher than DONE (and secondary filler) targets.
 obsmask:
-    - [UNOBS,       0, "unobserved"]
-    - [OBS,         1, "observed but no spectro pipeline answer yet"]
-    - [DONE,        2, "enough observations already obtained"]
-    - [MORE_ZWARN,  3, "ambiguous redshift; need more observations"]
-    - [MORE_ZGOOD,  4, "redshift known; need more observations"]
+    - [UNOBS,          0, "unobserved"]
+    - [OBS,            1, "observed but no spectro pipeline answer yet"]
+    - [DONE,           2, "enough observations already obtained"]
+    - [MORE_ZWARN,     3, "ambiguous redshift; need more observations"]
+    - [MORE_ZGOOD,     4, "redshift known; need more observations"]
+    - [MORE_MIDZQSO,   8, "z > 0.7 QSO; more observations at very low priority"]
     - [DONOTOBSERVE,  16, "Do not observe this target (possibly temporarily)"]
 
 #- Bits that can be set in TARGETID. See the schema at
@@ -175,9 +179,15 @@ targetid_mask:
 priorities:
     #- Dark Survey: priorities 3000 - 3999
     desi_mask:
-        ELG: {UNOBS: 3000, DONE: 2, OBS: 1, DONOTOBSERVE: 0}
-        LRG: {UNOBS: 3200, DONE: 2, OBS: 1, DONOTOBSERVE: 0}
-        QSO: {UNOBS: 3400, MORE_ZGOOD: 3350, MORE_ZWARN: 3300, DONE: 2, OBS: 1, DONOTOBSERVE: 0}
+# ADM safest to set MORE_ZGOOD for ELGs/LRGs to DONE PROVIDED they have NUMOBS=1 as
+# ADM they can match QSO targets that require multiple observations and trump those
+# ADM QSOs with a higher priority. There is a unit test to check NUMOBS=1 for ELGs/LRGs.
+        ELG: {UNOBS: 3000, MORE_ZGOOD: 2, DONE: 2, OBS: 1, DONOTOBSERVE: 0}
+        LRG: {UNOBS: 3200, MORE_ZGOOD: 2, DONE: 2, OBS: 1, DONOTOBSERVE: 0}
+# ADM The MORE_MIDZQSO priority is driven by secondary programs from Gontcho a Gontcho (1.4 < z < 2.1)
+# ADM and Weiner et al. (0.7 < z < 2.1) to reobserve confirmed quasars where possible. The priority
+# ADM of 100 should only be higher than DONE (and secondary filler) targets.
+        QSO: {UNOBS: 3400, MORE_ZGOOD: 3350, MORE_ZWARN: 3300, MORE_MIDZQSO: 100, DONE: 2, OBS: 1, DONOTOBSERVE: 0}
         # ADM don't prioritize a N/S target if it doesn't have other bits set
         LRG_NORTH: {UNOBS: 0, DONE: 0, OBS: 0, DONOTOBSERVE: 0}
         QSO_HIZ: SAME_AS_LRG_NORTH

py/desitarget/gaiamatch.py

@@ -7,7 +7,7 @@
 Useful Gaia matching and manipulation routines.
 
 .. _`Gaia Collaboration/Babusiaux et al. (2018)`: https://ui.adsabs.harvard.edu/abs/2018A%26A...616A..10G/abstract
-.. _`borrowed shamelessly from Sergey Koposov`: https://github.com/desihub/desispec/blob/e67142d5a3c5489ae861735b447443afc6211310/py/desispec/scripts/stdstars.py#L74
+.. _`borrowed shamelessly from Sergey Koposov`: https://github.com/desihub/desispec/blob/cd9af0dcc81c7c597aef2bc1c2a9454dcbc47e17/py/desispec/scripts/stdstars.py#L114
 """
 import os
 import sys
@@ -71,7 +71,7 @@
     ('PHOT_G_MEAN_MAG', '>f4'), ('PHOT_G_MEAN_FLUX_OVER_ERROR', '>f4'),
     ('PHOT_BP_MEAN_MAG', '>f4'), ('PHOT_BP_MEAN_FLUX_OVER_ERROR', '>f4'),
     ('PHOT_RP_MEAN_MAG', '>f4'), ('PHOT_RP_MEAN_FLUX_OVER_ERROR', '>f4'),
-    ('PHOT_BP_RP_EXCESS_FACTOR', '>f4'),
+    ('PHOT_BP_RP_EXCESS_FACTOR', '>f4'), ('PHOT_G_N_OBS', '>i4'),
     ('ASTROMETRIC_EXCESS_NOISE', '>f4'), ('ASTROMETRIC_EXCESS_NOISE_SIG', '>f4'),
     ('DUPLICATED_SOURCE', '?'), ('ASTROMETRIC_SIGMA5D_MAX', '>f4'),
     ('ASTROMETRIC_PARAMS_SOLVED', '>i1'), ('RUWE', '>f4'),
@@ -87,7 +87,7 @@
     ('EDR3_PHOT_G_MEAN_MAG', '>f4'), ('EDR3_PHOT_G_MEAN_FLUX_OVER_ERROR', '>f4'),
     ('EDR3_PHOT_BP_MEAN_MAG', '>f4'), ('EDR3_PHOT_BP_MEAN_FLUX_OVER_ERROR', '>f4'),
     ('EDR3_PHOT_RP_MEAN_MAG', '>f4'), ('EDR3_PHOT_RP_MEAN_FLUX_OVER_ERROR', '>f4'),
-    ('EDR3_PHOT_BP_RP_EXCESS_FACTOR', '>f4'),
+    ('EDR3_PHOT_BP_RP_EXCESS_FACTOR', '>f4'), ('EDR3_PHOT_G_N_OBS', '>i4'),
     ('EDR3_ASTROMETRIC_EXCESS_NOISE', '>f4'), ('EDR3_ASTROMETRIC_EXCESS_NOISE_SIG', '>f4'),
     ('EDR3_DUPLICATED_SOURCE', '?'), ('EDR3_ASTROMETRIC_SIGMA5D_MAX', '>f4'),
     ('EDR3_ASTROMETRIC_PARAMS_SOLVED', '>i1'), ('EDR3_RUWE', '>f4'),
@@ -236,15 +236,20 @@ def unextinct_gaia_mags(G, Bp, Rp, ebv, scaling=0.86):
 
     # ADM apply the extinction corrections in each band.
     gaia_a0 = 3.1 * ebv * scaling
-    bprp = Bp - Rp
-    for band in ['G', 'BP', 'RP']:
-        curp = gaia_poly_coeff[band]
-        dmag = (
-            np.poly1d(gaia_poly_coeff[band][:4][::-1])(bprp) +
-            curp[4]*gaia_a0 + curp[5]*gaia_a0**2 + curp[6]*bprp*gaia_a0
-        )*gaia_a0
-        # ADM populate the per-band extinction-corrected magnitudes.
-        outmags[band] = inmags[band] - dmag
+
+    for i in range(2):
+        if i == 0:
+            bprp = Bp - Rp
+        else:
+            bprp = outmags["BP"] - outmags["RP"]
+        for band in ['G', 'BP', 'RP']:
+            curp = gaia_poly_coeff[band]
+            dmag = (
+                np.poly1d(gaia_poly_coeff[band][:4][::-1])(bprp) +
+                curp[4]*gaia_a0 + curp[5]*gaia_a0**2 + curp[6]*bprp*gaia_a0
+            )*gaia_a0
+            # ADM populate the per-band extinction-corrected magnitudes.
+            outmags[band] = inmags[band] - dmag
 
     return outmags["G"], outmags["BP"], outmags["RP"]
 
@@ -437,7 +442,7 @@ def gaia_csv_to_fits(dr="dr2", numproc=32):
         os.makedirs(fitsdir)
 
     # ADM construct the list of input files.
-    infiles = glob("{}/*csv*".format(csvdir))
+    infiles = glob("{}/GaiaSource*csv*".format(csvdir))
     nfiles = len(infiles)
 
     # ADM the critical function to run on every file.
@@ -547,7 +552,7 @@ def gaia_fits_to_healpix(dr="dr2", numproc=32):
         - The environment variable $GAIA_DIR must be set.
         - if numproc==1, use the serial code instead of the parallel code.
         - Runs in 1-2 hours with numproc=32 for 61,234 Gaia DR2 files.
-        - Runs in 1-2 hours with numproc=32 for 3,386 Gaia EDR3 files.
+        - Runs in ~15 minutes with numproc=32 for 3,386 Gaia EDR3 files.
     """
     # ADM the resolution at which the Gaia HEALPix files should be stored.
     nside = _get_gaia_nside()