Merge ab1845d into e03f809

bin/make_imaging_weight_map

@@ -31,7 +31,7 @@ ap.add_argument("--nside", type=int,
                 help='The resolution (HEALPixel nside number) at which to build the map (defaults to 256)',
                 default="256")
 ap.add_argument("--gaialoc",
-                help='Directory of Gaia "chunks" files (to derive stellar density). A FITS file that already contains the Gaia stellar densities at nside can be passed as a speed-up')
+                help='A FITS file that already contains the Gaia stellar densities at nside to speed-up density calculations')
 
 ns = ap.parse_args()
 
@@ -43,12 +43,14 @@ if not os.path.exists(ns.targets):
     log.critical('Input directory does not exist: {}'.format(ns.targets))
     sys.exit(1)
 
-hdr = fitsio.read_header(ns.randoms,"RANDOMS")
-#ADM add HEALPixel information to the header
+hdr = fitsio.read_header(ns.randoms, "RANDOMS")
+#ADM add HEALPixel and gaialoc information to the header
+hdr['GAIALOC'] = ns.gaialoc
 hdr['HPXNSIDE'] = ns.nside
 hdr['HPXNEST'] = True
 
-pixmap = pixmap(ns.randoms, ns.targets, hdr["DENSITY"], nside=ns.nside, gaialoc=ns.gaialoc)
+pixmap, survey = pixmap(ns.randoms, ns.targets, hdr["DENSITY"], nside=ns.nside, gaialoc=ns.gaialoc)
+hdr["SURVEY"] = survey
 
 #ADM write out the map
 log.info('Writing pixel map to {}'.format(ns.dest))

bin/run_target_qa

@@ -7,10 +7,15 @@ import numpy as np
 import fitsio
 
 from desitarget.QA import make_qa_page, _parse_tcnames
+from desitarget.io import read_targets_in_box_header
 
 from desiutil.log import get_logger
 log = get_logger()
 
+import multiprocessing
+# ADM the rather large denominator guards against memory issues.
+nproc = multiprocessing.cpu_count() // 8
+
 from argparse import ArgumentParser
 ap = ArgumentParser()
 ap.add_argument("src", 
@@ -29,14 +34,30 @@ ap.add_argument('-t','--tcnames', default=None,
                 help="Names of specific target classes for which to make the QA pages. (e.g. QSO,LRG,ALL)")
 ap.add_argument('--nosystematics', action='store_true', 
                 help="Don't make systematics plots on the front page (perhaps because the weightmapfile was not passed)")
+ap.add_argument("--numproc", type=int,
+                help='number of concurrent processes to use when generating plots [defaults to {}]'.format(nproc),
+                default=nproc)
 
 ns = ap.parse_args()
 
-#ADM fail if systematics were requested but no systematics map was passed
+# ADM fail if systematics were requested but no systematics map was passed
 if ns.weightmapfile is None and ns.nosystematics is False:
     log.error("The weightmap file was not passed so systematics cannot be tracked. Try again sending --nosystematics.")
     raise IOError
 
+# ADM check that the passed pixweight file corresponds to the passed target file.
+if ns.weightmapfile is not None and not ns.mocks:
+    dsurv = read_targets_in_box_header(ns.src)["SURVEY"]
+    try:
+        psurv = fitsio.read_header(ns.weightmapfile, "PIXWEIGHTS")["SURVEY"]
+        if dsurv != psurv and 'cmx' not in dsurv:
+            msg = 'target file is type {} but pixweight file is type {}!!!'.format(dsurv, psurv)
+            log.critical(msg)
+            raise ValueError(msg)
+    # ADM we'll flag an error for early versions of files, which is fine.
+    except KeyError:
+        pass
+
 tcnames = ns.tcnames
 if tcnames is not None:
     tcnames = _parse_tcnames(tcstring=ns.tcnames, add_all=True)
@@ -50,7 +71,8 @@ if ns.nside:
 else:
     max_bin_area = 1.0
 
-make_qa_page(ns.src, qadir=ns.dest,
+log.info("generating plots on {} processors".format(ns.numproc))
+make_qa_page(ns.src, qadir=ns.dest, numproc=ns.numproc,
              mocks=ns.mocks, makeplots=makeplots, max_bin_area=max_bin_area, 
              imaging_map_file=ns.weightmapfile, tcnames=tcnames, systematics=systematics)
 log.info('Targeting QA pages and plots written to {}'.format(ns.dest))

bin/select_gfas

@@ -14,7 +14,7 @@ log = get_logger()
 
 import multiprocessing
 nproc = multiprocessing.cpu_count() // 2
-nside = 64 #ADM default HEALPix Nside used throughout desitarget
+nside = io.desitarget_nside()   # ADM default HEALPix Nside used throughout desitarget
 
 from argparse import ArgumentParser
 ap = ArgumentParser(description='Generates a file of GFA (Guide/Focus/Alignment) targets via matching to Gaia')
@@ -51,17 +51,14 @@ if ns.cmx:
     survey = 'cmx'
 
 infiles = io.list_sweepfiles(ns.surveydir)
-indir = ns.surveydir
 if ns.surveydir2 is not None:
     infiles2 = io.list_sweepfiles(ns.surveydir2)
     infiles += infiles2
-    indir = "{} {}".format(ns.surveydir, ns.surveydir2)
 if len(infiles) == 0:
     infiles = io.list_tractorfiles(ns.surveydir)
     if ns.surveydir2 is not None:
         infiles2 = io.list_tractorfiles(ns.surveydir2)
         infiles += infiles2
-        indir = "{} {}".format(ns.surveydir, ns.surveydir2)
 if len(infiles) == 0:
     log.critical('no sweep or tractor files found')
     sys.exit(1)
@@ -79,6 +76,6 @@ gfas = gfas[decgood]
 bbad = is_in_gal_box(gfas, [0., 360., -ns.mingalb, ns.mingalb])
 gfas = gfas[~bbad]
 
-io.write_gfas(ns.dest, gfas, indir=indir, nside=nside, survey=survey)
+io.write_gfas(ns.dest, gfas, indir=ns.surveydir, indir2=ns.surveydir2, nside=nside, survey=survey)
 
 log.info('{} GFAs written to {}...t = {:.1f} mins'.format(len(gfas), ns.dest, (time()-t0)/60.))

bin/select_randoms

@@ -53,10 +53,12 @@ ap.add_argument("--dustdir",
                 default=None)
 ap.add_argument("--noresolve", action='store_true',
                 help="Do NOT resolve into northern randoms in northern regions and southern randoms in southern regions")
+ap.add_argument("--aprad", type=float,
+                help="Radius of aperture in arcsec in which to generated sky background flux levels (defaults to 0.75; the DESI fiber radius)",
+                default=0.75)
 
 
 ns = ap.parse_args()
-
 # ADM parse the list of HEALPixels in which to run.
 pixlist = ns.healpixels
 if pixlist is not None:
@@ -90,12 +92,12 @@ if not 'dr5' in ns.surveydir and not 'dr6' in ns.surveydir:
             hdr[record] = rmhdr['_record_map'][record]
 
 randoms = select_randoms(ns.surveydir, density=ns.density, numproc=ns.numproc,
-                         nside=ns.nside, pixlist=pixlist,
+                         nside=ns.nside, pixlist=pixlist, aprad=ns.aprad,
                          bundlebricks=ns.bundlebricks, brickspersec=ns.brickspersec,
                          dustdir=ns.dustdir, resolverands=not(ns.noresolve))
 
 if ns.bundlebricks is None:
-    io.write_randoms(ns.dest, randoms, indir=ns.surveydir, 
+    io.write_randoms(ns.dest, randoms, indir=ns.surveydir, aprad=ns.aprad,
                      hdr=hdr, nside=nside, density=ns.density, resolve=not(ns.noresolve))
     log.info('wrote file of {} randoms to {}...t = {:.1f}s'
              .format(len(randoms), ns.dest,time()-start))

bin/select_skies

@@ -3,10 +3,13 @@
 from __future__ import print_function, division
 
 from desitarget.skyutilities.legacypipe.util import LegacySurveyData
-from desitarget.skyfibers import select_skies, density_of_sky_fibers
-from desitarget.io import write_skies
+from desitarget.skyfibers import select_skies, density_of_sky_fibers, get_brick_info
+from desitarget import io
+from desitarget.geomask import bundle_bricks
+from desitarget.targets import resolve
 
 import numpy as np
+import healpy as hp
 
 #import warnings
 #warnings.simplefilter('error')
@@ -15,7 +18,7 @@ import multiprocessing
 nproc = multiprocessing.cpu_count() // 2
 # ADM default HEALPix Nside used throughout desitarget.
 # ADM don't confuse this with the ns.nside parallelization input that is parsed below!!!
-nside = 64 
+nside = io.desitarget_nside()
 
 from desiutil.log import get_logger
 log = get_logger()
@@ -26,6 +29,9 @@ ap.add_argument("surveydir",
                 help="Base directory for a Legacy Surveys Data Release (e.g. '/global/project/projectdirs/cosmo/data/legacysurvey/dr6/' at NERSC")
 ap.add_argument("dest", 
                 help="Output sky targets file (e.g. /project/projectdirs/desi/target/catalogs/skies-dr4-0.20.0.fits)")
+ap.add_argument('-s2', "--surveydir2",
+                help='Additional Legacy Surveys directory (useful for combining, e.g., DR8 into one file of sky locations)',
+                default=None)
 ap.add_argument("--nskiespersqdeg", 
                 help="Number of sky locations to generate per sq. deg. (don't pass to read the default from desimodel.io with a 4x margin)",
                 default=None)
@@ -73,18 +79,42 @@ nskiesfloat = area*nskiespersqdeg
 nskies = (np.sqrt(nskiesfloat).astype('int16') + 1)**2
 log.info('Generating {} sky positions in each brick'.format(nskies))
 
-survey = LegacySurveyData(survey_dir=ns.surveydir)
+surveys = [LegacySurveyData(survey_dir=ns.surveydir)]
+if ns.surveydir2 is not None:
+    surveys.append(LegacySurveyData(survey_dir=ns.surveydir2))
+
+# ADM if bundlebricks is set, grab the HEALPixel number for each brick.
+if ns.bundlebricks is not None:
+    drdirs = [survey.survey_dir for survey in surveys]
+    brickdict = get_brick_info(drdirs, counts=True)
+    bra, bdec, _, _, _, _, cnts = np.vstack(brickdict.values()).T
+    theta, phi = np.radians(90-bdec), np.radians(bra)
+    pixnum = hp.ang2pix(ns.nside, theta, phi, nest=True)
+    # 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',
+                  surveydirs=drdirs, brickspersec=ns.brickspersec)
+else:
+    # ADM run the main sky selection code over the passed surveys.
+    skies = []
+    for survey in surveys:
+        skies.append(select_skies(
+            survey, numproc=ns.numproc, nskiespersqdeg=nskiespersqdeg,
+            bands=bands, apertures_arcsec=apertures,
+            nside=ns.nside, pixlist=pixlist, writebricks=ns.writebricks)
+        )
+    # ADM redact empty output (where there were no bricks in a survey).
+    skies = np.array(skies) 
+    ii = [sk is not None for sk in skies]
+    skies = np.concatenate(skies[ii])
 
-skies = select_skies(survey, numproc=ns.numproc, nskiespersqdeg=nskiespersqdeg, 
-                     bands=bands, apertures_arcsec=apertures,
-                     nside=ns.nside, pixlist=pixlist, writebricks=ns.writebricks,
-                     bundlebricks=ns.bundlebricks, brickspersec=ns.brickspersec)
+    # ADM resolve any duplicates between imaging data releases. 
+    skies = resolve(skies)
 
-# ADM if skies is None then the bundlebricks informational option was sent.
-if skies is not None:
     # ADM this correctly records the apertures in the output file header
     # ADM as well as adding HEALPixel information.
-    write_skies(ns.dest, skies, indir=ns.surveydir, nside=nside,
-                apertures_arcsec=apertures, nskiespersqdeg=nskiespersqdeg)
+    io.write_skies(ns.dest, skies, indir=ns.surveydir, indir2=ns.surveydir2,
+                   nside=nside, apertures_arcsec=apertures, nskiespersqdeg=nskiespersqdeg)
 
     log.info('{} skies written to {}'.format(len(skies), ns.dest))

bin/select_sv_targets

@@ -26,10 +26,13 @@ log = get_logger()
 
 from argparse import ArgumentParser
 ap = ArgumentParser(description='Generates DESI SV target bits from Legacy Surveys sweeps or tractor files')
-ap.add_argument("src", 
+ap.add_argument("sweepdir",
                 help="Tractor/sweeps file or root directory with tractor/sweeps files")
-ap.add_argument("dest", 
+ap.add_argument("dest",
                 help="Output target selection file")
+ap.add_argument('-s2', "--sweepdir2",
+                help='Additional Tractor/sweeps file or directory (useful for combining, e.g., DR8 into one file of targets)',
+                default=None)
 ap.add_argument('-m', "--mask", 
                 help="If sent then mask the targets, the name of the mask file should be supplied")
 ap.add_argument("--numproc", type=int,
@@ -57,11 +60,21 @@ ap.add_argument('--radecbox',
 ap.add_argument('--radecrad',
                 help="Only return targets in an RA/Dec circle/cap denoted by 'centerRA,centerDec,radius' in degrees (e.g. '140,150,0.5')",
                 default=None)
+ap.add_argument("--noresolve", action='store_true',
+                help="Do NOT resolve into northern targets in northern regions and southern targets in southern regions")
+
 
 ns = ap.parse_args()
-infiles = io.list_sweepfiles(ns.src)
+
+infiles = io.list_sweepfiles(ns.sweepdir)
+if ns.sweepdir2 is not None:
+    infiles2 = io.list_sweepfiles(ns.sweepdir2)
+    infiles += infiles2
 if len(infiles) == 0:
-    infiles = io.list_tractorfiles(ns.src)
+    infiles = io.list_tractorfiles(ns.sweepdir)
+    if ns.sweepdir2 is not None:
+        infiles2 = io.list_tractorfiles(ns.sweepdir2)
+        infiles += infiles2
 if len(infiles) == 0:
     log.critical('no sweep or tractor files found')
     sys.exit(1)
@@ -108,12 +121,12 @@ targets = select_targets(infiles, numproc=ns.numproc,
                          nside=ns.nside, pixlist=pixlist,
                          bundlefiles=ns.bundlefiles, filespersec=ns.filespersec,
                          radecbox=inlists[0], radecrad=inlists[1],
-                         tcnames=tcnames, survey=survey)
+                         tcnames=tcnames, survey=survey, resolvetargs=not(ns.noresolve))
 if ns.mask:
     targets = mask_targets(targets, inmaskfile=ns.mask, nside=nside)
 
 if ns.bundlefiles is None:
-    io.write_targets(ns.dest, targets, indir=ns.src, survey=survey,
-                     nsidefile=ns.nside, hpxlist=pixlist,
-                     qso_selection=survey, nside=nside)
+    io.write_targets(ns.dest, targets, indir=ns.sweepdir, indir2=ns.sweepdir2,
+                     survey=survey, nsidefile=ns.nside, hpxlist=pixlist,
+                     qso_selection=survey, nside=nside, resolve=not(ns.noresolve))
     log.info('{} targets written to {}...t={:.1f}s'.format(len(targets), ns.dest, time()-start))

bin/select_targets

@@ -77,17 +77,14 @@ ap.add_argument("--noresolve", action='store_true',
 ns = ap.parse_args()
 
 infiles = io.list_sweepfiles(ns.sweepdir)
-indir = ns.sweepdir
 if ns.sweepdir2 is not None:
     infiles2 = io.list_sweepfiles(ns.sweepdir2)
     infiles += infiles2
-    indir = "{} {}".format(ns.sweepdir, ns.sweepdir2)
 if len(infiles) == 0:
     infiles = io.list_tractorfiles(ns.sweepdir)
     if ns.sweepdir2 is not None:
         infiles2 = io.list_tractorfiles(ns.sweepdir2)
         infiles += infiles2
-        indir = "{} {}".format(ns.sweepdir, ns.sweepdir2)
 if len(infiles) == 0:
     log.critical('no sweep or tractor files found')
     sys.exit(1)
@@ -125,7 +122,7 @@ else:
         targets = mask_targets(targets, inmaskfile=ns.mask, nside=nside)
 
     if ns.bundlefiles is None:
-        io.write_targets(ns.dest, targets, resolve=not(ns.noresolve),
-                         indir=indir, survey="main", nsidefile=ns.nside, hpxlist=pixlist,
+        io.write_targets(ns.dest, targets, resolve=not(ns.noresolve), indir=ns.sweepdir,
+                         indir2=ns.sweepdir2, survey="main", nsidefile=ns.nside, hpxlist=pixlist,
                          qso_selection=ns.qsoselection, sandboxcuts=ns.sandbox, nside=nside)
         log.info('{} targets written to {}...t={:.1f}s'.format(len(targets), ns.dest, time()-start))

doc/changes.rst

@@ -5,6 +5,19 @@ desitarget Change Log
 0.29.2 (unreleased)
 -------------------
 
+* Further updates and enhancements for DR8 [`PR #490`_]. Includes:
+    * Resolve sky locations and SV targets in North/South regions.
+    * Update sky and SV slurming for DR8-style input (two directories).
+    * Write both of two input directories to output file headers.
+    * Parallelize plot production to speed-up QA by factors of 8.
+    * Add ``PSFSIZE`` to randoms, pixweight maps and QA plots.
+    * QA and pixweight maps work fully for SV-style files and bits.
+    * Pixweight code can now take HEALpixel-split targets as input.
+    * Add aperture-photometered background flux to randoms catalogs.
+    * Additional unit test module (:func:`test.test_geomask`).
+    * Deprecate `make_hpx_density_file`; use `make_imaging_weight_map`.
+    * :func:`io.read_targets_in_a_box` can now read headers.
+    * Update unit test data for new DR8 columns and functionality.
 * Update QSO targeting algorithms for DR8 [`PR #489`_]. Includes:
     * Update baseline quasar selection for the main survey.
     * Update QSO bits and selection algorithms for SV.
@@ -18,6 +31,7 @@ desitarget Change Log
 .. _`PR #484`: https://github.com/desihub/desitarget/pull/484
 .. _`PR #488`: https://github.com/desihub/desitarget/pull/488
 .. _`PR #489`: https://github.com/desihub/desitarget/pull/489
+.. _`PR #490`: https://github.com/desihub/desitarget/pull/490
 
 0.29.1 (2019-03-26)
 -------------------