Update masking to be "all-sky" using Gaia/Tycho/URAT

bin/make_bright_mask

@@ -1,49 +1,65 @@
 #!/usr/bin/env python
 
-from __future__ import print_function, division
-
 import os, sys
 import numpy as np
 
 from desitarget import io
-from desitarget.brightmask import make_bright_source_mask
+from desitarget.brightmask import make_bright_star_mask, get_mask_dir
+from desitarget.tychomatch import get_tycho_nside, get_tycho_dir
+from desitarget.uratmatch import get_urat_dir
+from desitarget.gaiamatch import get_gaia_dir
 
 #import warnings
 #warnings.simplefilter('error')
 
 import multiprocessing
 nproc = multiprocessing.cpu_count() // 2
+nside = get_tycho_nside()
 
 from desiutil.log import get_logger
 log = get_logger()
 
 from argparse import ArgumentParser
 ap = ArgumentParser()
-ap.add_argument("src", help="root directory with tractor/sweeps files from which to build mask")
-ap.add_argument("dest", help='Output file name to which to write mask')
-ap.add_argument("--bands", 
-                help='Bands to use to build the mask (e.g. GRZ)',
-                default="GRZ")
-ap.add_argument("--maglim",
-                help='Magnitude limits for building the mask (e.g. 10,10,10)',
-                default="10,10,10")
-### ap.add_argument('-b', "--bricklist", help='filename with list of bricknames to include')
+ap.add_argument("--maskdir",
+                help='Output directory to which to write mask [defaults to the $MASK_DIR environment variable]')
+ap.add_argument("--maglim", type=float,
+                help="Magnitude limit for building the mask [defaults to 12]",
+                default=12.)
+ap.add_argument("--matchrad", type=float,
+                help="Tycho stars are discarded if they match a Gaia star at this radius",
+                default=1.)
+ap.add_argument("--nside", type=int,
+                help="Write output files in HEALPixels at this nside [{}]".format(nside))
 ap.add_argument("--numproc", type=int,
-    help='number of concurrent processes to use [{}]'.format(nproc),
-    default=nproc)
+                help="number of concurrent processes to use [{}]".format(nproc),
+                default=nproc)
+ap.add_argument("--maskepoch", type=float,
+                help="Stars are shifted (by just proper motion) to this epoch to make the mask [2023.0]",
+                default=2023.0)
+ap.add_argument("--gaiaepoch", type=float,
+                help="Epoch of Gaia observations [2015.5]",
+                default=2015.5)
 
 ns = ap.parse_args()
-infiles = io.list_sweepfiles(ns.src)
-if len(infiles) == 0:
-    infiles = io.list_tractorfiles(ns.src)
-if len(infiles) == 0:
-    log.critical('no sweep or tractor files found')
-    sys.exit(1)
 
-maglim = [ float(ml) for ml in ns.maglim.split(',') ]
+# ADM if the MASK_DIR directory was passed, set it...
+maskdir = ns.maskdir
+if maskdir is None:
+    # ADM ...otherwise retrieve it from the environment variable.
+    maskdir = get_mask_dir()
+
+mask = make_bright_star_mask(
+    maglim=ns.maglim, matchrad=ns.matchrad, numproc=ns.numproc,
+    maskepoch=ns.maskepoch, gaiaepoch=ns.gaiaepoch)
 
-sourcemask = make_bright_source_mask(ns.bands,maglim,numproc=ns.numproc,
-                                   rootdirname=ns.src,outfilename=ns.dest)
+# ADM extra header keywords for the output fits file.
+extra = {k: v for k, v in zip(
+    ["matchrad", "gaiaepoc", "gaiadir", "tychodir", "uratdir"],
+    [ns.matchrad, ns.gaiaepoch, get_gaia_dir(), get_tycho_dir(), get_urat_dir()])}
 
-log.info('wrote a file of {} masks to {}'.format(len(sourcemask), ns.dest))
+# ADM write the mask file(s).
+nmasks, mxdir = io.write_masks(maskdir, mask, maglim=ns.maglim,
+                              maskepoch=ns.maskepoch, nside=nside, extra=extra)
 
+log.info('wrote {} total masks to files in {}'.format(nmasks, mxdir))

bin/select_skies

@@ -8,9 +8,11 @@ from desitarget import io
 from desitarget.io import check_both_set
 from desitarget.geomask import bundle_bricks
 from desitarget.targets import resolve
+from desitarget.brightmask import mask_targets, get_recent_mask_dir
 
 import numpy as np
 import healpy as hp
+from glob import glob
 
 #import warnings
 #warnings.simplefilter('error')
@@ -58,10 +60,26 @@ ap.add_argument("--brickspersec", type=float,
                 help="estimate of bricks completed per second by the (parallelized) code. Used with `bundlebricks` to guess run times (defaults to 1.6)",
                 default=1.6)
 ap.add_argument("--numproc", type=int,
-                help='number of concurrent processes to use [{}]'.format(nproc),
+                help="number of concurrent processes to use [{}]".format(nproc),
                 default=nproc)
+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 sky locations")
+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 build the list of command line arguments as
+# ADM bundlefiles potentially needs to know about them.
+extra = " --numproc {}".format(ns.numproc)
+nsdict = vars(ns)
+for nskey in "nskiespersqdeg", "bands", "apertures", "nomasking", "maskdir":
+    if isinstance(nsdict[nskey], bool):
+        if nsdict[nskey]:
+            extra += " --{}".format(nskey)
+    elif nsdict[nskey] is not None:
+        extra += " --{} {}".format(nskey, nsdict[nskey])
 
 for indir in ns.surveydir, ns.surveydir2:
     if indir is not None:
@@ -106,8 +124,9 @@ if ns.bundlebricks is not None:
     # ADM pixnum only contains unique bricks, need to add duplicates.
     allpixnum = np.concatenate([np.zeros(cnt, dtype=int)+pix
                         for cnt, pix in zip(cnts.astype(int), pixnum)])
-    bundle_bricks(allpixnum, ns.bundlebricks, ns.nside, prefix='skies',
-                  gather=False, surveydirs=drdirs, brickspersec=ns.brickspersec)
+    bundle_bricks(
+        allpixnum, ns.bundlebricks, ns.nside, prefix='skies', extra=extra,
+        gather=False, surveydirs=drdirs, brickspersec=ns.brickspersec)
 else:
     log.info("running on {} processors".format(ns.numproc))
     # ADM formally writing pixelized files requires both the nside
@@ -139,11 +158,20 @@ else:
     else:
         skies = resolved
 
+    # ADM mask the sky locations using a bright star mask.
+    if do_mask:
+        maskdir = get_recent_mask_dir(ns.maskdir)
+        skies = mask_targets(skies, maskdir, nside=ns.nside, pixlist=pixlist)
+
+    # ADM extra header keywords for the output fits file.
+    extra = {k: v for k, v in zip(["masked"],
+                                  [do_mask])}
+
     # ADM this correctly records the apertures in the output file header
     # ADM as well as adding HEALPixel information.
     nskies, outfile = io.write_skies(ns.dest, skies, indir=ns.surveydir,
                                      indir2=ns.surveydir2, nside=nside,
-                                     apertures_arcsec=apertures,
+                                     apertures_arcsec=apertures, extra=extra,
                                      nskiespersqdeg=nskiespersqdeg,
                                      nsidefile=ns.nside, hpxlist=pixlist)
 

bin/supplement_skies

@@ -3,6 +3,7 @@
 from desitarget.skyfibers import supplement_skies, density_of_sky_fibers
 from desitarget import io
 from desitarget.geomask import bundle_bricks, shares_hp, pixarea2nside
+from desitarget.brightmask import mask_targets, get_recent_mask_dir
 
 import os, sys
 import numpy as np
@@ -56,8 +57,14 @@ ap.add_argument("--mingalb", type=float,
 ap.add_argument("--radius", type=float,
                 help="Radius at which to avoid (all) Gaia sources (arcseconds; defaults to [2])",
                 default=2.)
+ap.add_argument("--nomasking", action='store_true',
+                help="Masking occurs by default. If this is set, do NOT use Tycho+Gaia+URAT bright star mask to mask the sky locations")
+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 check the input sky file is in the "official" format.
 official = io.is_sky_dir_official(ns.skydir)
@@ -81,8 +88,8 @@ if ns.bundlefiles is not None:
 gaiadir = ns.gaiadir
 if gaiadir is None:
     # ADM ...otherwise retrieve it from the environment variable.
-    from desitarget.gaiamatch import _get_gaia_dir
-    gaiadir = _get_gaia_dir()
+    from desitarget.gaiamatch import get_gaia_dir
+    gaiadir = get_gaia_dir()
 
 # ADM if needed, determine the minimum density of sky fibers to generate.
 nskiespersqdeg = ns.nskiespersqdeg
@@ -96,8 +103,12 @@ extra = " --numproc {}".format(ns.numproc)
 extra += " --nskiespersqdeg {}".format(nskiespersqdeg)
 extra += " --gaiadir {}".format(gaiadir)
 nsdict = vars(ns)
-for nskey in "mindec", "mingalb", "radius":
-    extra += " --{} {}".format(nskey, nsdict[nskey])
+for nskey in "mindec", "mingalb", "radius", "nomasking", "maskdir":
+    if isinstance(nsdict[nskey], bool):
+        if nsdict[nskey]:
+            extra += " --{}".format(nskey)
+    elif nsdict[nskey] is not None:
+        extra += " --{} {}".format(nskey, nsdict[nskey])
 
 # ADM only proceed if we're not writing a slurm script.
 if ns.bundlefiles is None:
@@ -135,15 +146,21 @@ if ns.bundlefiles is None:
                                   pixlist=pixlist, mindec=ns.mindec,
                                   mingalb=ns.mingalb)
 
+    # ADM mask the supplemental sky locations using a bright star mask.
+    if do_mask:
+        maskdir = get_recent_mask_dir(ns.maskdir)
+        supp_skies = mask_targets(supp_skies, maskdir,
+                                  nside=ns.nside, pixlist=pixlist)
+
     # ADM remove supplemental skies that share HEALPixels with skies.
     nside_resol = pixarea2nside(1./18000)//2
     ii, _ = shares_hp(nside_resol, supp_skies, skies)
     supp_skies = supp_skies[~ii]
     log.info("Removed {} supp skies that matched skies".format(np.sum(ii)))
 
     # ADM extra header keywords for the output fits file.
-    extra = {k: v for k, v in zip(["radius", "mindec", "mingalb"],
-                                  [ns.radius, ns.mindec, ns.mingalb])}
+    extra = {k: v for k, v in zip(["radius", "mindec", "mingalb", "masked"],
+                                  [ns.radius, ns.mindec, ns.mingalb, do_mask])}
 
     nskies, outfile = io.write_skies(ns.dest, supp_skies, supp=True, indir=gaiadir,
                                      nside=nside, nskiespersqdeg=nskiespersqdeg,

doc/changes.rst

@@ -5,6 +5,18 @@ desitarget Change Log
 0.40.1 (unreleased)
 -------------------
 
+* Update masking to be "all-sky" using Gaia/Tycho/URAT [`PR #625`_]:
+    * General desitarget functionality to work with Tycho files.
+    * Deprecate using the sweeps to mask bright objects as this is now
+      being done using MASKBITS from the imaging catalogs.
+    * Functionality to allow masks to be built at different epochs, via
+      careful treatment of Tycho/Gaia/URAT proper motions.
+    * Bright star masks are now explicitly written to a $MASK_DIR.
+    * The radius-magnitude relationship is now a single function.
+    * Refactoring of unit tests to be simpler and have more coverage.
+    * Skies and supplemental skies are now always masked by default.
+    * A lack of backward compatibility, which should be OK as the masking
+      formalism wasn't being extensively used.
 * Functionality for iterations of SV beyond sv1 [`PR #624`_]. Includes:
     * A script to create the necessary files for new iterations of SV.
     * Generalized mask/cuts handling for survey=svX, X being any integer.
@@ -21,6 +33,7 @@ desitarget Change Log
 
 .. _`PR #623`: https://github.com/desihub/desitarget/pull/623
 .. _`PR #624`: https://github.com/desihub/desitarget/pull/624
+.. _`PR #625`: https://github.com/desihub/desitarget/pull/625
 
 0.40.0 (2020-05-26)
 -------------------

etc/desitarget.module

@@ -76,8 +76,11 @@ setenv [string toupper $product] $PRODUCT_DIR
 #
 # Add any non-standard Module code below this point.
 #
-setenv GAIA_DIR /project/projectdirs/desi/target/gaia_dr2
-setenv CMX_DIR /project/projectdirs/desi/target/cmx_files
-setenv SCND_DIR /project/projectdirs/desi/target/secondary
-setenv URAT_DIR /project/projectdirs/desi/target/urat_dr1
-setenv TARG_DIR /project/projectdirs/desi/target/catalogs
+setenv TARG_DIR /global/cfs/cdirs/desi/target/catalogs
+setenv DUST_DIR /global/cfs/cdirs/cosmo/data/dust/v0_1
+setenv CMX_DIR /global/cfs/cdirs/desi/target/cmx_files
+setenv SCND_DIR /global/cfs/cdirs/desi/target/secondary
+setenv MASK_DIR /global/cfs/cdirs/desi/target/masks
+setenv GAIA_DIR /global/cfs/cdirs/desi/target/gaia_dr2
+setenv URAT_DIR /global/cfs/cdirs/desi/target/urat_dr1
+setenv TYCHO_DIR /global/cfs/cdirs/desi/target/tycho_dr2