Merge ded7a8e into 05db498

doc/changes.rst

@@ -11,6 +11,11 @@ desitarget Change Log
     * Change the overall density from 800/sq.deg. to 600/sq.deg.
     * Remove the LRG_LOWDENS target bit.
 * Add ``desispec.skybricks`` to lookup if ras,decs are blank sky [`PR #722`_].
+* Update MWS cuts to Gaia EDR3 [`PR #720`_]. Includes:
+   * AEN stellarity cut now 2 (previously 3).
+   * Parallax floor now 0.3 mas (previously 1 mas).
+   * RED/BROAD proper motion split now a function of magnitude.
+   * ASTROMETRIC_PARAMSS_SOVLED checks account for range of values in EDR3.
 * Add `MWS-FAINT-RED` and `MWS-FAINT-BLUE` target classes for main survey
   [`PR #719`_].
 * Prepare the Main Survey cuts and bit-masks [`PR #718`_]. Includes:
@@ -39,6 +44,7 @@ desitarget Change Log
 .. _`PR #717`: https://github.com/desihub/desitarget/pull/717
 .. _`PR #718`: https://github.com/desihub/desitarget/pull/718
 .. _`PR #719`: https://github.com/desihub/desitarget/pull/719
+.. _`PR #720`: https://github.com/desihub/desitarget/pull/720
 .. _`PR #722`: https://github.com/desihub/desitarget/pull/722
 .. _`PR #723`: https://github.com/desihub/desitarget/pull/723
 

py/desitarget/cuts.py

@@ -182,7 +182,7 @@ def isGAIA_STD(ra=None, dec=None, galb=None, gaiaaen=None, pmra=None, pmdec=None
                parallax=None, parallaxovererror=None, ebv=None, gaiabprpfactor=None,
                gaiasigma5dmax=None, gaiagmag=None, gaiabmag=None, gaiarmag=None,
                gaiadupsource=None, gaiaparamssolved=None,
-               primary=None, test=False, nside=2):
+               primary=None, test=False, nside=2, dr="dr2"):
     """Standards based solely on Gaia data.
 
     Parameters
@@ -195,6 +195,8 @@ def isGAIA_STD(ra=None, dec=None, galb=None, gaiaaen=None, pmra=None, pmdec=None
     nside : :class:`int`, optional, defaults to 2
         (NESTED) HEALPix nside, if targets are being parallelized.
         The default of 2 should be benign for serial processing.
+    dr : :class:`str`, optional, defaults to "dr2"
+        Name of a Gaia data release. Options are "dr2", "edr3".
 
     Returns
     -------
@@ -237,7 +239,7 @@ def isGAIA_STD(ra=None, dec=None, galb=None, gaiaaen=None, pmra=None, pmdec=None
                       gaiagmag=gaiagmag, gaiabmag=gaiabmag, gaiarmag=gaiarmag)
 
     # ADM restrict to point sources.
-    ispsf = gaia_psflike(gaiaaen, gaiagmag)
+    ispsf = gaia_psflike(gaiaaen, gaiagmag, dr=dr)
     std &= ispsf
 
     # ADM de-extinct the magnitudes before applying color cuts.
@@ -672,7 +674,7 @@ def isSTD_gaia(primary=None, gaia=None, astrometricexcessnoise=None,
 
     # ADM no obvious issues with the astrometry solution.
     std &= astrometricexcessnoise < 1
-    std &= paramssolved == 31
+    std &= paramssolved >= 31
 
     # ADM finite proper motions.
     std &= np.isfinite(pmra)
@@ -840,7 +842,7 @@ def isMWS_main(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
             ``True`` if the object is a ``MWS_FAINT_BLUE`` target.
 
     Notes:
-    - Current version (05/04/21) is version 253 on `the wiki`_.
+    - Current version (05/04/21) is version 256 on `the wiki`_.
     """
     if primary is None:
         primary = np.ones_like(gaia, dtype='?')
@@ -913,7 +915,7 @@ def isMWS_faint_colors(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=N
     # APC main targets are point-like based on DECaLS morphology
     # APC and GAIA_ASTROMETRIC_NOISE.
     faint &= _psflike(objtype)
-    faint &= gaiaaen < 3.0
+    faint &= gaiaaen < 2.0
 
     # APC faint targets are 19 <= r < 20
     faint &= rflux > 10**((22.5-20.0)/2.5)
@@ -935,23 +937,23 @@ def isMWS_faint_colors(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=N
     # ADM Turn off any NaNs for astrometric quantities to suppress
     # ADM warnings. Won't target these, using cuts on paramssolved
     # ADM (or will explicitly target them based on paramsssolved).
-    ii = paramssolved != 31
+    ii = paramssolved < 31
     parallax = parallax.copy()
     parallax[ii] = 0.
 
     # ADM MWS-FAINT-RED has g-r >= 0.7
     faint_red &= rflux >= gflux * 10**(0.7/2.5)  # (g-r)>=0.7
 
-    # APC MWS-FAINT-RED also has parallax < max(3parallax_err,1) mas
-    # APC and proper motion < 3 (lower than MAIN-MAIN-RED)
+    # APC MWS-FAINT-RED has samle plx cut as MAIN-RED
+    # APC and proper motion < 3 (lower than MAIN-RED)
     # ACP and all astrometric parameters are measured.
-    faint_red &= parallax < np.maximum(3*parallaxerr, 1)
+    faint_red &= parallax < (3 * parallaxerr + 0.3)
     faint_red &= pm < 3.
-    faint_red &= paramssolved == 31
+    faint_red &= paramssolved >= 31
 
-    # APC There is no equivalent of MWS-MAIN-BROAD for the faint extension, any
-    # APC stars failing the astrometry cuts or without the Gaia parameters are
-    # APC discarded.
+    # APC There is no equivalent of MWS-MAIN-BROAD for the faint sample.
+    # APC stars failing the astrometry cuts or without the Gaia parameters
+    # APC are discarded.
 
     return faint_red, faint_blue
 
@@ -970,7 +972,7 @@ def isMWS_main_colors(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=No
     # ADM main targets are point-like based on DECaLS morphology
     # ADM and GAIA_ASTROMETRIC_NOISE.
     mws &= _psflike(objtype)
-    mws &= gaiaaen < 3.0
+    mws &= gaiaaen < 2.0
 
     # ADM main targets are 16 <= r < 19
     mws &= rflux > 10**((22.5-19.0)/2.5)
@@ -992,27 +994,29 @@ def isMWS_main_colors(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=No
     # ADM Turn off any NaNs for astrometric quantities to suppress
     # ADM warnings. Won't target these, using cuts on paramssolved
     # ADM (or will explicitly target them based on paramsssolved).
-    ii = paramssolved != 31
+    ii = paramssolved < 31
     parallax = parallax.copy()
     parallax[ii], pm[ii] = 0., 0.
 
     # ADM MWS-RED and MWS-BROAD have g-r >= 0.7
     red &= rflux >= gflux * 10**(0.7/2.5)  # (g-r)>=0.7
     broad = red.copy()
 
-    # ADM MWS-RED also has parallax < max(3parallax_err,1)mas
-    # ADM and proper motion < 7
+    # ADM MWS-RED also has parallax < (3parallax_err+0.3)mas
+    # ADM and proper motion < 5 * sqrt(rflux/r_19)
     # ADM and all astrometric parameters are measured.
-    red &= parallax < np.maximum(3*parallaxerr, 1)
-    red &= pm < 7.
-    red &= paramssolved == 31
-
-    # ADM MWS-BROAD has parallax > max(3parallax_err,1)mas
-    # ADM OR proper motion > 7.
+    pmmax = 5 * np.sqrt(rflux/(10**((22.5-19)/2.5)))
+    plxmax = 3 * parallaxerr + 0.3
+    red &= parallax < plxmax
+    red &= pm < pmmax
+    red &= paramssolved >= 31
+
+    # ADM MWS-BROAD has parallax > (3parallax_err+0.3)mas
+    # ADM OR proper motion > max value
     # ADM OR astrometric parameters not measured.
-    broad &= ((parallax >= np.maximum(3*parallaxerr, 1)) |
-              (pm >= 7.)
-              | (paramssolved != 31))
+    broad &= ((parallax >= plxmax) |
+              (pm >= pmmax)
+              | (paramssolved < 31))
 
     return broad, red, blue
 
@@ -1057,7 +1061,7 @@ def isMWS_nearby(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None,
     # APC Gaia G mag of more than 16.
     mws &= gaiagmag > 16.
     # ADM all astrometric parameters are measured.
-    mws &= paramssolved == 31
+    mws &= paramssolved >= 31
     # ADM parallax cut corresponding to 100pc.
     mws &= (parallax + parallaxerr) > 10.   # NB: "+" is correct.
 
@@ -1136,8 +1140,8 @@ def isMWS_bhb(primary=None, objtype=None,
 
     # APC no gaia duplicated sources.
     mws &= ~gaiadupsource
-    # APC gaia astrometric excess noise < 3.
-    mws &= gaiaaen < 3.0
+    # APC gaia astrometric excess noise < 2.
+    mws &= gaiaaen < 2.0
 
     # APC BHB extinction-corrected color range -0.35 <= gmr <= -0.02.
     mws &= (gmr >= -0.35) & (gmr <= -0.02)
@@ -1208,7 +1212,7 @@ def isMWS_WD(primary=None, gaia=None, galb=None, astrometricexcessnoise=None,
     mws &= gaia
 
     # ADM and all astrometric parameters are measured.
-    mws &= paramssolved == 31
+    mws &= paramssolved >= 31
 
     # ADM Gaia G mag of less than 20
     mws &= gaiagmag < 20.
@@ -2785,6 +2789,9 @@ def apply_cuts_gaia(numproc=4, survey='main', nside=None, pixlist=None,
     if survey == 'main':
         import desitarget.cuts as targcuts
         from desitarget.targetmask import desi_mask, mws_mask
+
+        # APC Use Gaia EDR3 in main survey
+        dr = "edr3"
     elif survey[:2] == 'sv':
         try:
             targcuts = import_module("desitarget.{}.{}_cuts".format(survey,
@@ -2797,6 +2804,9 @@ def apply_cuts_gaia(numproc=4, survey='main', nside=None, pixlist=None,
             log.critical(msg)
             raise ModuleNotFoundError(msg)
         desi_mask, mws_mask = targmask.desi_mask, targmask.mws_mask
+
+        # We used DR2 in SV
+        dr = "dr2"
     else:
         msg = "survey must be either 'main'or 'svX', not {}!!!".format(survey)
         log.critical(msg)
@@ -2829,13 +2839,25 @@ def apply_cuts_gaia(numproc=4, survey='main', nside=None, pixlist=None,
 
     # ADM determine if a target is a Gaia-only standard.
     primary = np.ones_like(gaiaobjs, dtype=bool)
-    std_faint, std_bright, std_wd = targcuts.isGAIA_STD(
-        ra=ra, dec=dec, galb=galb, gaiaaen=gaiaaen, pmra=pmra, pmdec=pmdec,
-        parallax=parallax, parallaxovererror=parallaxovererror, ebv=ebv,
-        gaiabprpfactor=gaiabprpfactor, gaiasigma5dmax=gaiasigma5dmax,
-        gaiagmag=gaiagmag, gaiabmag=gaiabmag, gaiarmag=gaiarmag,
-        gaiadupsource=gaiadupsource, gaiaparamssolved=gaiaparamssolved,
-        primary=primary, nside=nside, test=test)
+
+    if survey == 'main':
+        # APC in this case have a dr option for isGAIA_STD
+        std_faint, std_bright, std_wd = targcuts.isGAIA_STD(
+            ra=ra, dec=dec, galb=galb, gaiaaen=gaiaaen, pmra=pmra, pmdec=pmdec,
+            parallax=parallax, parallaxovererror=parallaxovererror, ebv=ebv,
+            gaiabprpfactor=gaiabprpfactor, gaiasigma5dmax=gaiasigma5dmax,
+            gaiagmag=gaiagmag, gaiabmag=gaiabmag, gaiarmag=gaiarmag,
+            gaiadupsource=gaiadupsource, gaiaparamssolved=gaiaparamssolved,
+            primary=primary, nside=nside, test=test, dr=dr)
+    else:
+        # SV versions of isGAIA_STD don't have a dr option
+        std_faint, std_bright, std_wd = targcuts.isGAIA_STD(
+            ra=ra, dec=dec, galb=galb, gaiaaen=gaiaaen, pmra=pmra, pmdec=pmdec,
+            parallax=parallax, parallaxovererror=parallaxovererror, ebv=ebv,
+            gaiabprpfactor=gaiabprpfactor, gaiasigma5dmax=gaiasigma5dmax,
+            gaiagmag=gaiagmag, gaiabmag=gaiabmag, gaiarmag=gaiarmag,
+            gaiadupsource=gaiadupsource, gaiaparamssolved=gaiaparamssolved,
+            primary=primary, nside=nside, test=test)
 
     # ADM Construct the target flag bits.
     mws_target = backup_bright * mws_mask.BACKUP_BRIGHT

py/desitarget/gaiamatch.py

@@ -164,7 +164,7 @@ def get_gaia_nside_brick(bricksize=0.25):
     return pixarea2nside(bricksize*bricksize)
 
 
-def gaia_psflike(aen, g):
+def gaia_psflike(aen, g, dr="dr2"):
     """Whether an objects is PSF-like based on Gaia quantities.
 
     Parameters
@@ -173,6 +173,8 @@ def gaia_psflike(aen, g):
         Gaia Astrometric Excess Noise.
     g : :class:`array_like` or :class`float`
         Gaia-based g MAGNITUDE (not Galactic-extinction-corrected).
+    dr : :class:`str`, optional, defaults to "dr2"
+        Name of a Gaia data release. Options are "dr2", "edr3"
 
     Returns
     -------
@@ -184,10 +186,23 @@ def gaia_psflike(aen, g):
     -----
         - Input quantities are the same as in `the Gaia data model`_.
     """
-    psflike = np.logical_or(
-        (g <= 19.) * (aen < 10.**0.5),
-        (g >= 19.) * (aen < 10.**(0.5 + 0.2*(g - 19.)))
-    )
+    # ADM allowed Data Releases for input.
+    droptions = ["dr2", "edr3"]
+    if dr not in droptions:
+        msg = "input dr must be one of {}".format(droptions)
+        log.critical(msg)
+        raise IOError(msg)
+
+    if dr == "dr2":
+        psflike = np.logical_or(
+            (g <= 19.) * (aen < 10.**0.5),
+            (g >= 19.) * (aen < 10.**(0.5 + 0.2*(g - 19.)))
+        )
+    elif dr == "edr3":
+        psflike = np.logical_or(
+            (g <= 19.) * (aen < 10.**0.3),
+            (g >= 19.) * (aen < 10.**(0.3 + 0.2*(g - 19.)))
+        )
 
     return psflike
 
@@ -222,7 +237,8 @@ def sub_gaia_edr3(filename, objs=None):
 
     # ADM read the gaia sweep.
     cols = ["PARALLAX", "PARALLAX_IVAR",
-            "PMRA", "PMRA_IVAR", "PMDEC", "PMDEC_IVAR"]
+            "PMRA", "PMRA_IVAR", "PMDEC", "PMDEC_IVAR",
+            "GAIA_ASTROMETRIC_PARAMS_SOLVED"]
     gaiacols = ["EDR3_{}".format(col) for col in cols]
     gswobjs = fitsio.read(gsweepfn, "GAIA_SWEEP", columns=gaiacols)