From ff4d9c10b97a701bd7254595d53ea7e1a78e9660 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 2 Apr 2019 12:50:12 -0700 Subject: [PATCH 01/31] add resolve and multiple input sweeps options to sv_select_targets --- bin/select_sv_targets | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/bin/select_sv_targets b/bin/select_sv_targets index 1f106aa95..581734100 100755 --- a/bin/select_sv_targets +++ b/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,24 @@ 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) +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.src) + 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) @@ -108,7 +124,7 @@ 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) From 7c155bcee6ab94dfbb0b5e3a9c529c0985d06fc3 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 2 Apr 2019 12:56:33 -0700 Subject: [PATCH 02/31] code style clean-up --- py/desitarget/cuts.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/py/desitarget/cuts.py b/py/desitarget/cuts.py index 0ce2a992a..f099c402d 100644 --- a/py/desitarget/cuts.py +++ b/py/desitarget/cuts.py @@ -1204,17 +1204,15 @@ def isQSO_randomforest(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=N pcut[tmp_r_Reduced > 21.5] = 0.8125 - 0.15 * (tmp_r_Reduced[tmp_r_Reduced > 21.5] - 21.5) pcut[tmp_r_Reduced > 22.3] = 0.6925 - 0.70 * (tmp_r_Reduced[tmp_r_Reduced > 22.3] - 22.3) pcut_HighZ = np.where(tmp_r_Reduced > 20.5, - 0.55 - (tmp_r_Reduced - 20.5) * 0.025, 0.55) + 0.55 - (tmp_r_Reduced - 20.5) * 0.025, 0.55) else: pcut = np.where(tmp_r_Reduced > 20.8, 0.90 - (tmp_r_Reduced - 20.8) * 0.025, 0.90) pcut[tmp_r_Reduced > 21.5] = 0.8825 - 0.15 * (tmp_r_Reduced[tmp_r_Reduced > 21.5] - 21.5) pcut[tmp_r_Reduced > 22.3] = 0.7625 - 0.70 * (tmp_r_Reduced[tmp_r_Reduced > 22.3] - 22.3) pcut_HighZ = np.where(tmp_r_Reduced > 20.5, - 0.70 - (tmp_r_Reduced - 20.5) * 0.025, 0.70) + 0.70 - (tmp_r_Reduced - 20.5) * 0.025, 0.70) - - # Add rf proba test result to "qso" mask qso[colorsReducedIndex[tmpReleaseOK]] = \ (tmp_rf_proba >= pcut) | (tmp_rf_HighZ_proba >= pcut_HighZ) From bdaf977c36c76c33cb5141c6337989af564c9b59 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 2 Apr 2019 13:45:45 -0700 Subject: [PATCH 03/31] finalizing resolve for SV targets --- bin/select_sv_targets | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bin/select_sv_targets b/bin/select_sv_targets index 581734100..84d0ae068 100755 --- a/bin/select_sv_targets +++ b/bin/select_sv_targets @@ -129,7 +129,7 @@ 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, + io.write_targets(ns.dest, targets, indir=indir, survey=survey, nsidefile=ns.nside, hpxlist=pixlist, qso_selection=survey, nside=nside) log.info('{} targets written to {}...t={:.1f}s'.format(len(targets), ns.dest, time()-start)) From a6db7612c0bb4cb03e582eed392a2469fb8de6b3 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 2 Apr 2019 13:57:22 -0700 Subject: [PATCH 04/31] add full directory and resolve information to headers of output files --- bin/select_sv_targets | 10 +++++----- bin/select_targets | 8 ++++---- py/desitarget/io.py | 14 ++++++++------ 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/bin/select_sv_targets b/bin/select_sv_targets index 84d0ae068..703145a52 100755 --- a/bin/select_sv_targets +++ b/bin/select_sv_targets @@ -71,13 +71,13 @@ indir = ns.sweepdir if ns.sweepdir2 is not None: infiles2 = io.list_sweepfiles(ns.sweepdir2) infiles += infiles2 - indir = "{} {}".format(ns.sweepdir, ns.sweepdir2) + indir2 = 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) + indir2 = ns.sweepdir2 if len(infiles) == 0: log.critical('no sweep or tractor files found') sys.exit(1) @@ -129,7 +129,7 @@ if ns.mask: targets = mask_targets(targets, inmaskfile=ns.mask, nside=nside) if ns.bundlefiles is None: - io.write_targets(ns.dest, targets, indir=indir, survey=survey, - nsidefile=ns.nside, hpxlist=pixlist, - qso_selection=survey, nside=nside) + io.write_targets(ns.dest, targets, indir=indir, indir2=indir2, + 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)) diff --git a/bin/select_targets b/bin/select_targets index f2100c98c..1a7f64106 100755 --- a/bin/select_targets +++ b/bin/select_targets @@ -81,13 +81,13 @@ indir = ns.sweepdir if ns.sweepdir2 is not None: infiles2 = io.list_sweepfiles(ns.sweepdir2) infiles += infiles2 - indir = "{} {}".format(ns.sweepdir, ns.sweepdir2) + indir2 = 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) + indir2 = ns.sweepdir2 if len(infiles) == 0: log.critical('no sweep or tractor files found') sys.exit(1) @@ -125,7 +125,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=indir, + indir2=indir2, 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)) diff --git a/py/desitarget/io.py b/py/desitarget/io.py index 603fbd69d..6a7365650 100644 --- a/py/desitarget/io.py +++ b/py/desitarget/io.py @@ -475,9 +475,9 @@ def release_to_photsys(release): return r2p[release] -def write_targets(filename, data, indir=None, qso_selection=None, - sandboxcuts=False, nside=None, survey="?", - nsidefile=None, hpxlist=None, resolve=True): +def write_targets(filename, data, indir=None, indir2=None, + qso_selection=None, sandboxcuts=False, nside=None, + survey="?", nsidefile=None, hpxlist=None, resolve=True): """Write a target catalogue. Parameters @@ -486,9 +486,9 @@ def write_targets(filename, data, indir=None, qso_selection=None, output target selection file. data : :class:`~numpy.ndarray` numpy structured array of targets to save. - indir, qso_selection : :class:`str`, optional, default to `None` - If passed, note these as the input directory and - quasar selection method in the output file header. + indir, indir2, qso_selection : :class:`str`, optional, default to `None` + If passed, note these as the input directory, an additional input + directory, and the QSO selection method in the output file header. sandboxcuts : :class:`bool`, optional, defaults to ``False`` Written to the output file header as `sandboxcuts`. nside : :class:`int`, optional, defaults to `None` @@ -527,6 +527,8 @@ def write_targets(filename, data, indir=None, qso_selection=None, if indir is not None: depend.setdep(hdr, 'tractor-files', indir) + if indir2 is not None: + depend.setdep(hdr, 'tractor-files-2', indir2) if qso_selection is None: log.warning('qso_selection method not specified for output file') From 49230da3d918cb877b9d0dc6c78d8abe24229cac Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Tue, 2 Apr 2019 15:03:54 -0700 Subject: [PATCH 05/31] update slurm script production for select_targets, select_randoms and select_sv_targets to allow 2 input data directories for DR8 --- py/desitarget/cuts.py | 4 +++- py/desitarget/geomask.py | 26 +++++++++++++++++++------- py/desitarget/randoms.py | 2 +- 3 files changed, 23 insertions(+), 9 deletions(-) diff --git a/py/desitarget/cuts.py b/py/desitarget/cuts.py index f099c402d..fb49a51ed 100644 --- a/py/desitarget/cuts.py +++ b/py/desitarget/cuts.py @@ -2256,9 +2256,11 @@ def select_targets(infiles, numproc=4, qso_selection='randomforest', prefix = "targets" if survey != "main": prefix = "{}_targets".format(survey) + # ADM determine if one or two input directories were passed. + surveydirs = list(set([os.path.dirname(fn) for fn in infiles])) bundle_bricks(pixnum, bundlefiles, nside, brickspersec=filespersec, gather=False, - prefix=prefix, surveydir=os.path.dirname(infiles[0])) + prefix=prefix, surveydirs=surveydirs) return # ADM restrict to only input files in a set of HEALPixels, if requested. diff --git a/py/desitarget/geomask.py b/py/desitarget/geomask.py index d077ca319..a6aa6f8c7 100644 --- a/py/desitarget/geomask.py +++ b/py/desitarget/geomask.py @@ -604,8 +604,8 @@ def circle_boundaries(RAcens, DECcens, r, nloc): return np.hstack(ras), np.hstack(decs) -def bundle_bricks(pixnum, maxpernode, nside, brickspersec=1., prefix='targets', gather=True, - surveydir="/global/project/projectdirs/cosmo/data/legacysurvey/dr6"): +def bundle_bricks(pixnum, maxpernode, nside, brickspersec=1., prefix='targets', + gather=True, surveydirs=None): """Determine the optimal packing for bricks collected by HEALpixel integer. Parameters @@ -629,9 +629,11 @@ def bundle_bricks(pixnum, maxpernode, nside, brickspersec=1., prefix='targets', gather : :class:`bool`, optional, defaults to ``True`` If ``True`` then provide a final command for combining all of the HEALPix-split files into one large file. If ``False``, comment out that command. - surveydir : :class:`str`, optional, defaults to the DR6 directory at NERSC - The root directory pointing to a Data Release from the Legacy Surveys, - (e.g. "/global/project/projectdirs/cosmo/data/legacysurvey/dr6"). + surveydirs : :class:`list` + Root directories for a Legacy Surveys Data Release. The first element of the + list is interpreted as the main directory. IF the list is of length two + then the second directory is supplied as "-s2" in the output script. + (e.g. ["/global/project/projectdirs/cosmo/data/legacysurvey/dr6"]). Returns ------- @@ -643,6 +645,12 @@ def bundle_bricks(pixnum, maxpernode, nside, brickspersec=1., prefix='targets', ----- h/t https://stackoverflow.com/questions/7392143/python-implementations-of-packing-algorithm """ + # ADM interpret the passed directories. + surveydir = surveydirs[0] + surveydir2 = None + if len(surveydirs) == 2: + surveydir2 = surveydirs[1] + # ADM the number of pixels (numpix) in each pixel (pix) pix, numpix = np.unique(pixnum, return_counts=True) @@ -738,6 +746,10 @@ def bundle_bricks(pixnum, maxpernode, nside, brickspersec=1., prefix='targets', if prefix[0:2] == "sv": prefix2 = "sv_targets" + s2 = "" + if surveydir2 is not None: + s2 = "-s2 {}".format(surveydir2) + outfiles = [] for bin in bins: num = np.array(bin)[:, 0] @@ -750,8 +762,8 @@ def bundle_bricks(pixnum, maxpernode, nside, brickspersec=1., prefix='targets', # ADM the replace is to handle inputs that look like "sv1_targets". outfile = "$CSCRATCH/{}-dr{}-hp-{}.fits".format(prefix.replace("_", "-"), dr, strgoodpix) outfiles.append(outfile) - print("srun -N 1 select_{} {} {} --numproc 32 --nside {} --healpixels {} &" - .format(prefix2, surveydir, outfile, nside, strgoodpix)) + print("srun -N 1 select_{} {} {} {} --numproc 32 --nside {} --healpixels {} &" + .format(prefix2, surveydir, outfile, s2, nside, strgoodpix)) print("wait") print("") print("{}gather_targets '{}' $CSCRATCH/{}-dr{}.fits {}" diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index e8ced0fd3..a5d979be2 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -944,7 +944,7 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None, # ADM if the bundlebricks option was sent, call the packing code. if bundlebricks is not None: bundle_bricks(pixnum, bundlebricks, nside, brickspersec=brickspersec, - prefix='randoms', surveydir=drdir) + prefix='randoms', surveydirs=[drdir]) return # ADM restrict to only bricks in a set of HEALPixels, if requested. From 454860fe95f35a9a8d507e6efb6fc6bf35bfbd2f Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Wed, 3 Apr 2019 21:24:01 -0700 Subject: [PATCH 06/31] update skyfibers.py, select_skies and sky slurm script production to allow 2 input data directories --- bin/select_skies | 50 +++++++++++---- py/desitarget/randoms.py | 46 +++++--------- py/desitarget/skyfibers.py | 122 +++++++++++++++++++++++-------------- 3 files changed, 128 insertions(+), 90 deletions(-) diff --git a/bin/select_skies b/bin/select_skies index 225334d90..c75836747 100755 --- a/bin/select_skies +++ b/bin/select_skies @@ -3,10 +3,12 @@ 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 import numpy as np +import healpy as hp #import warnings #warnings.simplefilter('error') @@ -15,7 +17,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 +28,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 +78,39 @@ 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)) -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 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]) -# 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, nside=nside, + apertures_arcsec=apertures, nskiespersqdeg=nskiespersqdeg) log.info('{} skies written to {}'.format(len(skies), ns.dest)) diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index a5d979be2..98550d6c7 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -19,6 +19,7 @@ from desitarget.geomask import bundle_bricks, box_area from desitarget.targetmask import desi_mask, bgs_mask, mws_mask from desitarget.targets import resolve +from desitarget.skyfibers import get_brick_info # ADM the parallelization script from desitarget.internal import sharedmem @@ -907,43 +908,26 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None, dr7dir + 'coadd/111/1116p210/legacysurvey-1116p210-maskbits.fits.gz' EBV: E(B-V) at this location from the SFD dust maps """ - # ADM read in the survey bricks file, which lists the bricks of interest for this DR. - # ADM if this is pre-or-post-DR8 we need to find the correct directory or directories. + # ADM grab brick information for this data release. Depending on whether this + # ADM is pre-or-post-DR8 we need to find the correct directory or directories. drdirs = _pre_or_post_dr8(drdir) - bricknames = [] - for dd in drdirs: - sbfile = glob(dd+'/*bricks-dr*') - if len(sbfile) > 0: - sbfile = sbfile[0] - hdu = fits.open(sbfile) - bricknames.append(hdu[1].data['BRICKNAME']) - else: - # ADM hack for test bricks where we didn't generate the bricks file. - from desitarget.io import brickname_from_filename - fns = glob(os.path.join(dd, 'tractor', '*', '*fits')) - bricknames.append([brickname_from_filename(fn) for fn in fns]) - # ADM don't count bricks twice. - bricknames = np.unique(np.concatenate(bricknames)) - # ADM initialize the bricks class, retrieve the brick information look-up table, - # ADM turn it into a fast look-up dictionary. - from desiutil import brick - bricktable = brick.Bricks(bricksize=0.25).to_table() - brickdict = {} - for b in bricktable: - brickdict[b["BRICKNAME"]] = [b["RA"], b["DEC"], - b["RA1"], b["RA2"], b["DEC1"], b["DEC2"]] + brickdict = get_brick_info(drdirs, counts=True) + # ADM this is just the UNIQUE brick names across all surveys. + bricknames = np.array(list(brickdict.keys())) # ADM if the pixlist or bundlebricks option was sent, we'll need the HEALPixel # ADM information for each brick. if pixlist is not None or bundlebricks is not None: - bra, bdec, bramin, bramax, bdecmin, bdecmax = np.vstack( - [np.array(brickdict[bn]) for bn in bricknames]).T + bra, bdec, _, _, _, _, cnts = np.vstack(brickdict.values()).T theta, phi = np.radians(90-bdec), np.radians(bra) pixnum = hp.ang2pix(nside, theta, phi, nest=True) # ADM if the bundlebricks option was sent, call the packing code. if bundlebricks is not None: - bundle_bricks(pixnum, bundlebricks, nside, brickspersec=brickspersec, + # 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, bundlebricks, nside, brickspersec=brickspersec, prefix='randoms', surveydirs=[drdir]) return @@ -952,9 +936,9 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None, # ADM if an integer was passed, turn it into a list. if isinstance(pixlist, int): pixlist = [pixlist] - wbricks = np.where([pix in pixlist for pix in pixnum])[0] - bricknames = bricknames[wbricks] - if len(wbricks) == 0: + ii = [pix in pixlist for pix in pixnum] + bricknames = bricknames[ii] + if len(bricknames) == 0: log.warning('ZERO bricks in passed pixel list!!!') log.info("Processing bricks in (nside={}, pixel numbers={}) HEALPixels" .format(nside, pixlist)) @@ -971,7 +955,7 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None, def _get_quantities(brickname): '''wrapper on nobs_positions_in_a_brick_from_edges() given a brick name''' # ADM retrieve the edges for the brick that we're working on - bra, bdec, bramin, bramax, bdecmin, bdecmax = brickdict[brickname] + bra, bdec, bramin, bramax, bdecmin, bdecmax, _ = brickdict[brickname] # ADM populate the brick with random points, and retrieve the quantities # ADM of interest at those points. diff --git a/py/desitarget/skyfibers.py b/py/desitarget/skyfibers.py index a27bb97cb..e60f70f81 100644 --- a/py/desitarget/skyfibers.py +++ b/py/desitarget/skyfibers.py @@ -5,7 +5,7 @@ desitarget.skyfibers ==================== -Module dealing with the assignation of sky fibers at the pixel-level for target selection +Module to assign sky fibers at the pixel-level for target selection """ import os import sys @@ -15,6 +15,7 @@ from time import time import photutils import healpy as hp +from glob import glob from scipy.ndimage.morphology import binary_dilation, binary_erosion from scipy.ndimage.measurements import label, find_objects, center_of_mass from scipy.ndimage.filters import gaussian_filter @@ -23,9 +24,9 @@ from desitarget.skyutilities.astrometry.fits import fits_table from desitarget.skyutilities.legacypipe.util import find_unique_pixels -from desitarget.geomask import bundle_bricks from desitarget.targetmask import desi_mask, targetid_mask from desitarget.targets import encode_targetid, finalize +from desitarget.io import brickname_from_filename # ADM the parallelization script from desitarget.internal import sharedmem @@ -48,6 +49,64 @@ ('OBSCONDITIONS', '>i4')]) +def get_brick_info(drdirs, counts=False): + """Retrieve brick names and coordinates from Legacy Surveys directories. + + Parameters + ---------- + drdirs : :class:`list` + A list of strings, each of which corresponds to a directory pointing + to a Data Release from the Legacy Surveys. Can be of length one. + e.g. ['/global/project/projectdirs/cosmo/data/legacysurvey/dr7']. + counts : :class:`bool`, optional, defaults to ``False`` + If ``True`` also return a count of the number of times each brick + appears ([RAcen, DECcen, RAmin, RAmax, DECmin, DECmax, CNT]). + + Returns + ------- + :class:`dict` + UNIQUE bricks covered by the Data Release(s). Keys are brick names + and values are a list of the brick center and the brick corners + ([RAcen, DECcen, RAmin, RAmax, DECmin, DECmax]). + + Notes + ----- + - Tries a few different ways in case the survey bricks files have + not ywt been created. + """ + bricknames = [] + for dd in drdirs: + # ADM in the simplest case, read in the survey bricks file, which lists + # ADM the bricks of interest for this DR. + sbfile = glob(dd+'/*bricks-dr*') + if len(sbfile) > 0: + brickinfo = fitsio.read(sbfile[0]) + # ADM fitsio reads things in as bytes, so convert to unicode. + bricknames.append(brickinfo['brickname'].astype('U')) + else: + # ADM hack for test bricks where we didn't generate the bricks file. + fns = glob(os.path.join(dd, 'tractor', '*', '*fits')) + bricknames.append([brickname_from_filename(fn) for fn in fns]) + + # ADM don't count bricks twice. + bricknames, cnts = np.unique(np.concatenate(bricknames), return_counts=True) + + # ADM initialize the bricks class, retrieve the brick information look-up + # ADM table and turn it into a fast look-up dictionary. + from desiutil import brick + bricktable = brick.Bricks(bricksize=0.25).to_table() + brickdict = {} + for b in bricktable: + brickdict[b["BRICKNAME"]] = [b["RA"], b["DEC"], + b["RA1"], b["RA2"], + b["DEC1"], b["DEC2"]] + + # ADM only return the subset of the dictionary with bricks in the DR. + if counts: + return {bn: brickdict[bn] + [cnt] for bn, cnt in zip(bricknames, cnts)} + return {bn: brickdict[bn] for bn in bricknames} + + def density_of_sky_fibers(margin=1.5): """Use positioner patrol size to determine sky fiber density for DESI. @@ -649,8 +708,7 @@ def plot_good_bad_skies(survey, brickname, skies, def select_skies(survey, numproc=16, nskiespersqdeg=None, bands=['g', 'r', 'z'], - apertures_arcsec=[0.75], nside=2, pixlist=None, - writebricks=False, bundlebricks=None, brickspersec=1.6): + apertures_arcsec=[0.75], nside=2, pixlist=None, writebricks=False): """Generate skies in parallel for all bricks in a Legacy Surveys Data Release. Parameters @@ -680,17 +738,6 @@ def select_skies(survey, numproc=16, nskiespersqdeg=None, bands=['g', 'r', 'z'], from the input `survey` object and is in the form: `%(survey.survey_dir)/metrics/%(brick).3s/skies-%(brick)s.fits.gz` which is returned by `survey.find_file('skies')`. - bundlebricks : :class:`int`, defaults to None - If not None, then instead of selecting the skies, print, to screen, the slurm - script that will approximately balance the brick distribution at `bundlebricks` - bricks per node. So, for instance, if bundlebricks is 14000 (which as of - the latest git push works well to fit on the interactive nodes on Cori), then - commands would be returned with the correct pixlist values to pass to the code - to pack at about 14000 bricks per node across all of the bricks in `survey`. - brickspersec : :class:`float`, optional, defaults to 1.6 - The rough number of bricks processed per second by the code (parallelized across - a chosen number of nodes). Used in conjunction with `bundlebricks` for the code - to estimate time to completion when parallelizing across pixels. Returns ------- @@ -701,44 +748,25 @@ def select_skies(survey, numproc=16, nskiespersqdeg=None, bands=['g', 'r', 'z'], Notes ----- - Some core code in this module was initially written by Dustin Lang (@dstndstn). - - Returns nothing if bundlebricks is passed (and is not ``None``). """ - # ADM these comments were for debugging photutils/astropy dependencies - # ADM and they can be removed at any time -# import astropy -# print(astropy.version) -# print(astropy.version.version) -# print(photutils.version) -# print(photutils.version.version) - - # ADM read in the survey bricks file, which lists the bricks of interest for this DR - from glob import glob - sbfile = glob(survey.survey_dir+'/*bricks-dr*')[0] - brickinfo = fitsio.read(sbfile) - # ADM remember that fitsio reads things in as bytes, so convert to unicode - bricknames = brickinfo['brickname'].astype('U') - - # ADM if the pixlist or bundlebricks option was sent, we'll need the HEALPpixel - # ADM information for each brick - if pixlist is not None or bundlebricks is not None: - theta, phi = np.radians(90-brickinfo["dec"]), np.radians(brickinfo["ra"]) - pixnum = hp.ang2pix(nside, theta, phi, nest=True) - - # ADM if the bundlebricks option was sent, call the packing code - if bundlebricks is not None: - bundle_bricks(pixnum, bundlebricks, nside, prefix='skies', - surveydir=survey.survey_dir, brickspersec=brickspersec) - return + # ADM retrieve the bricks of interest for this DR. + brickdict = get_brick_info([survey.survey_dir]) + bricknames = np.array(list(brickdict.keys())) - # ADM restrict to only bricks in a set of HEALPixels, if requested + # ADM restrict to only bricks in a set of HEALPixels, if requested. if pixlist is not None: - # ADM if an integer was passed, turn it into a list + bra, bdec, _, _, _, _ = np.vstack(brickdict.values()).T + theta, phi = np.radians(90-bdec), np.radians(bra) + pixnum = hp.ang2pix(nside, theta, phi, nest=True) + # ADM if an integer was passed, turn it into a list. if isinstance(pixlist, int): pixlist = [pixlist] - wbricks = np.where([pix in pixlist for pix in pixnum])[0] - bricknames = bricknames[wbricks] - if len(wbricks) == 0: + ii = [pix in pixlist for pix in pixnum] + bricknames = bricknames[ii] + # ADM if there are no bricks to process, then die immediately. + if len(bricknames) == 0: log.warning('ZERO bricks in passed pixel list!!!') + return log.info("Processing bricks in (nside={}, pixel numbers={}) HEALPixels" .format(nside, pixlist)) From e9dc177a322292302040571b6b37cec950dc0f04 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 4 Apr 2019 10:07:13 -0700 Subject: [PATCH 07/31] update write to facilitate multiple directories for skies; shift bundling test to new geomask module --- bin/select_gfas | 7 ++----- bin/select_skies | 8 ++++---- bin/select_sv_targets | 5 +---- bin/select_targets | 7 ++----- py/desitarget/io.py | 24 ++++++++++++++---------- py/desitarget/randoms.py | 2 +- py/desitarget/skyfibers.py | 2 +- py/desitarget/test/test_skyfibers.py | 7 ------- 8 files changed, 25 insertions(+), 37 deletions(-) diff --git a/bin/select_gfas b/bin/select_gfas index f22ae1aea..ab15ea8a5 100755 --- a/bin/select_gfas +++ b/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.)) diff --git a/bin/select_skies b/bin/select_skies index c75836747..0612f0714 100755 --- a/bin/select_skies +++ b/bin/select_skies @@ -90,7 +90,7 @@ if ns.bundlebricks is not None: 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 + 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) @@ -99,7 +99,7 @@ else: skies = [] for survey in surveys: skies.append(select_skies( - survey, numproc=ns.numproc, nskiespersqdeg=nskiespersqdeg, + survey, numproc=ns.numproc, nskiespersqdeg=nskiespersqdeg, bands=bands, apertures_arcsec=apertures, nside=ns.nside, pixlist=pixlist, writebricks=ns.writebricks) ) @@ -110,7 +110,7 @@ else: # ADM this correctly records the apertures in the output file header # ADM as well as adding HEALPixel information. - io.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)) diff --git a/bin/select_sv_targets b/bin/select_sv_targets index 703145a52..ced130d8f 100755 --- a/bin/select_sv_targets +++ b/bin/select_sv_targets @@ -67,17 +67,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 - indir2 = 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 - indir2 = ns.sweepdir2 if len(infiles) == 0: log.critical('no sweep or tractor files found') sys.exit(1) @@ -129,7 +126,7 @@ if ns.mask: targets = mask_targets(targets, inmaskfile=ns.mask, nside=nside) if ns.bundlefiles is None: - io.write_targets(ns.dest, targets, indir=indir, indir2=indir2, + 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)) diff --git a/bin/select_targets b/bin/select_targets index 1a7f64106..118534927 100755 --- a/bin/select_targets +++ b/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 - indir2 = 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 - indir2 = 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, - indir2=indir2, 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)) diff --git a/py/desitarget/io.py b/py/desitarget/io.py index 6a7365650..ee439a3f3 100644 --- a/py/desitarget/io.py +++ b/py/desitarget/io.py @@ -569,8 +569,8 @@ def write_targets(filename, data, indir=None, indir2=None, fitsio.write(filename, data, extname='TARGETS', header=hdr, clobber=True) -def write_skies(filename, data, indir=None, apertures_arcsec=None, - nskiespersqdeg=None, nside=None): +def write_skies(filename, data, indir=None, indir2=None, + apertures_arcsec=None, nskiespersqdeg=None, nside=None): """Write a target catalogue of sky locations. Parameters @@ -579,9 +579,9 @@ def write_skies(filename, data, indir=None, apertures_arcsec=None, Output target selection file name data : :class:`~numpy.ndarray` Array of skies to write to file. - indir : :class:`str`, optional - Name of input Legacy Survey Data Release directory, write to header - of output file if passed (and if not None). + indir, indir2 : :class:`str`, optional + Name of input Legacy Survey Data Release directory or directories, + write to header of output file if passed (and if not None). apertures_arcsec : :class:`list` or `float`, optional list of aperture radii in arcseconds to write each aperture as an individual line in the header, if passed (and if not None). @@ -610,6 +610,8 @@ def write_skies(filename, data, indir=None, apertures_arcsec=None, # ADM be rewritten gracefully in the header. drstring = 'dr'+indir.split('dr')[-1][0] depend.setdep(hdr, 'photcat', drstring) + if indir2 is not None: + depend.setdep(hdr, 'input-data-release-2', indir2) if apertures_arcsec is not None: for i, ap in enumerate(apertures_arcsec): @@ -636,8 +638,8 @@ def write_skies(filename, data, indir=None, apertures_arcsec=None, fitsio.write(filename, data, extname='SKY_TARGETS', header=hdr, clobber=True) -def write_gfas(filename, data, indir=None, nside=None, survey="?", - gaiaepoch=None): +def write_gfas(filename, data, indir=None, indir2=None, nside=None, + survey="?", gaiaepoch=None): """Write a catalogue of Guide/Focus/Alignment targets. Parameters @@ -646,9 +648,9 @@ def write_gfas(filename, data, indir=None, nside=None, survey="?", Output file name. data : :class:`~numpy.ndarray` Array of GFAs to write to file. - indir : :class:`str`, optional, defaults to None. - Name of input Legacy Survey Data Release directory, write to header - of output file if passed (and if not None). + indir, indir2 : :class:`str`, optional, defaults to None. + Name of input Legacy Survey Data Release directory or directories, + write to header of output file if passed (and if not None). nside: :class:`int`, defaults to None. If passed, add a column to the GFAs array popluated with HEALPixels at resolution `nside`. @@ -673,6 +675,8 @@ def write_gfas(filename, data, indir=None, nside=None, survey="?", # ADM be rewritten gracefully in the header. drstring = 'dr'+indir.split('dr')[-1][0] depend.setdep(hdr, 'photcat', drstring) + if indir2 is not None: + depend.setdep(hdr, 'input-data-release-2', indir2) # ADM add HEALPix column, if requested by input. if nside is not None: diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index 98550d6c7..aebae6114 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -926,7 +926,7 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None, if 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)]) + for cnt, pix in zip(cnts.astype(int), pixnum)]) bundle_bricks(allpixnum, bundlebricks, nside, brickspersec=brickspersec, prefix='randoms', surveydirs=[drdir]) return diff --git a/py/desitarget/skyfibers.py b/py/desitarget/skyfibers.py index e60f70f81..215e4aa5c 100644 --- a/py/desitarget/skyfibers.py +++ b/py/desitarget/skyfibers.py @@ -88,7 +88,7 @@ def get_brick_info(drdirs, counts=False): fns = glob(os.path.join(dd, 'tractor', '*', '*fits')) bricknames.append([brickname_from_filename(fn) for fn in fns]) - # ADM don't count bricks twice. + # ADM don't count bricks twice, but record number of duplicate bricks. bricknames, cnts = np.unique(np.concatenate(bricknames), return_counts=True) # ADM initialize the bricks class, retrieve the brick information look-up diff --git a/py/desitarget/test/test_skyfibers.py b/py/desitarget/test/test_skyfibers.py index 41af58ee6..a7265b676 100644 --- a/py/desitarget/test/test_skyfibers.py +++ b/py/desitarget/test/test_skyfibers.py @@ -162,13 +162,6 @@ def test_target_bits(self): np.all(gskies[wbad]["DESI_TARGET"] == desi_mask.BAD_SKY) ) - def test_bundle_bricks(self): - """ - Test the bundle_bricks scripting code simply executes without bugs - """ - dum = skyfibers.bundle_bricks(1, 1, 1, surveydir=self.survey.survey_dir) - self.assertTrue(dum is None) - def test_select_skies(self): """ Test the wrapper function for batch selection of skies From 91252d885f942b4cd6901d400d136fd3107c46b5 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 4 Apr 2019 11:48:21 -0700 Subject: [PATCH 08/31] add resolve to skies --- bin/select_skies | 6 +++++- py/desitarget/skyfibers.py | 8 +++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/bin/select_skies b/bin/select_skies index 0612f0714..0187c03f4 100755 --- a/bin/select_skies +++ b/bin/select_skies @@ -6,6 +6,7 @@ from desitarget.skyutilities.legacypipe.util import LegacySurveyData 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 @@ -108,9 +109,12 @@ else: ii = [sk is not None for sk in skies] skies = np.concatenate(skies[ii]) + # ADM resolve any duplicates between imaging data releases. + skies = resolve(skies) + # ADM this correctly records the apertures in the output file header # ADM as well as adding HEALPixel information. - io.write_skies(ns.dest, skies, indir=ns.surveydir, indir2=ns.surveydir2, + 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)) diff --git a/py/desitarget/skyfibers.py b/py/desitarget/skyfibers.py index 215e4aa5c..ea9ea2abd 100644 --- a/py/desitarget/skyfibers.py +++ b/py/desitarget/skyfibers.py @@ -298,9 +298,11 @@ def make_skies_for_a_brick(survey, brickname, nskiespersqdeg=None, bands=['g', ' skies["BRICKNAME"] = skytable.brickname skies["BLOBDIST"] = skytable.blobdist - # ADM set the data release from the Legacy Surveys DR directory. - dr = int(survey.survey_dir.split('dr')[-1][0])*1000 - skies["RELEASE"] = dr + # ADM set the data release from an object in a Tractor file. + tfn = survey.find_file("tractor", brick=brickname) + # ADM this file should be guaranteed to exist, except for unit tests. + if os.path.exists(tfn): + skies["RELEASE"] = fitsio.read(tfn, rows=0, columns='RELEASE')[0] # ADM set the objid (just use a sequential number as setting skies # ADM to 1 in the TARGETID will make these unique. From 155529c96d1c77c38887733a3ee955903901f348 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 4 Apr 2019 11:50:06 -0700 Subject: [PATCH 09/31] deprecate make_hpx_density_file as it's been superseded by make_imaging_weight_map --- bin/make_hpx_density_file | 23 ----------------------- 1 file changed, 23 deletions(-) delete mode 100755 bin/make_hpx_density_file diff --git a/bin/make_hpx_density_file b/bin/make_hpx_density_file deleted file mode 100755 index c1bb892a9..000000000 --- a/bin/make_hpx_density_file +++ /dev/null @@ -1,23 +0,0 @@ -#!/usr/bin/env python - -from __future__ import print_function, division - -import os, sys -import numpy as np -import fitsio - -from desitarget.QA import HPX_info - -#import warnings -#warnings.simplefilter('error') - -from argparse import ArgumentParser -ap = ArgumentParser() -ap.add_argument("src", help="Input target file (e.g. /project/projectdirs/desi/target/catalogs/targets-dr3.1-0.14.0.fits)") -ap.add_argument("dest", help="Output file (e.g. /project/projectdirs/desi/target/catalogs/info-dr3.1-0.14.0.fits)") - -ns = ap.parse_args() - -HPX_info(ns.src,outfilename=ns.dest) -print('Density info file written to {}'.format(ns.dest)) - From 673461fa09315356856ae46881f1c388f770a9e1 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 4 Apr 2019 13:28:37 -0700 Subject: [PATCH 10/31] update weight map code to accept HEALPixel-split target files; deprecate some old plotting code --- bin/make_imaging_weight_map | 2 +- py/desitarget/randoms.py | 25 +++++++++---------------- 2 files changed, 10 insertions(+), 17 deletions(-) diff --git a/bin/make_imaging_weight_map b/bin/make_imaging_weight_map index c592b2e84..0a1c2fcc7 100755 --- a/bin/make_imaging_weight_map +++ b/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() diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index aebae6114..9b81b71f0 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -20,6 +20,7 @@ from desitarget.targetmask import desi_mask, bgs_mask, mws_mask from desitarget.targets import resolve from desitarget.skyfibers import get_brick_info +from desitarget.io import read_targets_in_box # ADM the parallelization script from desitarget.internal import sharedmem @@ -532,7 +533,7 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir, return qinfo -def pixweight(randoms, density, nobsgrz=[0, 0, 0], nside=256, outplot=None, outarea=True): +def pixweight(randoms, density, nobsgrz=[0, 0, 0], nside=256, outarea=True): """Fraction of area covered in HEALPixels by a random catalog Parameters @@ -552,9 +553,6 @@ def pixweight(randoms, density, nobsgrz=[0, 0, 0], nside=256, outplot=None, outa than -1) in r-band and z-band. nside : :class:`int`, optional, defaults to nside=256 (~0.0525 sq. deg. or "brick-sized") The resolution (HEALPixel NESTED nside number) at which to build the map. - outplot : :class:`str`, optional, defaults to not making a plot - Create a plot and write it to a file named `outplot` (this is passed to - the `savefig` routine from `matplotlib.pyplot`. outarea : :class:`boolean`, optional, defaults to True Print the total area of the survey for these values of `nobsgrz` to screen. @@ -573,8 +571,6 @@ def pixweight(randoms, density, nobsgrz=[0, 0, 0], nside=256, outplot=None, outa - `0 < WEIGHT < 1` for pixels that partially cover LS DR area with one or more observations. - The index of the array is the HEALPixel integer. """ - import matplotlib.pyplot as plt - # ADM if a file name was passed for the random catalog, read it in if isinstance(randoms, str): randoms = fitsio.read(randoms) @@ -602,12 +598,6 @@ def pixweight(randoms, density, nobsgrz=[0, 0, 0], nside=256, outplot=None, outa # ADM create a weight map based on the actual counts divided by the expected counts pix_weight = pix_cnt/expected_cnt - # ADM if outplot was passed, make a plot of the weights in Mollweide projection - if outplot is not None: - log.info('Plotting pixel map and writing to {}'.format(outplot)) - hp.mollview(pix_weight, nest=True) - plt.savefig(outplot) - # ADM if requested, print the total area of the survey to screen if outarea: area = np.sum(pix_weight*hp.nside2pixarea(nside, degrees=True)) @@ -673,7 +663,7 @@ def stellar_density(nside=256): def get_targ_dens(targets, nside=256): - """The density of targets in HEALPixels + """The density of targets in HEALPixels. Parameters ---------- @@ -739,8 +729,10 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None): A random catalog as made by, e.g., :func:`select_randoms()` or :func:`quantities_at_positions_in_a_brick()`, or the name of such a file. targets : :class:`~numpy.ndarray` or `str` - A corresponding (same Legacy Surveys Data Release) target catalog as made by, - e.g., :func:`desitarget.cuts.select_targets()`, or the name of such a file. + A corresponding (i.e. same Data Release) target catalog as made by, e.g., + :func:`desitarget.cuts.select_targets()`, or a file name of such targets, + or the name of a directory containing HEALPixel-split targets that can + be read by :func:`desitarget.io.read_targets_in_box()`. rand_density : :class:`int` The number of random points per sq. deg. At which the random catalog was generated (see also :func:`select_randoms()`). @@ -781,7 +773,8 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None): # ADM if a file name was passed for the targets catalog, read it in if isinstance(targets, str): log.info('Reading in target catalog...t = {:.1f}s'.format(time()-start)) - targets = fitsio.read(targets) + cols = ["RA", "DEC", "BGS_TARGET", "MWS_TARGET", "DESI_TARGET"] + targets = read_targets_in_box(targets, columns=cols) # ADM determine the areal coverage at of the randoms at this nside log.info('Determining footprint...t = {:.1f}s'.format(time()-start)) From e7c8b6d1bfb849f75a8c923915454682c280ad22 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Thu, 4 Apr 2019 17:42:50 -0700 Subject: [PATCH 11/31] New io function to read target column names from file header without reading in entire file --- py/desitarget/QA.py | 10 ++-------- py/desitarget/io.py | 32 +++++++++++++++++++++++++++++++- py/desitarget/randoms.py | 22 ++++++++++++---------- 3 files changed, 45 insertions(+), 19 deletions(-) diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py index 214b69b01..c2935c980 100644 --- a/py/desitarget/QA.py +++ b/py/desitarget/QA.py @@ -31,7 +31,7 @@ from desitarget.targetmask import desi_mask, bgs_mask, mws_mask from desitarget.cmx.cmx_targetmask import cmx_mask from desitarget.sv1.sv1_targetmask import desi_mask as sv1_desi_mask -from desitarget.io import read_targets_in_box +from desitarget.io import read_targets_in_box, target_columns_from_header # ADM fake the matplotlib display so it doesn't die on allocated nodes. import matplotlib matplotlib.use('Agg') @@ -358,13 +358,7 @@ def read_data(targfile, mocks=False): # ADM from the header of the input files, retrieve the appropriate # ADM names for the SV, main, or cmx _TARGET columns. - targcols = [] - fn = targfile - if os.path.isdir(targfile): - fn = next(iglob(os.path.join(targfile, '*fits'))) - hdr = fitsio.read_header(fn, "TARGETS") - allcols = np.array([hdr[name] if isinstance(hdr[name], str) else 'BLAT' for name in hdr]) - targcols = allcols[['_TARGET' in col for col in allcols]] + targcols = target_columns_from_header(targfile) # ADM limit to the data columns used by the QA code to save memory. colnames = ["RA", "DEC", "RELEASE", "PARALLAX", "PMRA", "PMDEC"] diff --git a/py/desitarget/io.py b/py/desitarget/io.py index ee439a3f3..dbc0abf36 100644 --- a/py/desitarget/io.py +++ b/py/desitarget/io.py @@ -16,7 +16,7 @@ from . import __version__ as desitarget_version import numpy.lib.recfunctions as rfn import healpy as hp -from glob import glob +from glob import glob, iglob from desiutil import depend from desitarget.geomask import hp_in_box, box_area, is_in_box @@ -1403,3 +1403,33 @@ def read_targets_in_cap(hpdirname, radecrad, columns=None): targets = rfn.drop_fields(targets[ii], addedcols) return targets + + +def target_columns_from_header(hpdirname): + """Grab the _TARGET column names from a target file or directory. + + Parameters + ---------- + hpdirname : :class:`str` + Full path to either a directory containing targets that + have been partitioned by HEALPixel (i.e. as made by + `select_targets` with the `bundle_files` option). Or the + name of a single file of targets. + + Returns + ------- + :class:`list` + The names of the _TARGET columns, notably whether they are + SV, main, or cmx _TARGET columns. + """ + # ADM determine whether we're dealing with a file or directory. + fn = hpdirname + if os.path.isdir(hpdirname): + fn = next(iglob(os.path.join(hpdirname, '*fits'))) + + # ADM read in the header and find any columns matching _TARGET. + hdr = fitsio.read_header(fn, "TARGETS") + allcols = np.array([hdr[name] if isinstance(hdr[name], str) else 'BLAT' for name in hdr]) + targcols = allcols[['_TARGET' in col for col in allcols]] + + return targcols diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index 9b81b71f0..5604d3530 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -20,7 +20,7 @@ from desitarget.targetmask import desi_mask, bgs_mask, mws_mask from desitarget.targets import resolve from desitarget.skyfibers import get_brick_info -from desitarget.io import read_targets_in_box +from desitarget.io import read_targets_in_box, target_columns_from_header # ADM the parallelization script from desitarget.internal import sharedmem @@ -773,28 +773,30 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None): # ADM if a file name was passed for the targets catalog, read it in if isinstance(targets, str): log.info('Reading in target catalog...t = {:.1f}s'.format(time()-start)) - cols = ["RA", "DEC", "BGS_TARGET", "MWS_TARGET", "DESI_TARGET"] + # ADM grab appropriate columns for an SV/cmx/main survey file. + targcols = target_columns_from_header(targets) + cols = np.concatenate([["RA", "DEC"], targcols]) targets = read_targets_in_box(targets, columns=cols) - # ADM determine the areal coverage at of the randoms at this nside + # ADM determine the areal coverage of the randoms at this nside. log.info('Determining footprint...t = {:.1f}s'.format(time()-start)) pw = pixweight(randoms, rand_density, nside=nside) npix = len(pw) - # ADM get the target densities + # ADM get the target densities. log.info('Calculating target densities...t = {:.1f}s'.format(time()-start)) targdens = get_targ_dens(targets, nside=nside) - # ADM set up the output array + # ADM set up the output array. datamodel = [('HPXPIXEL', '>i4'), ('FRACAREA', '>f4'), ('STARDENS', '>f4'), ('EBV', '>f4'), ('PSFDEPTH_G', '>f4'), ('PSFDEPTH_R', '>f4'), ('PSFDEPTH_Z', '>f4'), ('GALDEPTH_G', '>f4'), ('GALDEPTH_R', '>f4'), ('GALDEPTH_Z', '>f4')] datamodel += targdens.dtype.descr hpxinfo = np.zeros(npix, dtype=datamodel) - # ADM set initial values to -1 so that they can easily be clipped + # ADM set initial values to -1 so that they can easily be clipped. hpxinfo[...] = -1 - # ADM add the areal coverage, pixel information and target densities + # ADM add the areal coverage, pixel information and target densities. hpxinfo['HPXPIXEL'] = np.arange(npix) hpxinfo['FRACAREA'] = pw for col in targdens.dtype.names: @@ -812,20 +814,20 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None): .format(gaialoc, nside)) hpxinfo["STARDENS"] = sd - # ADM add the median values of all of the other systematics + # ADM add the median values of all of the other systematics. log.info('Calculating medians of systematics from random catalog...t = {:.1f}s' .format(time()-start)) ras, decs = randoms["RA"], randoms["DEC"] pixnums = hp.ang2pix(nside, np.radians(90.-decs), np.radians(ras), nest=True) - # ADM some sorting to order the values to extract the medians + # ADM some sorting to order the values to extract the medians. pixorder = np.argsort(pixnums) pixels, pixcnts = np.unique(pixnums, return_counts=True) pixcnts = np.insert(pixcnts, 0, 0) pixcnts = np.cumsum(pixcnts) # ADM work through the ordered pixels to populate the median for - # ADM each quantity of interest + # ADM each quantity of interest. cols = ['EBV', 'PSFDEPTH_G', 'GALDEPTH_G', 'PSFDEPTH_R', 'GALDEPTH_R', 'PSFDEPTH_Z', 'GALDEPTH_Z'] for i in range(len(pixcnts)-1): From ec5765acde948cf00155c6544274a38230b5fa25 Mon Sep 17 00:00:00 2001 From: Adam Myers Date: Fri, 5 Apr 2019 11:41:25 -0700 Subject: [PATCH 12/31] begin to refactor QA to work with main_cmx_or_sv framework. General QA clean-up --- py/desitarget/QA.py | 131 ++++++++++++++------------------------- py/desitarget/geomask.py | 4 +- py/desitarget/io.py | 2 +- py/desitarget/targets.py | 18 +++++- 4 files changed, 67 insertions(+), 88 deletions(-) diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py index c2935c980..5ee93e356 100644 --- a/py/desitarget/QA.py +++ b/py/desitarget/QA.py @@ -28,10 +28,10 @@ from desiutil import brick from desiutil.log import get_logger from desiutil.plots import init_sky, plot_sky_binned, plot_healpix_map, prepare_data -from desitarget.targetmask import desi_mask, bgs_mask, mws_mask -from desitarget.cmx.cmx_targetmask import cmx_mask -from desitarget.sv1.sv1_targetmask import desi_mask as sv1_desi_mask +from desitarget.targetmask import desi_mask, mws_mask, bgs_mask +from desitarget.targets import main_cmx_or_sv from desitarget.io import read_targets_in_box, target_columns_from_header +from desitarget.geomask import pixarea2nside # ADM fake the matplotlib display so it doesn't die on allocated nodes. import matplotlib matplotlib.use('Agg') @@ -164,8 +164,6 @@ def _load_targdens(tcnames=None, bit_mask=None): desi mask object, e.g., loaded as `from desitarget.targetmask import desi_mask`. Any bit names that contain "NORTH" or "SOUTH" or calibration bits will be removed. - mocks : :class:`boolean`, optional, default=False - If ``True``, also read the expected redshift distributions for each target class. Returns ------- @@ -416,8 +414,7 @@ def qaskymap(cat, objtype, qadir='.', upclip=None, weights=None, max_bin_area=1. A weight for each of the passed targets (e.g., to upweight each target in a partial pixel at the edge of the DESI footprint). max_bin_area : :class:`float`, optional, defaults to 1 degree - The bin size in the passed coordinates is chosen automatically to be as close as - possible to this value without exceeding it. + The bin size in RA/Dec in `targs` is chosen to be as close as possible to this value. fileprefix : :class:`str`, optional, defaults to ``"radec"`` for (RA/Dec) String to be added to the front of the output file name. @@ -561,15 +558,15 @@ def qasystematics_scatterplot(pixmap, syscolname, targcolname, qadir='.', downclip = -1e30 if upclip is None: upclip = 1e30 - w = np.where((pixmap['FRACAREA'] > 0.9) & - (pixmap[syscolname] >= downclip) & (pixmap[syscolname] < upclip))[0] - if len(w) > 0: - pixmapgood = pixmap[w] + ii = ((pixmap['FRACAREA'] > 0.9) & + (pixmap[syscolname] >= downclip) & (pixmap[syscolname] < upclip)) + if np.any(ii): + pixmapgood = pixmap[ii] else: log.error("Pixel map has no areas with >90% coverage for passed up/downclip") log.info("Proceeding without clipping systematics for {}".format(syscolname)) - w = np.where(pixmap['FRACAREA'] > 0.9) - pixmapgood = pixmap[w] + ii = pixmap['FRACAREA'] > 0.9 + pixmapgood = pixmap[ii] # ADM set up the x-axis as the systematic of interest. xx = pixmapgood[syscolname] @@ -631,8 +628,7 @@ def qahisto(cat, objtype, qadir='.', targdens=None, upclip=None, weights=None, m A weight for each of the passed targets (e.g., to upweight each target in a partial pixel at the edge of the DESI footprint). max_bin_area : :class:`float`, optional, defaults to 1 degree - The bin size in the passed coordinates is chosen automatically to be as close as - possible to this value without exceeding it. + The bin size in RA/Dec in `targs` is chosen to be as close as possible to this value. fileprefix : :class:`str`, optional, defaults to ``"histo"`` String to be added to the front of the output file name. catispix : :class:`boolean`, optional, defaults to ``False`` @@ -652,11 +648,7 @@ def qahisto(cat, objtype, qadir='.', targdens=None, upclip=None, weights=None, m import healpy as hp # ADM determine the nside for the passed max_bin_area. - for n in range(1, 25): - nside = 2 ** n - bin_area = hp.nside2pixarea(nside, degrees=True) - if bin_area <= max_bin_area: - break + nside = pixarea2nside(max_bin_area) # ADM the number of HEALPixels and their area at this nside. npix = hp.nside2npix(nside) @@ -1527,8 +1519,7 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True A dictionary of DESI target classes and the goal density for that class. Used to label the goal density on histogram plots. max_bin_area : :class:`float`, optional, defaults to 1 degree - The bin size in the passed coordinates is chosen automatically to be as close as - possible to this value without exceeding it. + The bin size in RA/Dec in `targs` is chosen to be as close as possible to this value. weight : :class:`boolean`, optional, defaults to True If this is set, weight pixels using the ``DESIMODEL`` HEALPix footprint file to ameliorate under dense pixels at the footprint edges. @@ -1585,16 +1576,12 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True else: # ADM if a filename was passed, read in the targets from that file. if isinstance(targs, str): - targs = fitsio.read(targs) + targs = read_targets_in_box(targs) log.info('Read in targets...t = {:.1f}s'.format(time()-start)) truths, objtruths = None, None # ADM determine the nside for the passed max_bin_area. - for n in range(1, 25): - nside = 2 ** n - bin_area = hp.nside2pixarea(nside, degrees=True) - if bin_area <= max_bin_area: - break + nside = pixarea2nside(max_bin_area) # ADM calculate HEALPixel numbers once, here, to avoid repeat calculations # ADM downstream. @@ -1619,10 +1606,10 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True fracarea = pixweight[pix] # ADM weight by 1/(the fraction of each pixel that is in the DESI footprint) # ADM except for zero pixels, which are all outside of the footprint. - w = np.where(fracarea == 0) - fracarea[w] = 1 # ADM to guard against division by zero warnings. + ii = fracarea == 0 + fracarea[ii] = 1 # ADM to guard against division by zero warnings. weights = 1./fracarea - weights[w] = 0 + weights[ii] = 0 # ADM if we have weights, then redetermine the total pix weight. totalpixweight = np.sum(pixweight[uniqpixset]) @@ -1637,7 +1624,7 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True # ADM Current goal target densities for DESI. if targdens is None: - targdens = _load_targdens(tcnames=tcnames, bit_mask=bit_mask, mocks=mocks) + targdens = _load_targdens(tcnames=tcnames, bit_mask=bit_mask) if mocks: dndz = _load_dndz() @@ -1659,59 +1646,59 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True for objtype in targdens: if 'ALL' in objtype: - w = np.arange(len(targs)) + ii = np.ones(len(targs)).astype('bool') else: if ('BGS' in objtype) and not('ANY' in objtype) and not(cmx): - w = np.where(targs["BGS_TARGET"] & bgs_mask[objtype])[0] + ii = targs["BGS_TARGET"] & bgs_mask[objtype] != 0 elif ('MWS' in objtype) and not('ANY' in objtype) and not(cmx): - w = np.where(targs["MWS_TARGET"] & mws_mask[objtype])[0] + ii = targs["MWS_TARGET"] & mws_mask[objtype] != 0 else: - w = np.where(targs["DESI_TARGET"] & main_mask[objtype])[0] + ii = targs["DESI_TARGET"] & main_mask[objtype] != 0 - if len(w) > 0: + if np.any(ii): # ADM make RA/Dec skymaps. - qaskymap(targs[w], objtype, qadir=qadir, upclip=upclipdict[objtype], - weights=weights[w], max_bin_area=max_bin_area) + qaskymap(targs[ii], objtype, qadir=qadir, upclip=upclipdict[objtype], + weights=weights[ii], max_bin_area=max_bin_area) log.info('Made sky map for {}...t = {:.1f}s' .format(objtype, time()-start)) # ADM make histograms of densities. We already calculated the correctly # ADM ordered HEALPixels and so don't need to repeat that calculation. - qahisto(pix[w], objtype, qadir=qadir, targdens=targdens, upclip=upclipdict[objtype], - weights=weights[w], max_bin_area=max_bin_area, catispix=True) + qahisto(pix[ii], objtype, qadir=qadir, targdens=targdens, upclip=upclipdict[objtype], + weights=weights[ii], max_bin_area=max_bin_area, catispix=True) log.info('Made histogram for {}...t = {:.1f}s' .format(objtype, time()-start)) # ADM make color-color plots - qacolor(targs[w], objtype, targs[w], qadir=qadir, fileprefix="color") + qacolor(targs[ii], objtype, targs[ii], qadir=qadir, fileprefix="color") log.info('Made color-color plot for {}...t = {:.1f}s' .format(objtype, time()-start)) # ADM make magnitude histograms - qamag(targs[w], objtype, qadir=qadir, fileprefix="nmag", area=totarea) + qamag(targs[ii], objtype, qadir=qadir, fileprefix="nmag", area=totarea) log.info('Made magnitude histogram plot for {}...t = {:.1f}s' .format(objtype, time()-start)) if truths is not None: # ADM make noiseless color-color plots - qacolor(truths[w], objtype, targs[w], qadir=qadir, mocks=True, + qacolor(truths[ii], objtype, targs[ii], qadir=qadir, mocks=True, fileprefix="mock-color", nodustcorr=True) log.info('Made (mock) color-color plot for {}...t = {:.1f}s' .format(objtype, time()-start)) # ADM make N(z) plots - mock_qanz(truths[w], objtype, qadir=qadir, area=totarea, dndz=dndz, + mock_qanz(truths[ii], objtype, qadir=qadir, area=totarea, dndz=dndz, fileprefixz="mock-nz", fileprefixzmag="mock-zvmag") log.info('Made (mock) redshift plots for {}...t = {:.1f}s' .format(objtype, time()-start)) # # ADM plot what fraction of each selected object is actually a contaminant - # mock_qafractype(truths[w], objtype, qadir=qadir, fileprefix="mock-fractype") + # mock_qafractype(truths[ii], objtype, qadir=qadir, fileprefix="mock-fractype") # log.info('Made (mock) classification fraction plots for {}...t = {:.1f}s'.format(objtype, time()-start)) # ADM make Gaia-based plots if we have Gaia columns - if "PARALLAX" in targs.dtype.names and np.sum(targs['PARALLAX'] != 0) > 0: - qagaia(targs[w], objtype, qadir=qadir, fileprefix="gaia") + if "PARALLAX" in targs.dtype.names and np.any(targs['PARALLAX'] != 0): + qagaia(targs[ii], objtype, qadir=qadir, fileprefix="gaia") log.info('Made Gaia-based plots for {}...t = {:.1f}s' .format(objtype, time()-start)) @@ -1736,8 +1723,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' makeplots : :class:`boolean`, optional, default=True If ``True``, then create the plots as well as the webpage. max_bin_area : :class:`float`, optional, defaults to 1 degree - The bin size in the passed coordinates is chosen automatically to be as close as - possible to this value without exceeding it. + The bin size in RA/Dec in `targs` is chosen to be as close as possible to this value. qadir : :class:`str`, optional, defaults to the current directory The output directory to which to write produced plots. clip2foot : :class:`boolean`, optional, defaults to False @@ -1792,33 +1778,13 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' truths, objtruths = None, None # ADM automatically detect whether we're running this for the main survey - # ADM or SV, etc. and change the column names accordingly. - svs = "DESI" # ADM this is to store sv iteration or cmx as a string. - colnames = np.array(targs.dtype.names) - svcolnames = colnames[['SV' in name or 'CMX' in name for name in colnames]] - # ADM set cmx flag to True if 'CMX_TARGET' is a column and rename that column. - cmx = 'CMX_TARGET' in svcolnames - targs = rfn.rename_fields(targs, {'CMX_TARGET': 'DESI_TARGET'}) - # ADM strip "SVX" off any columns (rfn.rename_fields forgives missing fields). - for field in svcolnames: - svs = field.split('_')[0] - targs = rfn.rename_fields(targs, {field: "_".join(field.split('_')[1:])}) - - # ADM use the commissioning mask bits/names if we have a CMX file. - # ADM note desi_mask was changed to the SV mask if SV columns exist. - if cmx: - main_mask = cmx_mask - elif svs == 'SV1': - main_mask = sv1_desi_mask - else: - main_mask = desi_mask + # ADM or SV, etc. and load the appropriate mask. + colnames, masks, svs, targs = main_cmx_or_sv(targs, rename=True) + svs = svs.upper() + cmx = svs == 'CMX' # ADM determine the working nside for the passed max_bin_area. - for n in range(1, 25): - nside = 2 ** n - bin_area = hp.nside2pixarea(nside, degrees=True) - if bin_area <= max_bin_area: - break + nside = pixarea2nside(max_bin_area) # ADM if requested, restrict the targets (and mock files) to the DESI footprint. if clip2foot: @@ -1840,11 +1806,10 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' # ADM Set up the names of the target classes and their goal densities using # ADM the goal target densities for DESI (read from the DESIMODEL defaults). - if svs == "DESI": + if svs == "MAIN": targdens = _load_targdens(tcnames=tcnames) - # ADM if we aren't looking at the main survey, pass the cmx or SV mask instead. else: - targdens = _load_targdens(tcnames=tcnames, bit_mask=main_mask) + targdens = _load_targdens(tcnames=tcnames, bit_mask=masks[0]) # ADM set up the html file and write preamble to it. htmlfile = makepath(os.path.join(qadir, 'index.html')) @@ -1872,7 +1837,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' # ADM html preamble. html.write('\n') - html.write('

DESI Targeting QA pages - {} ({})

\n'.format(objtype, DRs)) + html.write('

{} Targeting QA pages - {} ({})

\n'.format(svs, objtype, DRs)) # ADM Target Densities. html.write('

Target density

\n') @@ -2006,7 +1971,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' # ADM make the QA plots, if requested: if makeplots: - if svs == "DESI": + if svs == "MAIN": totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths, qadir=qadir, targdens=targdens, max_bin_area=max_bin_area, weight=weight, imaging_map_file=imaging_map_file, mocks=mocks) @@ -2014,7 +1979,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths, qadir=qadir, targdens=targdens, max_bin_area=max_bin_area, weight=weight, imaging_map_file=imaging_map_file, - cmx=cmx, bit_mask=main_mask, mocks=mocks) + cmx=cmx, bit_mask=masks[0], mocks=mocks) # ADM add a correlation matrix recording the overlaps between different target # ADM classes as a density. @@ -2022,7 +1987,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.' htmlmain.write('

Overlaps in target densities (per sq. deg.)

\n') htmlmain.write('
\n')
         # ADM only retain classes that are actually in the DESI target bit list.
-        settargdens = set(main_mask.names()).intersection(set(targdens))
+        settargdens = set(masks[0].names()).intersection(set(targdens))
         # ADM write out a list of the target categories.
         headerlist = list(settargdens)
         headerlist.sort()
@@ -2042,7 +2007,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
                     overlaps.append(" ")
                 else:
                     dt = targs["DESI_TARGET"]
-                    overlap = np.sum(((dt & main_mask[objtype1]) != 0) & ((dt & main_mask[objtype2]) != 0))/totarea
+                    overlap = np.sum(((dt & masks[0][objtype1]) != 0) & ((dt & masks[0][objtype2]) != 0))/totarea
                     overlaps.append("{:.1f}".format(overlap))
             htmlmain.write(" ".join([truncform.format(i) for i in overlaps])+'\n\n')
         # ADM close the matrix text output.
diff --git a/py/desitarget/geomask.py b/py/desitarget/geomask.py
index a6aa6f8c7..b7e13a06c 100644
--- a/py/desitarget/geomask.py
+++ b/py/desitarget/geomask.py
@@ -851,10 +851,10 @@ def hp_in_box(nside, radecbox, inclusive=True, fact=4):
     leeway = 1e-5
     if decmax > 90-leeway:
         decmax = 90-leeway
-        log.warning('Max Dec too close to pole; set to {}o'.format(decmax))
+        log.warning('Max Dec too close to pole; setting to {}o'.format(decmax))
     if decmin < -90+leeway:
         decmin = -90+leeway
-        log.warning('Min Dec too close to pole; set to {}o'.format(decmin))
+        log.warning('Min Dec too close to pole; setting to {}o'.format(decmin))
 
     # ADM area enclosed isn't well-defined if RA covers more than 180o.
     if np.abs(ramax-ramin) > 180.:
diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index dbc0abf36..50cc7512f 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -1426,7 +1426,7 @@ def target_columns_from_header(hpdirname):
     fn = hpdirname
     if os.path.isdir(hpdirname):
         fn = next(iglob(os.path.join(hpdirname, '*fits')))
-    
+
     # ADM read in the header and find any columns matching _TARGET.
     hdr = fitsio.read_header(fn, "TARGETS")
     allcols = np.array([hdr[name] if isinstance(hdr[name], str) else 'BLAT' for name in hdr])
diff --git a/py/desitarget/targets.py b/py/desitarget/targets.py
index 89cf40cb8..d92815eff 100644
--- a/py/desitarget/targets.py
+++ b/py/desitarget/targets.py
@@ -283,7 +283,7 @@ def decode_targetid(targetid):
     return objid, brickid, release, mock, sky
 
 
-def main_cmx_or_sv(targets):
+def main_cmx_or_sv(targets, rename=False):
     """determine whether a target array is main survey, commissioning, or SV
 
     Parameters
@@ -292,6 +292,9 @@ def main_cmx_or_sv(targets):
         An array of targets generated by, e.g., :mod:`~desitarget.cuts`
         must include at least (all of) the columns `DESI_TARGET`, `MWS_TARGET` and
         `BGS_TARGET` or the corresponding commissioning or SV columns.
+    rename : :class:`bool`, optional, defaults to ``False``
+        If ``True`` then also return a copy of `targets` with the input `_TARGET`
+        columns renamed to reflect the main survey format.
 
     Returns
     -------
@@ -306,9 +309,13 @@ def main_cmx_or_sv(targets):
     :class:`str`
         The string 'main', 'cmx' or 'svX' (where X = 1, 2, 3 etc.) for the main survey,
         commissioning and an iteration of SV. Specifies which type of file was sent.
+    :class:`~numpy.ndarray`, optional, if `rename` is ``True``
+        A copy of the input targets array with the `_TARGET` columns renamed to
+        `DESI_TARGET`, and (if they exist) `BGS_TARGET`, `MWS_TARGET`.
     """
     # ADM default to the main survey.
-    outcolnames = ["DESI_TARGET", "BGS_TARGET", "MWS_TARGET"]
+    maincolnames = ["DESI_TARGET", "BGS_TARGET", "MWS_TARGET"]
+    outcolnames = maincolnames.copy()
     masks = [desi_mask, bgs_mask, mws_mask]
     survey = 'main'
 
@@ -335,6 +342,13 @@ def main_cmx_or_sv(targets):
         log.critical(msg)
         raise ValueError(msg)
 
+    # ADM if requested, rename the columns.
+    if rename:
+        mapper = {}
+        for i, col in enumerate(outcolnames):
+            mapper[col] = maincolnames[i]
+        return outcolnames, masks, survey, rfn.rename_fields(targets, mapper)
+
     return outcolnames, masks, survey
 
 

From 4031503f65eb0934237ae74c43763e600a2b0b55 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Fri, 5 Apr 2019 13:24:09 -0700
Subject: [PATCH 13/31] update QA to plot BGS/MWS targets for SV

---
 py/desitarget/QA.py | 41 ++++++++++++++++++++++++-----------------
 1 file changed, 24 insertions(+), 17 deletions(-)

diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 5ee93e356..43581bfab 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -28,7 +28,7 @@
 from desiutil import brick
 from desiutil.log import get_logger
 from desiutil.plots import init_sky, plot_sky_binned, plot_healpix_map, prepare_data
-from desitarget.targetmask import desi_mask, mws_mask, bgs_mask
+from desitarget.targetmask import desi_mask, bgs_mask, mws_mask
 from desitarget.targets import main_cmx_or_sv
 from desitarget.io import read_targets_in_box, target_columns_from_header
 from desitarget.geomask import pixarea2nside
@@ -158,12 +158,12 @@ def _load_targdens(tcnames=None, bit_mask=None):
     tcnames : :class:`list`
         A list of strings, e.g. "['QSO','LRG','ALL'] If passed, return only a dictionary
         for those specific bits.
-    bit_mask : :class:`~numpy.array`, optional, defaults to ``None``
+    bit_mask : :class:`list` or `~numpy.array`, optional, defaults to ``None``
         If passed, load the bit names from this mask (with no associated expected
         densities) rather than loading the main survey bits and densities. Must be a
         desi mask object, e.g., loaded as `from desitarget.targetmask import desi_mask`.
         Any bit names that contain "NORTH" or "SOUTH" or calibration bits will be
-        removed.
+        removed. A list of serveral masks can be passed rather than a single mask.
 
     Returns
     -------
@@ -202,10 +202,14 @@ def _load_targdens(tcnames=None, bit_mask=None):
         targdens['MWS_WD'] = 0.
         targdens['MWS_NEARBY'] = 0.
     else:
-        # ADM this is the list of words contained in bits that we don't want to consider for QA.
-        badnames = ["NORTH", "SOUTH", "NO_TARGET", "SECONDARY", "BRIGHT_OBJECT", "SKY"]
-        names = [name for name in bit_mask.names() if not any(badname in name for badname in badnames)]
-        targdens = {k: 0. for k in names}
+        bit_masks = np.atleast_1d(bit_mask)
+        names = []
+        for bit_mask in bit_masks:
+            # ADM this is the list of words contained in bits that we don't want to consider for QA.
+            badnames = ["NORTH", "SOUTH", "NO_TARGET", "SECONDARY", "BRIGHT_OBJECT", "SKY"]
+            names.append([name for name in bit_mask.names()
+                          if not any(badname in name for badname in badnames)])
+        targdens = {k: 0. for k in np.concatenate(names)}
 
     if tcnames is None:
         return targdens
@@ -1538,9 +1542,9 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
     cmx : :class:`boolean`, defaults to ``False``
         Pass as ``True`` to operate on commissioning bits instead of SV or main survey
         bits. Commissioning files have no MWS or BGS columns.
-    bit_mask : :class:`~numpy.array`, optional, defaults to ``None``
-        Load the bit names from this passed mask (with zero density constraints)
-        instead of the main survey bits.
+    bit_mask : :class:`list` or `~numpy.array`, optional, defaults to ``None``
+        Load the bit names from this passed mask or list of masks (with zero density
+        constraints) instead of the main survey bits.
     mocks : :class:`boolean`, optional, default=False
         If ``True``, add plots that are only relevant to mocks at the bottom of the webpage.
 
@@ -1633,9 +1637,12 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
     # ADM by rejecting highly dense outliers.
     upclipdict = {k: 5000. for k in targdens}
     if bit_mask is not None:
-        main_mask = bit_mask
+        if cmx:
+            d_mask = bit_mask[0]
+        else:
+            d_mask, b_mask, m_mask = bit_mask
     else:
-        main_mask = desi_mask
+        d_mask, b_mask, m_mask = desi_mask, bgs_mask, mws_mask
         upclipdict = {'ELG': 4000, 'LRG': 1200, 'QSO': 400, 'ALL': 8000,
                       'STD_FAINT': 300, 'STD_BRIGHT': 300,
                       # 'STD_FAINT': 200, 'STD_BRIGHT': 50,
@@ -1649,11 +1656,11 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
             ii = np.ones(len(targs)).astype('bool')
         else:
             if ('BGS' in objtype) and not('ANY' in objtype) and not(cmx):
-                ii = targs["BGS_TARGET"] & bgs_mask[objtype] != 0
+                ii = targs["BGS_TARGET"] & b_mask[objtype] != 0
             elif ('MWS' in objtype) and not('ANY' in objtype) and not(cmx):
-                ii = targs["MWS_TARGET"] & mws_mask[objtype] != 0
+                ii = targs["MWS_TARGET"] & m_mask[objtype] != 0
             else:
-                ii = targs["DESI_TARGET"] & main_mask[objtype] != 0
+                ii = targs["DESI_TARGET"] & d_mask[objtype] != 0
 
         if np.any(ii):
             # ADM make RA/Dec skymaps.
@@ -1809,7 +1816,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
     if svs == "MAIN":
         targdens = _load_targdens(tcnames=tcnames)
     else:
-        targdens = _load_targdens(tcnames=tcnames, bit_mask=masks[0])
+        targdens = _load_targdens(tcnames=tcnames, bit_mask=masks)
 
     # ADM set up the html file and write preamble to it.
     htmlfile = makepath(os.path.join(qadir, 'index.html'))
@@ -1979,7 +1986,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
             totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths,
                                     qadir=qadir, targdens=targdens, max_bin_area=max_bin_area,
                                     weight=weight, imaging_map_file=imaging_map_file,
-                                    cmx=cmx, bit_mask=masks[0], mocks=mocks)
+                                    cmx=cmx, bit_mask=masks, mocks=mocks)
 
         # ADM add a correlation matrix recording the overlaps between different target
         # ADM classes as a density.

From 8aaebff757f3df12b229b0e8dafa59c5b39753aa Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Fri, 5 Apr 2019 15:17:48 -0700
Subject: [PATCH 14/31] pixel weight maps and full QA for SV. Won't fully work
 until we've finalized SV bits; basic header functionality for
 read_targets_in_a_box

---
 bin/make_imaging_weight_map |  6 ++--
 bin/run_target_qa           | 16 ++++++++-
 py/desitarget/QA.py         | 14 ++++----
 py/desitarget/io.py         | 68 ++++++++++++++++++++++++++++++++-----
 py/desitarget/randoms.py    | 21 +++++++++---
 5 files changed, 103 insertions(+), 22 deletions(-)

diff --git a/bin/make_imaging_weight_map b/bin/make_imaging_weight_map
index 0a1c2fcc7..ccc259ab5 100755
--- a/bin/make_imaging_weight_map
+++ b/bin/make_imaging_weight_map
@@ -44,11 +44,13 @@ if not os.path.exists(ns.targets):
     sys.exit(1)
 
 hdr = fitsio.read_header(ns.randoms,"RANDOMS")
-#ADM add HEALPixel information to the header
+#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))
diff --git a/bin/run_target_qa b/bin/run_target_qa
index e5add1f8a..62a536b54 100755
--- a/bin/run_target_qa
+++ b/bin/run_target_qa
@@ -7,6 +7,7 @@ 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()
@@ -32,11 +33,24 @@ ap.add_argument('--nosystematics', action='store_true',
 
 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)
diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 43581bfab..1594bf410 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -169,9 +169,15 @@ def _load_targdens(tcnames=None, bit_mask=None):
     -------
     :class:`dictionary`
         A dictionary where the keys are the bit names and the values are the densities.
+    
+    Notes
+    -----
+        If `bit_mask` happens to correpond to the main survey masks, then the default
+        behavior is triggered (as if `bit_mask=None`).
     """
+    bit_masks = np.atleast_1d(bit_mask)
 
-    if bit_mask is None:
+    if bit_mask is None or bit_masks[0]._name == 'desi_mask':
         from desimodel import io
         targdict = io.load_target_info()
 
@@ -202,7 +208,6 @@ def _load_targdens(tcnames=None, bit_mask=None):
         targdens['MWS_WD'] = 0.
         targdens['MWS_NEARBY'] = 0.
     else:
-        bit_masks = np.atleast_1d(bit_mask)
         names = []
         for bit_mask in bit_masks:
             # ADM this is the list of words contained in bits that we don't want to consider for QA.
@@ -1813,10 +1818,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
 
     # ADM Set up the names of the target classes and their goal densities using
     # ADM the goal target densities for DESI (read from the DESIMODEL defaults).
-    if svs == "MAIN":
-        targdens = _load_targdens(tcnames=tcnames)
-    else:
-        targdens = _load_targdens(tcnames=tcnames, bit_mask=masks)
+    targdens = _load_targdens(tcnames=tcnames, bit_mask=masks)
 
     # ADM set up the html file and write preamble to it.
     htmlfile = makepath(os.path.join(qadir, 'index.html'))
diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index 50cc7512f..c7e2ca096 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -1229,7 +1229,8 @@ def check_hp_target_dir(hpdirname):
     return nside[0], pixdict
 
 
-def read_targets_in_hp(hpdirname, nside, pixlist, columns=None):
+def read_targets_in_hp(hpdirname, nside, pixlist, columns=None, 
+                       header=False):
     """Read in targets in a set of HEALPixels.
 
     Parameters
@@ -1245,11 +1246,19 @@ def read_targets_in_hp(hpdirname, nside, pixlist, columns=None):
         Return targets in these HEALPixels at the passed `nside`.
     columns : :class:`list`, optional
         Only read in these target columns.
+    header : :class:`bool`, optional, defaults to ``False``
+        If ``True`` then return the header of either the `hpdirname`
+        file, or the last file read from the `hpdirname` directory.
 
     Returns
     -------
     :class:`~numpy.ndarray`
         An array of targets in the passed pixels.
+
+    Notes
+    -----
+        - If `header` is ``True``, then a second output (the file
+          header is returned).
     """
     # ADM we'll need RA/Dec for final cuts, so ensure they're read.
     addedcols = []
@@ -1280,22 +1289,28 @@ def read_targets_in_hp(hpdirname, nside, pixlist, columns=None):
         # ADM read in the files and concatenate the resulting targets.
         targets = []
         for infile in infiles:
-            targets.append(fitsio.read(infile, columns=columnscopy))
-        targets = np.concatenate(targets)
+            targs, hdr = fitsio.read(infile, 'TARGETS',
+                                     columns=columnscopy, header=True)
+            targets.append(targs)
+        targets = np.concatenate(targets)            
     # ADM ...otherwise just read in the targets.
     else:
-        targets = fitsio.read(hpdirname, columns=columnscopy)
+        targets, hdr = fitsio.read_header(hpdirname, 'TARGETS',
+                                          columns=columnscopy,
+                                          header=True)
 
     # ADM restrict the targets to the actual requested HEALPixels...
     ii = is_in_hp(targets, nside, pixlist)
     # ADM ...and remove RA/Dec columns if we added them.
     targets = rfn.drop_fields(targets[ii], addedcols)
 
+    if header:
+        return targets, hdr
     return targets
 
 
 def read_targets_in_box(hpdirname, radecbox=[0., 360., -90., 90.],
-                        columns=None):
+                        columns=None, header=False):
     """Read in targets in an RA/Dec box.
 
     Parameters
@@ -1310,11 +1325,19 @@ def read_targets_in_box(hpdirname, radecbox=[0., 360., -90., 90.],
         forming the edges of a box in RA/Dec (degrees).
     columns : :class:`list`, optional
         Only read in these target columns.
+    header : :class:`bool`, optional, defaults to ``False``
+        If ``True`` then return the header of either the `hpdirname`
+        file, or the last file read from the `hpdirname` directory.
 
     Returns
     -------
     :class:`~numpy.ndarray`
         An array of targets in the passed RA/Dec box.
+
+    Notes
+    -----
+        - If `header` is ``True``, then a second output (the file
+          header is returned).
     """
     # ADM we'll need RA/Dec for final cuts, so ensure they're read.
     addedcols = []
@@ -1336,17 +1359,22 @@ def read_targets_in_box(hpdirname, radecbox=[0., 360., -90., 90.],
         pixlist = hp_in_box(nside, radecbox)
 
         # ADM read in targets in these HEALPixels.
-        targets = read_targets_in_hp(hpdirname, nside, pixlist,
-                                     columns=columnscopy)
+        targets, hdr = read_targets_in_hp(hpdirname, nside, pixlist,
+                                          columns=columnscopy,
+                                          header=True)
     # ADM ...otherwise just read in the targets.
     else:
-        targets = fitsio.read(hpdirname, columns=columnscopy)
+        targets, hdr = fitsio.read(hpdirname, 'TARGETS',
+                                   columns=columnscopy,
+                                   header=True)
 
     # ADM restrict only to targets in the requested RA/Dec box...
     ii = is_in_box(targets, radecbox)
     # ADM ...and remove RA/Dec columns if we added them.
     targets = rfn.drop_fields(targets[ii], addedcols)
 
+    if header:
+        return targets, hdr
     return targets
 
 
@@ -1405,6 +1433,30 @@ def read_targets_in_cap(hpdirname, radecrad, columns=None):
     return targets
 
 
+def read_targets_in_box_header(hpdirname):
+    """Read in targets in an RA/Dec box.
+
+    Parameters
+    ----------
+    hpdirname : :class:`str`
+        Full path to either a directory containing targets that
+        have been partitioned by HEALPixel (i.e. as made by
+        `select_targets` with the `bundle_files` option). Or the
+        name of a single file of targets.
+
+    Returns
+    -------
+    :class:`FITSHDR`
+        The header of `hpdirname` if it is a file, or the header
+        of the first file encountered in `hpdirname`
+    """
+    if os.path.isdir(hpdirname):
+        gen = iglob(os.path.join(hpdirname, '*fits'))
+        hpdirname = next(gen)
+
+    return fitsio.read_header(hpdirname, 'TARGETS')
+
+
 def target_columns_from_header(hpdirname):
     """Grab the _TARGET column names from a target file or directory.
 
diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py
index 5604d3530..e18b324f0 100644
--- a/py/desitarget/randoms.py
+++ b/py/desitarget/randoms.py
@@ -18,7 +18,7 @@
 from desitarget.gaiamatch import _get_gaia_dir
 from desitarget.geomask import bundle_bricks, box_area
 from desitarget.targetmask import desi_mask, bgs_mask, mws_mask
-from desitarget.targets import resolve
+from desitarget.targets import resolve, main_cmx_or_sv
 from desitarget.skyfibers import get_brick_info
 from desitarget.io import read_targets_in_box, target_columns_from_header
 
@@ -662,7 +662,7 @@ def stellar_density(nside=256):
     return pixout/pixarea
 
 
-def get_targ_dens(targets, nside=256):
+def get_targ_dens(targets, nside=256, bit_mask=None):
     """The density of targets in HEALPixels.
 
     Parameters
@@ -672,6 +672,12 @@ def get_targ_dens(targets, nside=256):
         e.g., :func:`desitarget.cuts.select_targets()`, or the name of such a file.
     nside : :class:`int`, optional, defaults to nside=256 (~0.0525 sq. deg. or "brick-sized")
         The resolution (HEALPixel nside number) at which to build the map (NESTED scheme).
+    bit_mask : :class:`list` or `~numpy.array`, optional, defaults to ``None``
+        If passed, load the bit names from this mask (with no associated expected
+        densities) rather than loading the main survey bits and densities. Must be a
+        desi mask object, e.g., loaded as `from desitarget.targetmask import desi_mask`.
+        Any bit names that contain "NORTH" or "SOUTH" or calibration bits will be
+        removed. A list of serveral masks can be passed rather than a single mask.
 
     Returns
     -------
@@ -695,7 +701,7 @@ def get_targ_dens(targets, nside=256):
 
     # ADM retrieve the bit names of interest
     from desitarget.QA import _load_targdens
-    bitnames = np.array(list(_load_targdens().keys()))
+    bitnames = np.array(list(_load_targdens(bit_mask=bit_mask).keys()))
 
     # ADM and set up an array to hold the output target densities
     targdens = np.zeros(npix, dtype=[(bitname, 'f4') for bitname in bitnames])
@@ -761,6 +767,9 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
                                 the median of GALDEPTH values in the passed random catalog.
             - One column for every bit returned by :func:`desitarget.QA._load_targdens()`.
               Each column contains the density of targets in pixels at the passed `nside`
+    :class:`str`
+        The type of survey to which `targets` corresponds, e.g., 'main', 'sv1', etc.
+
     Notes
     -----
         - If `gaialoc` is ``None`` then the environment variable $GAIA_DIR must be set.
@@ -777,6 +786,8 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
         targcols = target_columns_from_header(targets)
         cols = np.concatenate([["RA", "DEC"], targcols])
         targets = read_targets_in_box(targets, columns=cols)
+    # ADM change target column names, and retrieve associated survey information.
+    _, bit_mask, survey, targets = main_cmx_or_sv(targets, rename=True)
 
     # ADM determine the areal coverage of the randoms at this nside.
     log.info('Determining footprint...t = {:.1f}s'.format(time()-start))
@@ -785,7 +796,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
 
     # ADM get the target densities.
     log.info('Calculating target densities...t = {:.1f}s'.format(time()-start))
-    targdens = get_targ_dens(targets, nside=nside)
+    targdens = get_targ_dens(targets, nside=nside, bit_mask=bit_mask)
 
     # ADM set up the output array.
     datamodel = [('HPXPIXEL', '>i4'), ('FRACAREA', '>f4'), ('STARDENS', '>f4'), ('EBV', '>f4'),
@@ -838,7 +849,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
 
     log.info('Done...t = {:.1f}s'.format(time()-start))
 
-    return hpxinfo
+    return hpxinfo, survey
 
 
 def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None,

From 10da98a85e84de0d00d218062dbb919e9db3f2c0 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Fri, 5 Apr 2019 19:23:01 -0700
Subject: [PATCH 15/31] pixel weights and QA now have full functionality for SV

---
 bin/make_imaging_weight_map |  2 +-
 py/desitarget/QA.py         | 10 ++++-----
 py/desitarget/io.py         |  4 ++--
 py/desitarget/randoms.py    | 42 ++++++++++++++++++++-----------------
 4 files changed, 31 insertions(+), 27 deletions(-)

diff --git a/bin/make_imaging_weight_map b/bin/make_imaging_weight_map
index ccc259ab5..0468a06d6 100755
--- a/bin/make_imaging_weight_map
+++ b/bin/make_imaging_weight_map
@@ -43,7 +43,7 @@ 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")
+hdr = fitsio.read_header(ns.randoms, "RANDOMS")
 #ADM add HEALPixel and gaialoc information to the header
 hdr['GAIALOC'] = ns.gaialoc
 hdr['HPXNSIDE'] = ns.nside
diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 1594bf410..e8a948b80 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -163,13 +163,13 @@ def _load_targdens(tcnames=None, bit_mask=None):
         densities) rather than loading the main survey bits and densities. Must be a
         desi mask object, e.g., loaded as `from desitarget.targetmask import desi_mask`.
         Any bit names that contain "NORTH" or "SOUTH" or calibration bits will be
-        removed. A list of serveral masks can be passed rather than a single mask.
+        removed. A list of several masks can be passed rather than a single mask.
 
     Returns
     -------
     :class:`dictionary`
         A dictionary where the keys are the bit names and the values are the densities.
-    
+
     Notes
     -----
         If `bit_mask` happens to correpond to the main survey masks, then the default
@@ -211,7 +211,7 @@ def _load_targdens(tcnames=None, bit_mask=None):
         names = []
         for bit_mask in bit_masks:
             # ADM this is the list of words contained in bits that we don't want to consider for QA.
-            badnames = ["NORTH", "SOUTH", "NO_TARGET", "SECONDARY", "BRIGHT_OBJECT", "SKY"]
+            badnames = ["NORTH", "SOUTH", "NO_TARGET", "SECONDARY", "BRIGHT_OBJECT", "SKY", "KNOWN"]
             names.append([name for name in bit_mask.names()
                           if not any(badname in name for badname in badnames)])
         targdens = {k: 0. for k in np.concatenate(names)}
@@ -1660,9 +1660,9 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
         if 'ALL' in objtype:
             ii = np.ones(len(targs)).astype('bool')
         else:
-            if ('BGS' in objtype) and not('ANY' in objtype) and not(cmx):
+            if ('BGS' in objtype) and not('S_ANY' in objtype) and not(cmx):
                 ii = targs["BGS_TARGET"] & b_mask[objtype] != 0
-            elif ('MWS' in objtype) and not('ANY' in objtype) and not(cmx):
+            elif ('MWS' in objtype) and not('S_ANY' in objtype) and not(cmx):
                 ii = targs["MWS_TARGET"] & m_mask[objtype] != 0
             else:
                 ii = targs["DESI_TARGET"] & d_mask[objtype] != 0
diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index c7e2ca096..c4d1401a6 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -1229,7 +1229,7 @@ def check_hp_target_dir(hpdirname):
     return nside[0], pixdict
 
 
-def read_targets_in_hp(hpdirname, nside, pixlist, columns=None, 
+def read_targets_in_hp(hpdirname, nside, pixlist, columns=None,
                        header=False):
     """Read in targets in a set of HEALPixels.
 
@@ -1292,7 +1292,7 @@ def read_targets_in_hp(hpdirname, nside, pixlist, columns=None,
             targs, hdr = fitsio.read(infile, 'TARGETS',
                                      columns=columnscopy, header=True)
             targets.append(targs)
-        targets = np.concatenate(targets)            
+        targets = np.concatenate(targets)
     # ADM ...otherwise just read in the targets.
     else:
         targets, hdr = fitsio.read_header(hpdirname, 'TARGETS',
diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py
index e18b324f0..8d89f96a7 100644
--- a/py/desitarget/randoms.py
+++ b/py/desitarget/randoms.py
@@ -17,7 +17,6 @@
 from glob import glob
 from desitarget.gaiamatch import _get_gaia_dir
 from desitarget.geomask import bundle_bricks, box_area
-from desitarget.targetmask import desi_mask, bgs_mask, mws_mask
 from desitarget.targets import resolve, main_cmx_or_sv
 from desitarget.skyfibers import get_brick_info
 from desitarget.io import read_targets_in_box, target_columns_from_header
@@ -662,7 +661,7 @@ def stellar_density(nside=256):
     return pixout/pixarea
 
 
-def get_targ_dens(targets, nside=256, bit_mask=None):
+def get_targ_dens(targets, Mx, nside=256):
     """The density of targets in HEALPixels.
 
     Parameters
@@ -670,14 +669,11 @@ def get_targ_dens(targets, nside=256, bit_mask=None):
     targets : :class:`~numpy.ndarray` or `str`
         A corresponding (same Legacy Surveys Data Release) target catalog as made by,
         e.g., :func:`desitarget.cuts.select_targets()`, or the name of such a file.
+    Mx : :class:`list` or `~numpy.array`
+        The targeting bitmasks associated with the passed targets, assumed to be 
+        a desi, bgs and mws mask in that order (for either SV or the main survey).
     nside : :class:`int`, optional, defaults to nside=256 (~0.0525 sq. deg. or "brick-sized")
         The resolution (HEALPixel nside number) at which to build the map (NESTED scheme).
-    bit_mask : :class:`list` or `~numpy.array`, optional, defaults to ``None``
-        If passed, load the bit names from this mask (with no associated expected
-        densities) rather than loading the main survey bits and densities. Must be a
-        desi mask object, e.g., loaded as `from desitarget.targetmask import desi_mask`.
-        Any bit names that contain "NORTH" or "SOUTH" or calibration bits will be
-        removed. A list of serveral masks can be passed rather than a single mask.
 
     Returns
     -------
@@ -691,6 +687,14 @@ def get_targ_dens(targets, nside=256, bit_mask=None):
         log.info('Reading in target catalog...t = {:.1f}s'.format(time()-start))
         targets = fitsio.read(targets)
 
+    # ADM retrieve the bitmasks.
+    if Mx[0]._name == 'cmx_mask':
+        msg = 'generating target densities does NOT work for CMX files!!!'
+        log.critical(msg)
+        raise ValueError(msg)
+    else:
+        desi_mask, bgs_mask, mws_mask = Mx
+
     # ADM the number of pixels and the pixel area at the passed nside
     npix = hp.nside2npix(nside)
     pixarea = hp.nside2pixarea(nside, degrees=True)
@@ -701,26 +705,26 @@ def get_targ_dens(targets, nside=256, bit_mask=None):
 
     # ADM retrieve the bit names of interest
     from desitarget.QA import _load_targdens
-    bitnames = np.array(list(_load_targdens(bit_mask=bit_mask).keys()))
+    bitnames = np.array(list(_load_targdens(bit_mask=Mx).keys()))
 
     # ADM and set up an array to hold the output target densities
     targdens = np.zeros(npix, dtype=[(bitname, 'f4') for bitname in bitnames])
 
     for bitname in bitnames:
         if 'ALL' in bitname:
-            wbit = np.arange(len(targets))
+            ii = np.ones(len(targets)).astype('bool')
         else:
-            if ('BGS' in bitname) & ~('ANY' in bitname):
-                wbit = np.where(targets["BGS_TARGET"] & bgs_mask[bitname])[0]
-            elif ('MWS' in bitname) & ~('ANY' in bitname):
-                wbit = np.where(targets["MWS_TARGET"] & mws_mask[bitname])[0]
+            if ('BGS' in bitname) and not('S_ANY' in bitname):
+                ii = targets["BGS_TARGET"] & bgs_mask[bitname] != 0
+            elif ('MWS' in bitname) and not('S_ANY' in bitname):
+                ii = targets["MWS_TARGET"] & mws_mask[bitname] != 0
             else:
-                wbit = np.where(targets["DESI_TARGET"] & desi_mask[bitname])[0]
+                ii = targets["DESI_TARGET"] & desi_mask[bitname] != 0
 
-        if len(wbit) > 0:
+        if np.any(ii):
             # ADM calculate the number of objects in each pixel for the
             # ADM targets of interest
-            pixnum, pixcnt = np.unique(pixnums[wbit], return_counts=True)
+            pixnum, pixcnt = np.unique(pixnums[ii], return_counts=True)
             targdens[bitname][pixnum] = pixcnt/pixarea
 
     return targdens
@@ -787,7 +791,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
         cols = np.concatenate([["RA", "DEC"], targcols])
         targets = read_targets_in_box(targets, columns=cols)
     # ADM change target column names, and retrieve associated survey information.
-    _, bit_mask, survey, targets = main_cmx_or_sv(targets, rename=True)
+    _, Mx, survey, targets = main_cmx_or_sv(targets, rename=True)
 
     # ADM determine the areal coverage of the randoms at this nside.
     log.info('Determining footprint...t = {:.1f}s'.format(time()-start))
@@ -796,7 +800,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
 
     # ADM get the target densities.
     log.info('Calculating target densities...t = {:.1f}s'.format(time()-start))
-    targdens = get_targ_dens(targets, nside=nside, bit_mask=bit_mask)
+    targdens = get_targ_dens(targets, Mx, nside=nside)
 
     # ADM set up the output array.
     datamodel = [('HPXPIXEL', '>i4'), ('FRACAREA', '>f4'), ('STARDENS', '>f4'), ('EBV', '>f4'),

From d6a7c7dbe9b5cce6b690653666ecbbc0330a89ac Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Sat, 6 Apr 2019 06:33:42 -0700
Subject: [PATCH 16/31] bundling test now in new test_geomask module

---
 py/desitarget/test/test_geomask.py | 43 ++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100644 py/desitarget/test/test_geomask.py

diff --git a/py/desitarget/test/test_geomask.py b/py/desitarget/test/test_geomask.py
new file mode 100644
index 000000000..90c341d9f
--- /dev/null
+++ b/py/desitarget/test/test_geomask.py
@@ -0,0 +1,43 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+# -*- coding: utf-8 -*-
+"""Test desitarget.geomask.
+"""
+import unittest
+from pkg_resources import resource_filename
+import numpy as np
+import os
+
+from desitarget import geomask
+
+
+class TestGEOMASK(unittest.TestCase):
+
+    @classmethod
+    def setUp(self):
+        drdir = '/global/project/projectdirs/cosmo/work/legacysurvey/dr8b/'
+        self.surveydir = os.path.join(drdir, 'decam')
+        self.surveydir2 = os.path.join(drdir, '90prime-mosaic')
+
+    def test_bundle_bricks(self):
+        """
+        Test the bundle_bricks scripting code executes without bugs
+        """
+        blat = geomask.bundle_bricks(1, 1, 1,
+                                     surveydirs=[self.surveydir])
+        self.assertTrue(blat is None)
+
+        foo = geomask.bundle_bricks(1, 1, 1,
+                                    surveydirs=[self.surveydir, self.surveydir2])
+        self.assertTrue(foo is None)
+
+
+if __name__ == '__main__':
+    unittest.main()
+
+
+def test_suite():
+    """Allows testing of only this module with the command:
+
+        python setup.py test -m desitarget.test.test_geomask
+    """
+    return unittest.defaultTestLoader.loadTestsFromName(__name__)

From ec57cc817db0c967cbe1245b560e4bc017b7b039 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Sat, 6 Apr 2019 07:22:02 -0700
Subject: [PATCH 17/31] add PSFSIZEs to random catalogs and pixweight maps

---
 py/desitarget/randoms.py | 67 +++++++++++++++++++---------------------
 1 file changed, 31 insertions(+), 36 deletions(-)

diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py
index 8d89f96a7..55e2d4682 100644
--- a/py/desitarget/randoms.py
+++ b/py/desitarget/randoms.py
@@ -246,7 +246,7 @@ def dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
 
 
 def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
-    """Return NOBS, GALDEPTH, PSFDEPTH (per-band) at positions in one brick of the Legacy Surveys
+    """Observational quantities (per-band) at positions in a Legacy Surveys brick.
 
     Parameters
     ----------
@@ -263,8 +263,9 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
     Returns
     -------
     :class:`dictionary`
-       The number of observations (NOBS_X), PSF depth (PSFDEPTH_X) and Galaxy depth (GALDEPTH_X)
-       at each passed position in the Legacy Surveys in each band X. In addition, the MASKBITS
+       The number of observations (NOBS_X), PSF depth (PSFDEPTH_X)
+       Galaxy depth (GALDEPTH_X) and PSF size (PSFSIZE_X) at each
+       at each passed position in each band X. Plus, , the MASKBITS
        information at each passed position for the brick.
 
     Notes
@@ -291,9 +292,9 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
     for filt in ['g', 'r', 'z']:
         # ADM the input file labels, and output column names and output formats
         # ADM for each of the quantities of interest.
-        qnames = zip(['nexp', 'depth', 'galdepth'],
-                     ['nobs', 'psfdepth', 'galdepth'],
-                     ['i2', 'f4', 'f4'])
+        qnames = zip(['nexp', 'depth', 'galdepth', 'psfsize'],
+                     ['nobs', 'psfdepth', 'galdepth', 'psfsize'],
+                     ['i2', 'f4', 'f4', 'f4'])
         for qin, qout, qform in qnames:
             fn = os.path.join(
                 rootdir, 'legacysurvey-{}-{}-{}.fits.{}'.format(brickname, qin, filt, extn)
@@ -474,21 +475,15 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
     -------
     :class:`~numpy.ndarray`
         a numpy structured array with the following columns:
-            RA: Right Ascension of a random point
-            DEC: Declination of a random point
-            BRICKNAME: Passed brick name
-            NOBS_G: Number of observations at this location in the g-band
-            NOBS_R: Number of observations at this location in the r-band
-            NOBS_Z: Number of observations at this location in the z-band
-            PSFDEPTH_G: PSF depth at this location in the g-band
-            PSFDEPTH_R: PSF depth at this location in the r-band
-            PSFDEPTH_Z: PSF depth at this location in the z-band
-            GALDEPTH_G: Galaxy depth at this location in the g-band
-            GALDEPTH_R: Galaxy depth at this location in the r-band
-            GALDEPTH_Z: Galaxy depth at this location in the z-band
+            RA, DEC: Right Ascension, Declination of a random location.
+            BRICKNAME: Passed brick name.
+            NOBS_G, R, Z: Number of observations at this location in g, r, z-band.
+            PSFDEPTH_G, R, Z: PSF depth at this location in g, r, z.
+            GALDEPTH_G, R, Z: Galaxy depth in g, r, z.
+            PSFSIZE_G, R, Z: Weighted average PSF FWHM (arcsec) in g, r, z.
             MASKBITS: Extra mask bits info as stored in the header of e.g.,
               dr7dir + 'coadd/111/1116p210/legacysurvey-1116p210-maskbits.fits.gz'
-            EBV: E(B-V) at this location from the SFD dust maps
+            EBV: E(B-V) at this location from the SFD dust maps.
     """
     # ADM this is only intended to work on one brick, so die if a larger array is passed.
     if not isinstance(brickname, str):
@@ -516,6 +511,7 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
                             ('NOBS_G', 'i2'), ('NOBS_R', 'i2'), ('NOBS_Z', 'i2'),
                             ('PSFDEPTH_G', 'f4'), ('PSFDEPTH_R', 'f4'), ('PSFDEPTH_Z', 'f4'),
                             ('GALDEPTH_G', 'f4'), ('GALDEPTH_R', 'f4'), ('GALDEPTH_Z', 'f4'),
+                            ('PSFSIZE_G', 'f4'), ('PSFSIZE_R', 'f4'), ('PSFSIZE_Z', 'f4'),
                             ('MASKBITS', 'i2'), ('EBV', 'f4'), ('PHOTSYS', '|S1')])
     # ADM store each quantity of interest in the structured array
     # ADM remembering that the dictionary keys are in lower case text.
@@ -769,6 +765,8 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
                                 the median of PSFDEPTH values in the passed random catalog.
             - GALDEPTH_G, R, Z: The galaxy depth in g, r, z-band in the pixel, derived from
                                 the median of GALDEPTH values in the passed random catalog.
+            - PSFSIZE_G, R, Z: The weighted average PSF FWHM, in arcsec, in g, r, z in the pixel, 
+                               from the median of PSFSIZE values in the passed random catalog.
             - One column for every bit returned by :func:`desitarget.QA._load_targdens()`.
               Each column contains the density of targets in pixels at the passed `nside`
     :class:`str`
@@ -805,7 +803,8 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
     # ADM set up the output array.
     datamodel = [('HPXPIXEL', '>i4'), ('FRACAREA', '>f4'), ('STARDENS', '>f4'), ('EBV', '>f4'),
                  ('PSFDEPTH_G', '>f4'), ('PSFDEPTH_R', '>f4'), ('PSFDEPTH_Z', '>f4'),
-                 ('GALDEPTH_G', '>f4'), ('GALDEPTH_R', '>f4'), ('GALDEPTH_Z', '>f4')]
+                 ('GALDEPTH_G', '>f4'), ('GALDEPTH_R', '>f4'), ('GALDEPTH_Z', '>f4'),
+                 ('PSFSIZE_G', '>f4'), ('PSFSIZE_R', '>f4'), ('PSFSIZE_Z', '>f4')]
     datamodel += targdens.dtype.descr
     hpxinfo = np.zeros(npix, dtype=datamodel)
     # ADM set initial values to -1 so that they can easily be clipped.
@@ -843,8 +842,10 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
 
     # ADM work through the ordered pixels to populate the median for
     # ADM each quantity of interest.
-    cols = ['EBV', 'PSFDEPTH_G', 'GALDEPTH_G', 'PSFDEPTH_R', 'GALDEPTH_R',
-            'PSFDEPTH_Z', 'GALDEPTH_Z']
+    cols = ['EBV',
+            'PSFDEPTH_G', 'GALDEPTH_G', 'PSFSIZE_G',
+            'PSFDEPTH_R', 'GALDEPTH_R', 'PSFSIZE_R',
+            'PSFDEPTH_Z', 'GALDEPTH_Z', 'PSFSIZE_Z']
     for i in range(len(pixcnts)-1):
         inds = pixorder[pixcnts[i]:pixcnts[i+1]]
         pix = pixnums[inds][0]
@@ -902,21 +903,15 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None,
     -------
     :class:`~numpy.ndarray`
         a numpy structured array with the following columns:
-            RA: Right Ascension of a random point
-            DEC: Declination of a random point
-            BRICKNAME: Passed brick name
-            NOBS_G: Number of observations at this location in the g-band
-            NOBS_R: Number of observations at this location in the r-band
-            NOBS_Z: Number of observations at this location in the z-band
-            PSFDEPTH_G: PSF depth at this location in the g-band
-            PSFDEPTH_R: PSF depth at this location in the r-band
-            PSFDEPTH_Z: PSF depth at this location in the z-band
-            GALDEPTH_G: Galaxy depth at this location in the g-band
-            GALDEPTH_R: Galaxy depth at this location in the r-band
-            GALDEPTH_Z: Galaxy depth at this location in the z-band
+            RA, DEC: Right Ascension, Declination of a random location.
+            BRICKNAME: Passed brick name.
+            NOBS_G, R, Z: Number of observations at this location in g, r, z-band.
+            PSFDEPTH_G, R, Z: PSF depth at this location in g, r, z.
+            GALDEPTH_G, R, Z: Galaxy depth in g, r, z.
+            PSFSIZE_G, R, Z: Weighted average PSF FWHM (arcsec) in g, r, z.
             MASKBITS: Extra mask bits info as stored in the header of e.g.,
-              dr7dir + 'coadd/111/1116p210/legacysurvey-1116p210-maskbits.fits.gz'
-            EBV: E(B-V) at this location from the SFD dust maps
+              dr7dir + 'coadd/111/1116p210/legacysurvey-1116p210-maskbits.fits.gz'.
+            EBV: E(B-V) at this location from the SFD dust maps.
     """
     # ADM grab brick information for this data release. Depending on whether this
     # ADM is pre-or-post-DR8 we need to find the correct directory or directories.

From 15c3e6b1ac3556a919a2728d52554690108e7161 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 08:35:55 -0700
Subject: [PATCH 18/31] add sky locations to dictionary of meta-data

---
 py/desitarget/randoms.py | 55 +++++++++++++++++++++++++++++-----------
 1 file changed, 40 insertions(+), 15 deletions(-)

diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py
index 55e2d4682..743ac24f2 100644
--- a/py/desitarget/randoms.py
+++ b/py/desitarget/randoms.py
@@ -14,6 +14,7 @@
 from time import time
 import healpy as hp
 import fitsio
+import photutils
 from glob import glob
 from desitarget.gaiamatch import _get_gaia_dir
 from desitarget.geomask import bundle_bricks, box_area
@@ -245,7 +246,8 @@ def dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
     return qcombine
 
 
-def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
+def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir,
+                                       aprad=0.75):
     """Observational quantities (per-band) at positions in a Legacy Surveys brick.
 
     Parameters
@@ -259,7 +261,9 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
     drdir : :class:`str`
        The root directory pointing to a Data Release from the Legacy Surveys
        e.g. /global/project/projectdirs/cosmo/data/legacysurvey/dr7.
-
+    aprad : :class:`float`, optional, defaults to 0.75
+        Radii in arcsec of aperture for which to derive sky fluxes
+        defaults to the DESI fiber radius.
     Returns
     -------
     :class:`dictionary`
@@ -287,18 +291,17 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
     instrum = None
 
     rootdir = os.path.join(drdir, 'coadd', brickname[:3], brickname)
+    fileform = os.path.join(rootdir, 'legacysurvey-{}-{}-{}.fits.{}')
     # ADM loop through each of the filters and store the number of observations at the
     # ADM RA and Dec positions of the passed points.
     for filt in ['g', 'r', 'z']:
         # ADM the input file labels, and output column names and output formats
         # ADM for each of the quantities of interest.
-        qnames = zip(['nexp', 'depth', 'galdepth', 'psfsize'],
-                     ['nobs', 'psfdepth', 'galdepth', 'psfsize'],
-                     ['i2', 'f4', 'f4', 'f4'])
+        qnames = zip(['nexp', 'depth', 'galdepth', 'psfsize', 'image'],
+                     ['nobs', 'psfdepth', 'galdepth', 'psfsize', 'sky'],
+                     ['i2', 'f4', 'f4', 'f4', 'f4'])
         for qin, qout, qform in qnames:
-            fn = os.path.join(
-                rootdir, 'legacysurvey-{}-{}-{}.fits.{}'.format(brickname, qin, filt, extn)
-                )
+            fn = fileform.format(brickname, qin, filt, extn)
             # ADM only process the WCS if there is a file corresponding to this filter.
             if os.path.exists(fn):
                 img = fits.open(fn)
@@ -310,14 +313,36 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
                     iswcs = True
                 # ADM determine the quantity of interest at each passed location
                 # ADM and store in a dictionary with the filter and quantity name.
-                qdict[qout+'_'+filt] = img[extn_nb].data[y.astype("int"), x.astype("int")]
-                # log.info('Determined {} using WCS for {}...t = {:.1f}s'
-                #          .format(qout+'_'+filt,fn,time()-start))
+                if qout == 'sky':
+                    # ADM special treatment to photometer sky. Read in the ivar image.
+                    fnivar = fileform.format(brickname, 'invvar', filt, extn)
+                    ivar = fits.open(fnivar)
+                    with np.errstate(divide='ignore', invalid='ignore'):
+                        # ADM transform ivars to errors, guarding against 1/0.
+                        imsigma = 1./np.sqrt(ivar)
+                        imsigma[ivar == 0] = 0
+                    # ADM perform aperture photometry at the requested radius (aprad).
+                    apxy = np.vstack((x, y)).T
+                    aper = photutils.CircularAperture(apxy, aprad)
+                    p = photutils.aperture_photometry(img, aper, error=imsigma)
+                    # ADM store the results.
+                    qdict[qout+'apflux_'+filt] = p.field('aperture_sum') 
+                    err = p.field('aperture_sum_err')
+                    with np.errstate(divide='ignore', invalid='ignore'):
+                        # ADM transform errors to ivars, guarding against 1/0.
+                        ivar = 1./err**2.
+                        ivar[err == 0] = 0.
+                    qdict[qout+'apflux_ivar_'+filt] = ivar
+                else:
+                    qdict[qout+'_'+filt] = img[extn_nb].data[y.astype("int"), x.astype("int")]
+                
             else:
-                # log.info('no {} file at {}...t = {:.1f}s'
-                #          .format(qin+'_'+filt,fn,time()-start))
                 # ADM if the file doesn't exist, set the relevant quantities to zero.
-                qdict[qout+'_'+filt] = np.zeros(npts, dtype=qform)
+                if qout == 'sky':
+                    qdict[qout+'apflux_'+filt] = np.zeros(npts, dtype=qform)
+                    qdict[qout+'apflux_ivar_'+filt] = np.zeros(npts, dtype=qform)
+                else:
+                    qdict[qout+'_'+filt] = np.zeros(npts, dtype=qform)
 
     # ADM add the mask bits information.
     fn = os.path.join(rootdir,
@@ -765,7 +790,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
                                 the median of PSFDEPTH values in the passed random catalog.
             - GALDEPTH_G, R, Z: The galaxy depth in g, r, z-band in the pixel, derived from
                                 the median of GALDEPTH values in the passed random catalog.
-            - PSFSIZE_G, R, Z: The weighted average PSF FWHM, in arcsec, in g, r, z in the pixel, 
+            - PSFSIZE_G, R, Z: Weighted average PSF FWHM, in arcsec, in g, r, z in the pixel,
                                from the median of PSFSIZE values in the passed random catalog.
             - One column for every bit returned by :func:`desitarget.QA._load_targdens()`.
               Each column contains the density of targets in pixels at the passed `nside`

From fbb3b4cc0899187d88d1e4ba1d0123480780c1b8 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 10:55:11 -0700
Subject: [PATCH 19/31] debugging and finalizing writing of sky fluxes to
 randoms

---
 bin/select_randoms       |  6 ++-
 py/desitarget/randoms.py | 91 +++++++++++++++++++++-------------------
 2 files changed, 53 insertions(+), 44 deletions(-)

diff --git a/bin/select_randoms b/bin/select_randoms
index eda226348..0236808a4 100755
--- a/bin/select_randoms
+++ b/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,7 +92,7 @@ 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))
 
diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py
index 743ac24f2..d764955c1 100644
--- a/py/desitarget/randoms.py
+++ b/py/desitarget/randoms.py
@@ -216,7 +216,8 @@ def _pre_or_post_dr8(drdir):
     return drdirs
 
 
-def dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
+def dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir,
+                                           aprad=0.75):
     """Wrapper on `quantities_at_positions_in_a_brick` for DR8 imaging and beyond.
 
     Notes
@@ -233,7 +234,7 @@ def dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir):
     # ADM determine the dictionary of quantities for one or two directories.
     qall = []
     for dd in drdirs:
-        q = quantities_at_positions_in_a_brick(ras, decs, brickname, dd)
+        q = quantities_at_positions_in_a_brick(ras, decs, brickname, dd, aprad=aprad)
         # ADM don't count bricks where we never read a file header.
         if q is not None:
             qall.append(q)
@@ -267,16 +268,22 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir,
     Returns
     -------
     :class:`dictionary`
-       The number of observations (NOBS_X), PSF depth (PSFDEPTH_X)
-       Galaxy depth (GALDEPTH_X) and PSF size (PSFSIZE_X) at each
-       at each passed position in each band X. Plus, , the MASKBITS
+       The number of observations (`nobs_x`), PSF depth (`psfdepth_x`)
+       Galaxy depth (`galdepth_x`), PSF size (`psfsize_x`), and sky
+       background (`apflux_x`) and inverse variance (`apflux_ivar_x`)
+       at each passed position in each band X. Plus, the `maskbits`
        information at each passed position for the brick.
 
     Notes
     -----
         - First version copied shamelessly from Anand Raichoor.
     """
+    # ADM guard against too low a density of random locations.
     npts = len(ras)
+    if npts == 0:
+        msg = 'brick {} is empty. Increase the density of random points!'.format(brickname)
+        log.critical(msg)
+        raise ValueError(msg)
 
     # ADM determine whether the coadd files have extension .gz or .fz based on the DR directory.
     extn, extn_nb = dr_extension(drdir)
@@ -298,25 +305,25 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir,
         # ADM the input file labels, and output column names and output formats
         # ADM for each of the quantities of interest.
         qnames = zip(['nexp', 'depth', 'galdepth', 'psfsize', 'image'],
-                     ['nobs', 'psfdepth', 'galdepth', 'psfsize', 'sky'],
+                     ['nobs', 'psfdepth', 'galdepth', 'psfsize', 'apflux'],
                      ['i2', 'f4', 'f4', 'f4', 'f4'])
         for qin, qout, qform in qnames:
             fn = fileform.format(brickname, qin, filt, extn)
             # ADM only process the WCS if there is a file corresponding to this filter.
             if os.path.exists(fn):
-                img = fits.open(fn)
+                img = fits.open(fn)[extn_nb]
                 if not iswcs:
                     # ADM also store the instrument name, if it isn't yet stored.
-                    instrum = img[extn_nb].header["INSTRUME"].lower().strip()
-                    w = WCS(img[extn_nb].header)
+                    instrum = img.header["INSTRUME"].lower().strip()
+                    w = WCS(img.header)
                     x, y = w.all_world2pix(ras, decs, 0)
                     iswcs = True
                 # ADM determine the quantity of interest at each passed location
                 # ADM and store in a dictionary with the filter and quantity name.
-                if qout == 'sky':
+                if qout == 'apflux':
                     # ADM special treatment to photometer sky. Read in the ivar image.
                     fnivar = fileform.format(brickname, 'invvar', filt, extn)
-                    ivar = fits.open(fnivar)
+                    ivar = fits.open(fnivar)[extn_nb].data
                     with np.errstate(divide='ignore', invalid='ignore'):
                         # ADM transform ivars to errors, guarding against 1/0.
                         imsigma = 1./np.sqrt(ivar)
@@ -324,41 +331,38 @@ def quantities_at_positions_in_a_brick(ras, decs, brickname, drdir,
                     # ADM perform aperture photometry at the requested radius (aprad).
                     apxy = np.vstack((x, y)).T
                     aper = photutils.CircularAperture(apxy, aprad)
-                    p = photutils.aperture_photometry(img, aper, error=imsigma)
+                    p = photutils.aperture_photometry(img.data, aper, error=imsigma)
                     # ADM store the results.
-                    qdict[qout+'apflux_'+filt] = p.field('aperture_sum') 
+                    qdict[qout+'_'+filt] = np.array(p.field('aperture_sum'))
                     err = p.field('aperture_sum_err')
                     with np.errstate(divide='ignore', invalid='ignore'):
                         # ADM transform errors to ivars, guarding against 1/0.
                         ivar = 1./err**2.
                         ivar[err == 0] = 0.
-                    qdict[qout+'apflux_ivar_'+filt] = ivar
+                    qdict[qout+'_ivar_'+filt] = np.array(ivar)
                 else:
-                    qdict[qout+'_'+filt] = img[extn_nb].data[y.astype("int"), x.astype("int")]
-                
+                    qdict[qout+'_'+filt] = img.data[y.astype("int"), x.astype("int")]
+            # ADM if the file doesn't exist, set the relevant quantities to zero.
             else:
-                # ADM if the file doesn't exist, set the relevant quantities to zero.
-                if qout == 'sky':
-                    qdict[qout+'apflux_'+filt] = np.zeros(npts, dtype=qform)
-                    qdict[qout+'apflux_ivar_'+filt] = np.zeros(npts, dtype=qform)
-                else:
-                    qdict[qout+'_'+filt] = np.zeros(npts, dtype=qform)
+                if qout == 'apflux':
+                    qdict['apflux_ivar_'+filt] = np.zeros(npts, dtype=qform)
+                qdict[qout+'_'+filt] = np.zeros(npts, dtype=qform)
 
     # ADM add the mask bits information.
     fn = os.path.join(rootdir,
                       'legacysurvey-{}-maskbits.fits.{}'.format(brickname, extn))
     # ADM only process the WCS if there is a file corresponding to this filter.
     if os.path.exists(fn):
-        img = fits.open(fn)
+        img = fits.open(fn)[extn_nb]
         # ADM use the WCS calculated for the per-filter quantities above, if it exists.
         if not iswcs:
             # ADM also store the instrument name, if it isn't yet stored.
-            instrum = img[extn_nb].header["INSTRUME"].lower().strip()
-            w = WCS(img[extn_nb].header)
+            instrum = img.header["INSTRUME"].lower().strip()
+            w = WCS(img.header)
             x, y = w.all_world2pix(ras, decs, 0)
             iswcs = True
         # ADM add the maskbits to the dictionary.
-        qdict['maskbits'] = img[extn_nb].data[y.astype("int"), x.astype("int")]
+        qdict['maskbits'] = img.data[y.astype("int"), x.astype("int")]
     else:
         # ADM if there is no maskbits file, populate with zeros.
         qdict['maskbits'] = np.zeros(npts, dtype='i2')
@@ -471,7 +475,7 @@ def get_dust(ras, decs, scaling=1, dustdir=None):
 
 
 def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
-                              density=100000, dustdir=None):
+                              density=100000, dustdir=None, aprad=0.75):
     """NOBS, DEPTHS etc. (per-band) for random points in a brick of the Legacy Surveys
 
     Parameters
@@ -495,6 +499,9 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
     dustdir : :class:`str`, optional, defaults to $DUST_DIR+'/maps'
         The root directory pointing to SFD dust maps. If not
         sent the code will try to use $DUST_DIR+'/maps' before failing.
+    aprad : :class:`float`, optional, defaults to 0.75
+        Radii in arcsec of aperture for which to derive sky fluxes
+        defaults to the DESI fiber radius.
 
     Returns
     -------
@@ -506,6 +513,8 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
             PSFDEPTH_G, R, Z: PSF depth at this location in g, r, z.
             GALDEPTH_G, R, Z: Galaxy depth in g, r, z.
             PSFSIZE_G, R, Z: Weighted average PSF FWHM (arcsec) in g, r, z.
+            APFLUX_G, R, Z: Sky background extracted in `aprad` in g, r, z.
+            APFLUX_IVAR_G, R, Z: Inverse variance of sky background in g, r, z.
             MASKBITS: Extra mask bits info as stored in the header of e.g.,
               dr7dir + 'coadd/111/1116p210/legacysurvey-1116p210-maskbits.fits.gz'
             EBV: E(B-V) at this location from the SFD dust maps.
@@ -519,7 +528,7 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
     ras, decs = randoms_in_a_brick_from_edges(ramin, ramax, decmin, decmax, density=density)
 
     # ADM retrieve the dictionary of quantities for each random point.
-    qdict = dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir)
+    qdict = dr8_quantities_at_positions_in_a_brick(ras, decs, brickname, drdir, aprad=aprad)
 
     # ADM if we detected different cameras corresponding to 2 Data Releases
     # ADM for this brick then we need to duplicate the ras, decs as well.
@@ -537,6 +546,8 @@ def get_quantities_in_a_brick(ramin, ramax, decmin, decmax, brickname, drdir,
                             ('PSFDEPTH_G', 'f4'), ('PSFDEPTH_R', 'f4'), ('PSFDEPTH_Z', 'f4'),
                             ('GALDEPTH_G', 'f4'), ('GALDEPTH_R', 'f4'), ('GALDEPTH_Z', 'f4'),
                             ('PSFSIZE_G', 'f4'), ('PSFSIZE_R', 'f4'), ('PSFSIZE_Z', 'f4'),
+                            ('APFLUX_G', 'f4'), ('APFLUX_R', 'f4'), ('APFLUX_Z', 'f4'),
+                            ('APFLUX_IVAR_G', 'f4'), ('APFLUX_IVAR_R', 'f4'), ('APFLUX_IVAR_Z', 'f4'),
                             ('MASKBITS', 'i2'), ('EBV', 'f4'), ('PHOTSYS', '|S1')])
     # ADM store each quantity of interest in the structured array
     # ADM remembering that the dictionary keys are in lower case text.
@@ -691,7 +702,7 @@ def get_targ_dens(targets, Mx, nside=256):
         A corresponding (same Legacy Surveys Data Release) target catalog as made by,
         e.g., :func:`desitarget.cuts.select_targets()`, or the name of such a file.
     Mx : :class:`list` or `~numpy.array`
-        The targeting bitmasks associated with the passed targets, assumed to be 
+        The targeting bitmasks associated with the passed targets, assumed to be
         a desi, bgs and mws mask in that order (for either SV or the main survey).
     nside : :class:`int`, optional, defaults to nside=256 (~0.0525 sq. deg. or "brick-sized")
         The resolution (HEALPixel nside number) at which to build the map (NESTED scheme).
@@ -884,7 +895,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None):
 
 def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None,
                    bundlebricks=None, brickspersec=2.5,
-                   dustdir=None, resolverands=True):
+                   dustdir=None, resolverands=True, aprad=0.75):
     """NOBS, DEPTHs, MASKBITs (per-band) for random points in a Legacy Surveys DR.
 
     Parameters
@@ -923,20 +934,16 @@ def select_randoms(drdir, density=100000, numproc=32, nside=4, pixlist=None,
     resolverands : :class:`boolean`, optional, defaults to ``True``
         If ``True``, resolve randoms into northern randoms in northern regions
         and southern randoms in southern regions.
+    aprad : :class:`float`, optional, defaults to 0.75
+        Radii in arcsec of aperture for which to derive sky fluxes
+        defaults to the DESI fiber radius.
 
     Returns
     -------
     :class:`~numpy.ndarray`
-        a numpy structured array with the following columns:
-            RA, DEC: Right Ascension, Declination of a random location.
-            BRICKNAME: Passed brick name.
-            NOBS_G, R, Z: Number of observations at this location in g, r, z-band.
-            PSFDEPTH_G, R, Z: PSF depth at this location in g, r, z.
-            GALDEPTH_G, R, Z: Galaxy depth in g, r, z.
-            PSFSIZE_G, R, Z: Weighted average PSF FWHM (arcsec) in g, r, z.
-            MASKBITS: Extra mask bits info as stored in the header of e.g.,
-              dr7dir + 'coadd/111/1116p210/legacysurvey-1116p210-maskbits.fits.gz'.
-            EBV: E(B-V) at this location from the SFD dust maps.
+        a numpy structured array with the same columns as returned by
+        :func:`~desitarget.randoms.get_quantities_in_a_brick`.
+
     """
     # ADM grab brick information for this data release. Depending on whether this
     # ADM is pre-or-post-DR8 we need to find the correct directory or directories.
@@ -989,8 +996,8 @@ def _get_quantities(brickname):
 
         # ADM populate the brick with random points, and retrieve the quantities
         # ADM of interest at those points.
-        return get_quantities_in_a_brick(bramin, bramax, bdecmin, bdecmax, brickname,
-                                         drdir, density=density, dustdir=dustdir)
+        return get_quantities_in_a_brick(bramin, bramax, bdecmin, bdecmax, brickname, drdir,
+                                         density=density, dustdir=dustdir, aprad=aprad)
 
     # ADM this is just to count bricks in _update_status
     nbrick = np.zeros((), dtype='i8')

From 9030834e1c07ca63cfc914d7b51dcbbfede98aff Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 11:20:00 -0700
Subject: [PATCH 20/31] add aperture radius for photometering sky positions to
 output file header

---
 bin/select_randoms  | 2 +-
 py/desitarget/io.py | 9 +++++++--
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/bin/select_randoms b/bin/select_randoms
index 0236808a4..1faa67c5b 100755
--- a/bin/select_randoms
+++ b/bin/select_randoms
@@ -97,7 +97,7 @@ randoms = select_randoms(ns.surveydir, density=ns.density, numproc=ns.numproc,
                          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))
diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index c4d1401a6..f5574733a 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -700,7 +700,7 @@ def write_gfas(filename, data, indir=None, indir2=None, nside=None,
 
 
 def write_randoms(filename, data, indir=None, hdr=None, nside=None,
-                  density=None, resolve=True):
+                  density=None, resolve=True, aprad=None):
     """Write a catalogue of randoms and associated pixel-level information.
 
     Parameters
@@ -722,7 +722,8 @@ def write_randoms(filename, data, indir=None, hdr=None, nside=None,
         write to header of the output file if not None.
     resolve : :class:`bool`, optional, defaults to ``True``
         Written to the output file header as `RESOLVE`.
-
+    aprad : :class:`float, optional, defaults to ``None``
+        If passes, written to the outpue header as `APRAD`.
     """
     # ADM create header to include versions, etc. If a `hdr` was
     # ADM passed, then use it, if not then create a new header.
@@ -751,6 +752,10 @@ def write_randoms(filename, data, indir=None, hdr=None, nside=None,
     if density is not None:
         hdr['DENSITY'] = density
 
+    # ADM add aperture radius (in arcsec) if requested by input.
+    if aprad is not None:
+        hdr['APRAD'] = aprad
+
     # ADM add whether or not the randoms were resolved to the header.
     hdr["RESOLVE"] = resolve
 

From c53bd39da8131a65546c4254f1567bc5b12f0a67 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 12:43:00 -0700
Subject: [PATCH 21/31] parallelize plot production in QA

---
 bin/run_target_qa   |  9 ++++++++-
 py/desitarget/QA.py | 47 ++++++++++++++++++++++++++++++++++++---------
 py/desitarget/io.py |  2 +-
 3 files changed, 47 insertions(+), 11 deletions(-)

diff --git a/bin/run_target_qa b/bin/run_target_qa
index 62a536b54..7dcf2eaa0 100755
--- a/bin/run_target_qa
+++ b/bin/run_target_qa
@@ -12,6 +12,9 @@ from desitarget.io import read_targets_in_box_header
 from desiutil.log import get_logger
 log = get_logger()
 
+import multiprocessing
+nproc = multiprocessing.cpu_count() // 2
+
 from argparse import ArgumentParser
 ap = ArgumentParser()
 ap.add_argument("src", 
@@ -30,6 +33,9 @@ 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()
 
@@ -64,7 +70,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))
diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index e8a948b80..667d46a08 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -28,6 +28,7 @@
 from desiutil import brick
 from desiutil.log import get_logger
 from desiutil.plots import init_sky, plot_sky_binned, plot_healpix_map, prepare_data
+from desitarget.internal import sharedmem
 from desitarget.targetmask import desi_mask, bgs_mask, mws_mask
 from desitarget.targets import main_cmx_or_sv
 from desitarget.io import read_targets_in_box, target_columns_from_header
@@ -1514,7 +1515,7 @@ def _in_desi_footprint(targs):
 
 def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True,
                   imaging_map_file=None, truths=None, objtruths=None, tcnames=None,
-                  cmx=False, bit_mask=None, mocks=False):
+                  cmx=False, bit_mask=None, mocks=False, numproc=32):
     """Make DESI targeting QA plots given a passed set of targets.
 
     Parameters
@@ -1552,6 +1553,8 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
         constraints) instead of the main survey bits.
     mocks : :class:`boolean`, optional, default=False
         If ``True``, add plots that are only relevant to mocks at the bottom of the webpage.
+    numproc : :class:`int`, optional, defaults to 32
+        The number of parallel processes to use to generate plots.
 
     Returns
     -------
@@ -1564,7 +1567,6 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
           target densities.
         - On execution, a set of .png plots for target QA are written to `qadir`.
     """
-
     # ADM set up the default logger from desiutil.
     log = get_logger()
 
@@ -1656,7 +1658,19 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
                       'MWS_ANY': 2000, 'MWS_BROAD': 2000, 'MWS_WD': 50, 'MWS_NEARBY': 50,
                       'MWS_MAIN_RED': 2000, 'MWS_MAIN_BLUE': 2000}
 
-    for objtype in targdens:
+    nbits = len(targdens)
+    nbit = np.ones((), dtype='i8')
+    t0 = time()
+    def _update_status(result):
+        """wrapper function for the critical reduction operation,
+        that occurs on the main parallel process"""
+        log.info('Done {}/{} bitnames...t = {:.1f}s'.format(nbit, nbits, time()-t0))
+        nbit[...] += 1    # this is an in-place modification                                                                                                                                                
+        return result
+
+    def _generate_plots(objtype):
+        """Make relevant plots for each bit name in objtype"""
+
         if 'ALL' in objtype:
             ii = np.ones(len(targs)).astype('bool')
         else:
@@ -1714,13 +1728,21 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
                 log.info('Made Gaia-based plots for {}...t = {:.1f}s'
                          .format(objtype, time()-start))
 
+    if numproc > 1:
+        pool = sharedmem.MapReduce(np=numproc)
+        with pool:
+            pool.map(_generate_plots, list(targdens.keys()), reduce=_update_status)
+    else:
+        for objtype in targdens:
+            _update_status(_generate_plots(objtype))
+
     log.info('Made QA plots...t = {:.1f}s'.format(time()-start))
     return totarea
 
 
 def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.',
                  clip2foot=False, weight=True, imaging_map_file=None,
-                 tcnames=None, systematics=True):
+                 tcnames=None, systematics=True, numproc=32):
     """Create a directory containing a webpage structure in which to embed QA plots.
 
     Parameters
@@ -1755,6 +1777,8 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
         for those specific bits. A useful speed-up when testing
     systematics : :class:`boolean`, optional, defaults to ``True``
         If sent, then add plots of systematics to the front page.
+    numproc : :class:`int`, optional, defaults to 32
+        The number of parallel processes to use to generate plots.
 
     Returns
     -------
@@ -1981,11 +2005,11 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
     # ADM make the QA plots, if requested:
     if makeplots:
         if svs == "MAIN":
-            totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths,
+            totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths, numproc=numproc,
                                     qadir=qadir, targdens=targdens, max_bin_area=max_bin_area,
                                     weight=weight, imaging_map_file=imaging_map_file, mocks=mocks)
         else:
-            totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths,
+            totarea = make_qa_plots(targs, truths=truths, objtruths=objtruths, numproc=numproc,
                                     qadir=qadir, targdens=targdens, max_bin_area=max_bin_area,
                                     weight=weight, imaging_map_file=imaging_map_file,
                                     cmx=cmx, bit_mask=masks, mocks=mocks)
@@ -2000,17 +2024,22 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
         # ADM write out a list of the target categories.
         headerlist = list(settargdens)
         headerlist.sort()
+        # ADM edit SV target categories to their initial letters to squeeze space.
+        hl = headerlist.copy()
+        if svs == 'SV':
+            for tc in 'QSO', 'ELG', 'LRG', 'BGS', 'MWS':
+                hl = [h.replace(tc, tc[0]) for h in hl]
         # ADM truncate the bit names at "trunc" characters to pack them more easily.
-        trunc = 9
+        trunc = 8
         truncform = '{:>'+str(trunc)+'s}'
-        headerwrite = [bitname[:trunc] for bitname in headerlist]
+        headerwrite = [bitname[:trunc] for bitname in hl]
         headerwrite.insert(0, " ")
         header = " ".join([truncform.format(i) for i in headerwrite])+'\n\n'
         htmlmain.write(header)
         # ADM for each pair of target classes, determine how many targets per unit area
         # ADM have the relevant target bit set for both target classes in the pair.
         for i, objtype1 in enumerate(headerlist):
-            overlaps = [objtype1[:trunc]]
+            overlaps = [hl[i][:trunc]]
             for j, objtype2 in enumerate(headerlist):
                 if j < i:
                     overlaps.append(" ")
diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index f5574733a..316c73738 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -723,7 +723,7 @@ def write_randoms(filename, data, indir=None, hdr=None, nside=None,
     resolve : :class:`bool`, optional, defaults to ``True``
         Written to the output file header as `RESOLVE`.
     aprad : :class:`float, optional, defaults to ``None``
-        If passes, written to the outpue header as `APRAD`.
+        If passed, written to the outpue header as `APRAD`.
     """
     # ADM create header to include versions, etc. If a `hdr` was
     # ADM passed, then use it, if not then create a new header.

From 2180dc1e83ae6d3c78ce0cccd07f18efdd8f2c65 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 12:47:56 -0700
Subject: [PATCH 22/31] add PSFSIZE plots to QA output pages

---
 py/desitarget/QA.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 667d46a08..4bd027157 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -109,6 +109,9 @@ def _load_systematics():
     sysdict['GALDEPTH_G'] = [63., 6300., 'Galaxy Depth in g-band']
     sysdict['GALDEPTH_R'] = [25., 2500., 'Galaxy Depth in r-band']
     sysdict['GALDEPTH_Z'] = [4., 400., 'Galaxy Depth in z-band']
+    sysdict['PSFSIZE_G'] = [0., 3., 'PSF Size in g-band']
+    sysdict['PSFSIZE_R'] = [0., 3., 'PSF Size in r-band']
+    sysdict['PSFSIZE_Z'] = [0., 3., 'PSF Size in z-band']
 
     return sysdict
 
@@ -1664,8 +1667,8 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
     def _update_status(result):
         """wrapper function for the critical reduction operation,
         that occurs on the main parallel process"""
-        log.info('Done {}/{} bitnames...t = {:.1f}s'.format(nbit, nbits, time()-t0))
-        nbit[...] += 1    # this is an in-place modification                                                                                                                                                
+        log.info('Done {}/{} bit names...t = {:.1f}s'.format(nbit, nbits, time()-t0))
+        nbit[...] += 1    # this is an in-place modification.
         return result
 
     def _generate_plots(objtype):

From 900fd6075a20c5988d4ff2a15cb773e698047c90 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 15:30:08 -0700
Subject: [PATCH 23/31] update data files for unit tests

---
 py/desitarget/QA.py                 |   1 +
 py/desitarget/io.py                 |   2 +-
 py/desitarget/test/t/pixweight.fits | Bin 8640 -> 8640 bytes
 py/desitarget/test/t/targets.fits   |  14 +++++++-------
 4 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 4bd027157..52800be8c 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -1664,6 +1664,7 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
     nbits = len(targdens)
     nbit = np.ones((), dtype='i8')
     t0 = time()
+
     def _update_status(result):
         """wrapper function for the critical reduction operation,
         that occurs on the main parallel process"""
diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index 316c73738..36623eb91 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -341,7 +341,7 @@ def read_tractor(filename, header=False, columns=None):
         for col in dr7datamodel.dtype.names:
             readcolumns.append(col)
         # ADM deal with some custom files I made that don't contain WISEMASK.
-        if 'WISEMASK_W1' not in fxcolnames:
+        if ('WISEMASK_W1' not in fxcolnames) and ('wisemask_w1' not in fxcolnames):
             readcolumns.remove('WISEMASK_W1')
             readcolumns.remove('WISEMASK_W2')
     # ADM if BRIGHTBLOB exists (it does for DR8, not for DR7) add it and
diff --git a/py/desitarget/test/t/pixweight.fits b/py/desitarget/test/t/pixweight.fits
index 9c8851de28cda55a515ae93ea2e4fd717ce6f75f..72a2b8a763785d01fd19ba4637122d7cc4fed15d 100644
GIT binary patch
delta 199
zcmX@$e872w4-b>6(Plp$dqzg%&3(LgS-Ap&-GV)%T;tulq6@Obh}vFtbh!KLG(SvyKhj1(Qn^D3g;DZIj>?pt1ts0keP>fCsZc
PAbK??.?EM޵?S?a"?n??}8S?~J@u#^{
 ngɃX>?%o?@v5@6?S?a#?n@=?}8i?~J@u"RtX0yw@BC~X@@΁A-A+AvN?Sm?a?n2?}6:?~H@H@@u"~hȿxGuX@z?AW#AۯXBgA(l?St?a?n5?}6?~I@u#"~\|X>}!>I?^k@s>?;,B365?SL?a?n<8?}7?~I@u#SpizKX>e?El?1A9ZA@{?S?a?n;?}7?~I@u"p\X>>z?é@6w?F=?Su?ak?n6?}6?~I
@u#i]i5AdX?X?>?@@b?Sb?a?n? '?}??S{4?ag?n8\?}7#?~II@u }	ۿkX>l>"?aؾʾV?S.?`+?n?}.?~D@u ]i~߹lX>z>i?@[e@7@)[?R?`?n.?}.?~C@u ΙbfAoX?V@xEAyB$TA.?Rݻ?`s?mv?},T?~B	@u"|	bv#XABkAYBBrBH[zB!a?St?a?n5?}6?~I@u"襞Mc0
 X?>1A(A7B$`A€?Su?a?n6?}6?~I	@u!F
N@XC8C
DJD8C}?S ?`u?n?}0?~Et@u!Mҿ<]XX@SAn%B4Bg`µ?Ԉ?z*
W?S,?`?n?}1?~E@u">S?/tIX@
-AVkAkAWA
-?So?a
?n35?}6N?~H@u#ʰRX>?wAB"$A?S{a?a?n8o?}7&?~IK	@u#7I/LqihlX>>J?r?lf?Sq$?a0?n4?}6s?~H@u#UO엓|X>J>}?Ohڰоda\?S{?aQ?n8O?}7!?~IH@u!d7޿K:hXA(>A>"3AuBjBZ?S@rwn@u"pKڿJ)ؐX>Ч>	4?W4@RL>?SY6?`?n)?}4?~G@u 1'IXJqX?	Ҝ?/??W?R?`?m?}+?~B)@u"^uBEBYBX?S.?`n?n?}1?~F@u!̓
+AVkAkAWA
-?So?a
?n35?}6N?~H@u#ʰRX>?wAB"$A?S{a?a?n8o?}7&?~IK	@u#7I/LqihlX>>J?r?lf?Sq$?a0?n4?}6s?~H@u#UO엓|X>J>}?Ohڰоda\?S{?aQ?n8O?}7!?~IH@u!d7޿K:hXA(>A>"3AuBjBZ?S@rwn@u"pKڿJ)ؐX>Ч>	4?W4@RL>?SY6?`?n)?}4?~G@u 1'IXJqX?	Ҝ?/??W?R?`?m?}+?~B)@u!}q ;2#X@oA_B.jBBB" ?S,?`ݬ?n?}1?~E@u!522CX@PYA{>uBEBYBX?S.?`n?n?}1?~F@u!̓
 /6 XAaA׶"B$GAA?S.?`߫?n?}1?~F@u 3EzqXA6ABAaAIZZ?RR?`O?m?},?~B@u 	A'b Q-X?F?@9_@+@x8?R]?`|?m}?}(t?~@;@u#Hn+}X>??3y?;=l?Sc=?a8?n.0?}5?~HF@u!8'E#sXA~AB~oB=
B
 i?R(?`n?n?}.?~D	@u"E%79X?N#@hUAJB:jAu?SU?`?n(0?}4?~G	@u#XS$2SKX??0A?@fk@sn?SZ?`?n*?}4?~G@u"Rʇ íOX?'?y@]`
 "?SM,?`e?n$?}3?~GW@u 32ܿk^BX@AD"B&'BoBq7&?R?`x?m?}(
?~?	@u tٿ
@@ -40,7 +40,7 @@ xdX>
 A
A?A|R?YC?e7X?p?}Z?~@u"yFS0{X?{?%v???65?Y?e?p?}+?~F@u#TLx@X>z>3?fҿSpP?Yk?eS?p?}?~@u 0F߼BݰX>%?<$??L?V?W?d?p?}?~v@u _,*ځXۄXA'гAArmAAA03?YX?e?p?}?~[@u 䦹'ڔX@ж.@ATAͲB/6?Y?e?ph?}?~c>ZSR@u!.9/X?&"?`??L?Y?e$?p?}?~(@u ®ą
XAǶAŢBM?@B.AW?Y?eN?pS?}?~ @u#>`,qX>>?2z:?Y?>?e4?p?}
 ?~@u _?ʗ"X>@A=+5B49A?Y?es?p;?}?~{	@u"n'lX>H?&?C?Z@?Y?ep?pb?}?~i@u o91pX>G}>? ?@
FD@?Y?e1?pn?}?~b@u"Ér'}ɡX>^R>>?}M?_t?Y\?eq?p?}?~w@u"MyvQ+ZX>>?su/@HA@o,?Y?e}?p?}?~&@u""_@XX>>Ꞹ?c@K@3?Y?eu?p?}(?~@u zF`Kz:XASTABAA(?YL?e=?pJ?}?~q@u"aaX?ڱ?Q@?HJY?Y?ecJ?p[?}|?~@u#$݇@X>?b?
 =]>پ?Yw.?e[?p?}?~/@u #(e:X?B[?
-@E7?Y?etF?pS?}?~@u Ӆ%w
qX>?}??]O?g?Yx?et?p˥?}!?~@u"gzkX>/?+?9@D@?Y?ez?p.?}?~@u#pTdԾX>O>?eݾDS|@R?Y/?e)"?p4?}?~@@u#f- /F4X?xQ?j?6?h??Y?e?p?}E?~@u!SD2NymX>l>(?q<@AF?Yz?e}Y?pЧ?}?~!@u#+~#sa`XC}DDDD??Yj?eS	?p?}?~@u T2X?AA/BB,t?YQ?er?pT?}?~	@u T-DݗX@]A1TAnP?&??@?Y?es?p?}?~@u"<Su3X?I?$t?qZ?+?L?Y?e|?p?}?~
@u 
P(X@AuB)رB|͎B8$?Y?ek?pe?}J?~@u#BoWX>i`@5xAXB{AjI?X?e
?p?}??~	@u 7XX>9>S?@&X?Y?eb|?p?}h?~@u#TB@X@&>AB0BHA?Y#?e _?p?}1?~@u"D :+XAeABKIB[nBR?Y?ez>?p?}?~@u ϙRw/+X>R'>L1?> _?Yz?eh?pď?}?~@u#
+@E7?Y?etF?pS?}?~@u Ӆ%w
qX>?}??]O?g?Yx?et?p˥?}!?~@u"gzkX>/?+?9@D@?Y?ez?p.?}?~@u#pTdԾX>O>?eݾDS|@R?Y/?e)"?p4?}?~@@u#f- /F4X?xQ?j?6?h??Y?e?p?}E?~@u!SD2NymX>l>(?q<@AF?Yz?e}Y?pЧ?}?~!@u#+~#sa`XC}DDDD??Yj?eS	?p?}?~@u T2X?AA/BB,t?YQ?er?pT?}?~	@u T-DݗX@]A1TAnP?&??@?Y?es?p?}?~@u 
P(X@AuB)رB|͎B8$?Y?ek?pe?}J?~@u#BoWX>i`@5xAXB{AjI?X?e
?p?}??~	@u 7XX>9>S?@&X?Y?eb|?p?}h?~@u#TB@X@&>AB0BHA?Y#?e _?p?}1?~@u"D :+XAeABKIB[nBR?Y?ez>?p?}?~@u ϙRw/+X>R'>L1?> _?Yz?eh?pď?}?~@u#
 0:/>Q_X?C?bR?޴@~N@?Y\?e?p+?}?~_@u#V ItXA$AAz5B#B+!?Y&?e"?pX?}f?~@u#j{@GX?g?@l}@Ǧ@[?Y(?e#?p(?}?~@u#X9X>?>D?qX4G+?Y"?e?p?}&?~@u"xG
 XAOX@AyNAצAAEk?Y?ecW?pb?}}?~@u #"\rX>??ӽQU?XEU?dD?p=#?}M?~|@u l9hXX?+Qw?;{T?;nh>?@[B@RF ?Yo?eV?p?}I?~@u#<X>:{?
 Π?G5HIc?Yl?eTm?p?}?~@u"-x7n2~fmX?6?d$?@l{-0&]?Y-?ev?p[?}>?~@u!^a2X>5>?6@@D`?Y?ep"?p?}?~]@u"-h7\>)X>[>cS?_ϭ>w?Y?en?p?}?~H@u K)ZmXAPaABskBFA6?Y\a?eH?p?}?~@u 
kŪ>X?f?@H@l?l?W_?dNS?p?}g?~y@u"{>+eX>/?"?!^ft?YP?ed_?p?}?~@u#1'5X>e?t.?@@?Y#/?e ?p?}*?~@u#+>Q\0X>T?`?̇Ͻf?X^?d9?p=?}?~@u!A54y%X?,$?[?:?@/x?Yjd?eR?p?}?~@u#L(X>6>Y$?@eN@Y?YXZ?eE?p?}?~@u#J}a匤XA:Ag%A?ADv?z4?Y;?e1{?p?}?~@u"6-D9X?De>=?B???i?YB?eb3?p?}a?~@u#ΎexQpX>r}>ь/?#,D?[U?Xs?d?pB?}?~@u!iw.X?~?fA
@@ -105,22 +105,22 @@ p
 ?fZ?nm?v?~u'?
;@u&16(
 X>?mk?o4?ˋ??fI?na?u{?~t	?@u#ޱpJOXAYk9A`B6AeA6?eu?m?uP?~f?L@u#}0;X?u@LAB\)AZ?e?m"?u?~j?	@u$9k6#|X?H??O? @#d@
 ?e?n\?ui?~k?@u$J1ӳ@X>z>)?{@'n1F?e?mQ?uM?~j?@u$%:tGX@AOAAA#?e?n?uf?~k?O@u%넬ǪX@h-@%A=8AijAb?eЏ?n
u?u?~lg?	@u'rnZ|X?6D?2`?D@@A2&?f?n2?us?~o?	@u#tok(0 X>r?t??@Υ?e?m?u?~`?@u%8̀[X>$>)?d@LU@|p?e?m*?u?~g?@u$ZgdvTUX@X;?u??!??e?m?u?~k&?@u#$vTwcX?a	?=@@a@?d?m3?uN?~X?~@u%+ X@F<]ASBmBF3?e?nN?uw?~m?H	@u$Ag,#JX?AA)A=V?ea?mf?u?~j}? @u$bșQqX@A	
-BFAA?eB?mٚ?u?~g?	@u%O}׮ƾLX@Ah;B(B)B6 ?e?me?uH?~g?@u%ؘ`X??Au:AEBJ3A?eI?m?u?~h?#	@u#oX@(@tmR@LA	A|??e#?m?u@?~ap?@u#Zd8X?'G6?3?c@}8?olL?e3z?m?u?~bo?@u$G#b8KaX@vAB&/BB+ܩ?e?m|?u?~g?	@u$
GoA؋9X>?\W?]@,Q?ui?e{?m6?u?~g?@u&QԸ՟=X?z@jxABn!B-?e	?n&?uۡ?~n?	D	@u$OtwwdPX?6+?k?@8I{?Z?eUK?mk?u:?~d?@u&bҦX>a?!]?j2?f+?nL?u?~r*?d@u#yCUsvX>$7?-u??j O?eQ-?m?u?~dR?@u&2pX@<ArB?BrPB)?f?n0l?u ?~o?		@u&KuydCX?L?@B_@+~@}?ez?n	?u?~l?@u$ȶ-X?fǟ@)@(0BBC(D?en?mȴ?u7?~f'?@u%ѤO^b?XAJ}B|BbB{(B.2?ei?m3?u0?~e?@u&VsjXAbE"A^B"AAA?e?n
+BFAA?eB?mٚ?u?~g?	@u%O}׮ƾLX@Ah;B(B)B6 ?e?me?uH?~g?@u%ؘ`X??Au:AEBJ3A?eI?m?u?~h?#	@u#oX@(@tmR@LA	A|??e#?m?u@?~ap?@u#Zd8X?'G6?3?c@}8?olL?e3z?m?u?~bo?@u$G#b8KaX@vAB&/BB+ܩ?e?m|?u?~g?	@u$
GoA؋9X>?\W?]@,Q?ui?e{?m6?u?~g?@u&QԸ՟=X?z@jxABn!B-?e	?n&?uۡ?~n?	D	@u$OtwwdPX?6+?k?@8I{?Z?eUK?mk?u:?~d?@u&bҦX>a?!]?j2?f+?nL?u?~r*?d@u#yCUsvX>$7?-u??j O?eQ-?m?u?~dR?@u&2pX@<ArB?BrPB)?f?n0l?u ?~o?		@u&KuydCX?L?@B_@+~@}?ez?n	?u?~l?@u%ѤO^b?XAJ}B|BbB{(B.2?ei?m3?u0?~e?@u&VsjXAbE"A^B"AAA?e?n
 L?u?~l?@u&ueU_X>a>?;|QQ?f
?n87?u?~pI?
 <@u#]=ê:WX?-r@AaBbAO#?e6?m%?u?~b?	@u#̓T&䋈(gX?R3Ȓ>3X?/@0@ځ?e`?ml?u?~eO?z@u# 6-UX?:@iA9=B4gAպ?d?m:#?uR?~Y-?~	@u$<=X?~?0?ُV@?N?eD?m?u?~c?c@u&_g3TX>9>?-@SA]X?e?n?u?~ks?B@u%X"M14X?d@Az/B=*A?e?m??u?~gz?	@u#
=CNXAIAA>NAאA?eb?m?u?~a?@u$q;IO#@fX>I>?A	V@C`M?e,?m?u?~a?j@u&#FO`uX>g?-m?׼D_rX?e?n>?u?~kK?)@u#
 s]IX?Tx?(?e@g>	@?d!?m.`?uK?~X?~Y@u&"Dn̆?X@@ADSAAC]@c?e?m?u%?~j?@u&<JX?I?c
-?A9bArL?e?n?uɂ?~k?@u&u
EX@UGAAzA$B1CAA"+?e?n?u?~k?X@u#o3kX@A1lAȼHBgBf?d?msD?us?~^a?~6@u%Ӟ&LX?O"@FAYB*A$9?ect?mN?u?~e{?	@u%Wrg#kX?.#?cv!??)?,~?ev?m?u?~f?X@u$ bU4{^X>:?RB?>A:AK?e?mf?u@?~a?@u$)b.6̶X>CJ>$?gsFԾ?d?mo2?uqz?~^?~@u&aƐ`OF	X@ ABwnB1A?e?n?u>?~kY?2@u$JcбXX>_>?\@@(A?e
+?A9bArL?e?n?uɂ?~k?@u&u
EX@UGAAzA$B1CAA"+?e?n?u?~k?X@u#o3kX@A1lAȼHBgBf?d?msD?us?~^a?~6@u%Ӟ&LX?O"@FAYB*A$9?ect?mN?u?~e{?	@u%Wrg#kX?.#?cv!??)?,~?ev?m?u?~f?X@u$ bU4{^X>:?RB?>A:AK?e?mf?u@?~a?@u$)b.6̶X>CJ>$?gsFԾ?d?mo2?uqz?~^?~@u&aƐ`OF	X@ ABwnB1A?e?n?u>?~kY?2@u$JcбXX>_>?\@@(A?e
 A?m?u|?~_?@u$L[N^??pvp>B*=?eb?m?u}d?~_?#@u#h$ՐY]:I X@]xA#A4BA?dR?mb?uja?~\?~K@u#^hlKXC4DyD")C#KBvo?d?m_?uh?~\?~!?	+go@u#(ut**X??]%?@o>?!#@Y0?-?e?m~?uzr?~_l?~@u#gX>>%(?aq?e?m|I?uy?~_3?~@u&:&mQ%X?H@9AWAeA4?e?mb?u~?~j?	@u$k:G)~01X?O?K@]A=A|*?dة?m`]?uh?~\?~'@u&h*gsX@YAArBAZ?eL?ni?u?~kO?+@u$;qeVX>>?Jɘ?6@1?d?m[?uf.?~\;?~@u$M}>	X>>Ԡ?:)`@@x?dژ?ma?ui?~\?~:@u$)u|/NX?X@vAB/B%'Ae?d?mk?uol?~]?~@u&1G4پuX?#A?E??Ҋ@@'h?e?mP?ut?~j?@u#ڻ.5||bX?X?WP?amg,'?dZ?mh?um?~]g?~@u#pb.m,X>G>~?O:=4?d?mlc?uo?~]?~@u%)3X>˸;>	?	V9A 
?e ?m?u8?~iu?@u%^[*`X@yAaB0B̦Bo?e ?m?u?~a??@u&}.>OX>ɝq>?&'@-Gav?e?n?u?~k?Y@u%	A]bn^X@fA%jAAAe?e
?m?u?~i?@u#gj	` X?@A6B#
Am?d߯?meF?uk?~]?~m	@u&kXb)XBBvCFCB?e?m?u?~hk?d@u&Ii2MDX?O?W?lj	?h.?ea?m?uN?~i(?@u#4'X@/A/A[B4B }?d7?mas?ui?~\?~6@u&!ds'oX>1?P?qHn.Y?e?m?ux?~k?@u&BU]1~X?M@"UA$lB
 AH?e&?m?u±?~j?	@u#Qf5|`X>k?c?ѻ?+?an?d"?m&?uG?~Wg?~@u&ϬOca2XANBB1CB2BDM?en?m[?u?~f6?@u&Jn.ޓ9KX>]??0awl?50?e?n?ul?~k?O@u&?~4X>o?	?͉Bc$N?e
-?mۋ?u?~g?
@u$Gg6X>>]?77
߇>m?d?mq=?ur?~^2?~@u#EF(X?@{AAB8DA ?d?m@v?uVs?~Y?~]	@u&fvr)nTX?m?&@@@#2A9F?eL?n?u3?~k?{@u&˚L+XBeOCC1CCBN?eh?m?u?~e?@u%x՞2X>iJ??V@ܪ}@?e3?m?u?~hR?U@u&]l(3%gȦX>Ĉ>r?^/?@?e?n?u?~k?x@u#HX>ç	>]?^?"e?df?mOi?u_?~[?~3@u$ʗ-X?@???i2@.j@X[?dj?mP?u_?~[>?~G@u%fY;T5pX>>>?Z"'?q?e>?m'?u?~c?@u$>/U#bUoX>?]?nB?8
+?mۋ?u?~g?
@u$Gg6X>>]?77
߇>m?d?mq=?ur?~^2?~@u#EF(X?@{AAB8DA ?d?m@v?uVs?~Y?~]	@u&˚L+XBeOCC1CCBN?eh?m?u?~e?@u%x՞2X>iJ??V@ܪ}@?e3?m?u?~hR?U@u&]l(3%gȦX>Ĉ>r?^/?@?e?n?u?~k?x@u#HX>ç	>]?^?"e?df?mOi?u_?~[?~3@u$ʗ-X?@???i2@.j@X[?dj?mP?u_?~[>?~G@u%fY;T5pX>>>?Z"'?q?e>?m'?u?~c?@u$>/U#bUoX>?]?nB?8
 ??dd?mC?uX8?~Z?~@u$@"S@ֵX>Υ>Q.>>M@H?d%?mEs?uYV?~Z5?~@u$Aa!X>=?&?V?h@zW?d?mFi?uY?~ZK?~@u%Q\jmXBXBwCg57C30B&>?eF?m?uh?~h?-@u#-*&J+X@A7A%ANAVXC?dl?m?u=_?~U?~@u#-(WJVX@6AEA AgoA^?dU?m?u4?~TT?~@u$zʇI͸X>?C?)@L2@?d?ms?ut+?~^o?~>@u&'\#}TX?8o@GAdcB0A?eĻ?n4?u+?~k?b	@u#ͶgX??6?@@n?d?m)$?uH?~W?~@u#,/%7X>B?6?k[ɿQ?2e?d[?4?cD?]?d|?m7?uP?~X?~@u$
2(C$DX>;:???S??do?mB(?uWn?~Y?~u@u&x+
 7X>>?8\?
 H{?e
?n?u͉?~l?@u&*>-XA";BArB?BnB[Q?e?m;?u?~gy?@u&ZޯȻjdX>z)>U?d+R>|@1f?eΦ?n ?u*?~lH?@u#\]X?eX>ki>Aha?dD?l!?u-.?~S@?~\@u#r󸸀m?54X?7?@Cs@_+}@M?dCc?l?u,?~S&?~L@u&mtLHcX?}?V@~7<4ļ?eo?m?u?~fG?@u&`CjrX?
 ?.?¿]d?1Q?eT?n?uͦ?~l?@u&ŭqÒX?>|?D-?e?m?u?~hn?f@u&x;9[pX?{B?e@Dv@ļ@k?e?n?u?~l=?@u%F&qnX>b@EAAAroAl]h?e?m?uc?~h?w @u%hq[ÛX?J8?w@#.@N@7?eh?m*?u+?~e?@u&7hX@NA9AAA?e?n?u?~m-?R@u%zpRX>>?ʗ?zV@?eh?m
 ?u?~i?3@u&[LX>?)/?i@?@5?ei?mŏ?uf?~e?@u$Y>:?B%*X>*?=?/?&>]?d|?m A?uC?~V?~@u%N
ѣk+pEX>r@':A^BFAgy?eV?m?u?~d?	@u/OnUOtxXAwB:O)BdA)A`?[?f?q^?}q?~@u/In|g)bX?XIN?@8J-Y@߳E?[?f?q?}~?~@u,mdp:TX>4-@PA9+B6A?[q?g$?qȲ?}ɏ?~	@u,ajDX@AMNzA{?A{A%?[?gK?q?}?~/@u,?SO"XB'@B௶C`=CM!B?[?g?q?f?o@ii
 ?[?g?q?}?~w@u,bX@NAHȵAA;A?[?gv?{n@ѿ/AdžAd?[?g?q?}ȫ?~ @u.93|XA0A3B,9B)1dAb?[?f?q&?}?~U@u/c~׼5X@AAPAIJ"@?[U?f?q?}?~
-@u/*,ׯ
KX??S2?|@#i??[z?fU?q??}þ?~@u./hי4X>r?y?_~?[?f?qh?}?~-@u-m·X>>0?ݿj	i?[&?g?q?}?~@u,+kחM]YX>,>	K?pS
	?[?g?qa?}?~@u,x"צKX>=x>V`?@?%?[?g?q4?}?~@u/Yc+׆zX>F>?,(@I@k0?[_?fx?q?}l?~@u.^"ףWwX>??@"z?e?[?f?q?}\?~X@u/ׂ&۟JX@6AixBB	B]*?[?f?q?}?~\@u-TyH?]X>@MAiXBnAp?[]?g
2?q?}_?~4	@u-$2ipX0X@thAB?tAAUAX}?[F?g?q:?}?~@u-޿WzTX@,KAF2AEBVbBhg?[?g?q?}Ʈ?~@u-ĿVTiX>??/@NԐ?[a?g ?q?}ƶ?~@u.4%o_s 4X>˖?5&?S@
'pVS?[?fA?q~?}ų?~+@u,#SlvX>?7??E??[h?g?qx?}ȑ?~@u-K$IXAuBh/wBBO.B?[|?g?q?}ǥ?~_@u,#,ڿ1wvX??D'?	H@?s?[?g?qW?}?~@u.gӉ6$@X@03X@?6@T5@>@~?[S?f?q?}A?~@u,UP,nX??Q@VM@ī@Y?[?g?q?}?~@u/V#Ч#XA2AAB2A]Aq/?[?fh?q?}¹?~V@u-9E#f]k#X?Ui@rAJ~BF
+@u/*,ׯ
KX??S2?|@#i??[z?fU?q??}þ?~@u./hי4X>r?y?_~?[?f?qh?}?~-@u-m·X>>0?ݿj	i?[&?g?q?}?~@u,+kחM]YX>,>	K?pS
	?[?g?qa?}?~@u,x"צKX>=x>V`?@?%?[?g?q4?}?~@u/Yc+׆zX>F>?,(@I@k0?[_?fx?q?}l?~@u.^"ףWwX>??@"z?e?[?f?q?}\?~X@u/ׂ&۟JX@6AixBB	B]*?[?f?q?}?~\@u-TyH?]X>@MAiXBnAp?[]?g
2?q?}_?~4	@u-$2ipX0X@thAB?tAAUAX}?[F?g?q:?}?~@u-޿WzTX@,KAF2AEBVbBhg?[?g?q?}Ʈ?~@u-ĿVTiX>??/@NԐ?[a?g ?q?}ƶ?~@u.4%o_s 4X>˖?5&?S@
'pVS?[?fA?q~?}ų?~+@u,#SlvX>?7??E??[h?g?qx?}ȑ?~@u-K$IXAuBh/wBBO.B?[|?g?q?}ǥ?~_@u,#,ڿ1wvX??D'?	H@?s?[?g?qW?}?~@u.gӉ6$@X@03X@?6@T5@>@~?[S?f?q?}A?~@u,UP,nX??Q@VM@ī@Y?[?g?q?}?~@u/V#Ч#XA2AAB2A]Aq/?[?fh?q?}¹?~V@u-9E#f]k#X?Ui@rAJ~BF
 A'?[ ?g?q?}9?~	@u-I`LJwT3X>?#?S(@Mb.^;?[!?g	}?q?}?~@u,6y` X>]?A!BA?[?gb?q4?}?~	@u/ָB 3IX>Z??Y>@7C?[?f?q.?}L?~@u/|jpX?+?D?AWA@`?[?fE?q?}?~@u-i]	X>?@$?-}??s?릥?[?fO?q?}?~@u,BgX? ?_I???[!?g?q?}ǰ?~e@u-ǀUZ"X?C@_4AVBAF?[#?f?q?}?~k	@u/H&0\X?@?R?NA"A|?[?f8?q?}?~@u.@Tq\ӯX>U>g?%<>0V?[?f?q??}?~@u/~X?V?~>?@<@7]2?[?f?q?}’?~>@u/mY!1ֿS_@2?c?w>!(?[ ?f?q?}?~@u.j鶿X?
 5A?x?kU="FB?[Z?f?q+?}ć?~s@u.H0ٿ֚UX>@7ƻA	@u/f
 Tg֢QQX?v&.@ɕAzeB3	A?[(?f&?q6?}B?~	@u-šTtX>>o?9
H[A1?[?fw?q?}p?~@u-diÿ֤cbX>N*?3]o?@7H??[R?g?q?}Ɣ?~@u.'\֍hX@ٿA6hA0B!A?[,?f~?q?}R?~R@u-Yop'!V(X>?w??~>7?[?g?q?}m?~@u.;֖.hX>?^@A~@?[?fƁ?a)AATA9?[Ǐ?f?q4?}!?~p @u-31{W5yX?f@U/AtBPnB?[?f?q?}?~^	@u-2Dd8:X?
@@ -144,4 +144,4 @@ Z
 ?@
 ?>@?[@?f?q{a?}8?~@u-Bm=ڮXA@SAUVA|A|B6d?[:?f3?qx?}?~e),@u/u'DEX@AhB
8B'B?["v?f?qn?}3?~@u,s2/X>
 =E8>@5@BNA'xB?[?fC?qi?}L?~%@u,Mp,ꦦXCZM%C'DGD<?[ ?f>?qm?}?~
\ No newline at end of file

From 02e6e0280b427e967f086c26e0f17b3355f63242 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 15:45:57 -0700
Subject: [PATCH 24/31] update changes documentation

---
 doc/changes.rst     | 14 ++++++++++++++
 py/desitarget/QA.py |  2 +-
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/doc/changes.rst b/doc/changes.rst
index 2b30a3055..3dad9921d 100644
--- a/doc/changes.rst
+++ b/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 by 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-16.
+    * 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)
 -------------------
diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 52800be8c..0928477a3 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -2030,7 +2030,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
         headerlist.sort()
         # ADM edit SV target categories to their initial letters to squeeze space.
         hl = headerlist.copy()
-        if svs == 'SV':
+        if svs[0:2] == 'SV':
             for tc in 'QSO', 'ELG', 'LRG', 'BGS', 'MWS':
                 hl = [h.replace(tc, tc[0]) for h in hl]
         # ADM truncate the bit names at "trunc" characters to pack them more easily.

From e3478d859f8ab25d8e95fee98e1b3ad414fca704 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 16:27:09 -0700
Subject: [PATCH 25/31] update default number of processors for making QA plots

---
 bin/run_target_qa | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/bin/run_target_qa b/bin/run_target_qa
index 7dcf2eaa0..83fd6cb86 100755
--- a/bin/run_target_qa
+++ b/bin/run_target_qa
@@ -13,7 +13,8 @@ from desiutil.log import get_logger
 log = get_logger()
 
 import multiprocessing
-nproc = multiprocessing.cpu_count() // 2
+# ADM the rather large denominator guards against memory issues.
+nproc = multiprocessing.cpu_count() // 8
 
 from argparse import ArgumentParser
 ap = ArgumentParser()

From 5190c8962929fb2dd8311854c029eb871ac7fc37 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 16:28:19 -0700
Subject: [PATCH 26/31] further update default number of processors for making
 QA plots

---
 py/desitarget/QA.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/py/desitarget/QA.py b/py/desitarget/QA.py
index 0928477a3..281733943 100644
--- a/py/desitarget/QA.py
+++ b/py/desitarget/QA.py
@@ -1518,7 +1518,7 @@ def _in_desi_footprint(targs):
 
 def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True,
                   imaging_map_file=None, truths=None, objtruths=None, tcnames=None,
-                  cmx=False, bit_mask=None, mocks=False, numproc=32):
+                  cmx=False, bit_mask=None, mocks=False, numproc=8):
     """Make DESI targeting QA plots given a passed set of targets.
 
     Parameters
@@ -1556,7 +1556,7 @@ def make_qa_plots(targs, qadir='.', targdens=None, max_bin_area=1.0, weight=True
         constraints) instead of the main survey bits.
     mocks : :class:`boolean`, optional, default=False
         If ``True``, add plots that are only relevant to mocks at the bottom of the webpage.
-    numproc : :class:`int`, optional, defaults to 32
+    numproc : :class:`int`, optional, defaults to 8
         The number of parallel processes to use to generate plots.
 
     Returns
@@ -1746,7 +1746,7 @@ def _generate_plots(objtype):
 
 def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.',
                  clip2foot=False, weight=True, imaging_map_file=None,
-                 tcnames=None, systematics=True, numproc=32):
+                 tcnames=None, systematics=True, numproc=8):
     """Create a directory containing a webpage structure in which to embed QA plots.
 
     Parameters
@@ -1781,7 +1781,7 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
         for those specific bits. A useful speed-up when testing
     systematics : :class:`boolean`, optional, defaults to ``True``
         If sent, then add plots of systematics to the front page.
-    numproc : :class:`int`, optional, defaults to 32
+    numproc : :class:`int`, optional, defaults to 8
         The number of parallel processes to use to generate plots.
 
     Returns

From ad47ff7ffd8759059a5f4cb435ac318b933c1e28 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Mon, 8 Apr 2019 16:40:23 -0700
Subject: [PATCH 27/31] one clarification in changes file [ci skip]

---
 doc/changes.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changes.rst b/doc/changes.rst
index 3dad9921d..dc1a90aad 100644
--- a/doc/changes.rst
+++ b/doc/changes.rst
@@ -9,7 +9,7 @@ desitarget Change Log
     * Resolve sky locations and SV targets by 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-16.
+    * 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.

From b5003bb563853e4269796d61ff51ffb2093c93fe Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Sun, 21 Apr 2019 00:51:14 -0700
Subject: [PATCH 28/31] check that unit tests still pass after hiatus

---
 doc/changes.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changes.rst b/doc/changes.rst
index dc1a90aad..0a0ff17db 100644
--- a/doc/changes.rst
+++ b/doc/changes.rst
@@ -6,7 +6,7 @@ desitarget Change Log
 -------------------
 
 * Further updates and enhancements for DR8 [`PR #490`_]. Includes:
-    * Resolve sky locations and SV targets by North/South regions.
+    * Resolve sky locations and SV targets into 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.

From ca704728395289d6c9b80dbd9890338be51dd323 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Fri, 26 Apr 2019 08:33:25 -0700
Subject: [PATCH 29/31] minor update to re-trigger Travis

---
 doc/changes.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changes.rst b/doc/changes.rst
index 0a0ff17db..dc1a90aad 100644
--- a/doc/changes.rst
+++ b/doc/changes.rst
@@ -6,7 +6,7 @@ desitarget Change Log
 -------------------
 
 * Further updates and enhancements for DR8 [`PR #490`_]. Includes:
-    * Resolve sky locations and SV targets into North/South regions.
+    * Resolve sky locations and SV targets by 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.

From 437a98a56fed3c32223cc6747e2e294b727d6c22 Mon Sep 17 00:00:00 2001
From: Adam Myers 
Date: Fri, 26 Apr 2019 15:18:05 -0700
Subject: [PATCH 30/31] attempt to be simultaneously compatible with fitsio
 v0.9 and v1.0 to pass both NERSC and Travis unit tests

---
 py/desitarget/io.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/py/desitarget/io.py b/py/desitarget/io.py
index 36623eb91..5e819b2b0 100644
--- a/py/desitarget/io.py
+++ b/py/desitarget/io.py
@@ -277,7 +277,13 @@ def add_photsys(indata):
     # ADM only add the PHOTSYS column if RELEASE exists.
     if 'RELEASE' in indata.dtype.names:
         # ADM add PHOTSYS to the data model.
-        pdt = [('PHOTSYS', '|S1')]
+        # 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'):
+            pdt = [('PHOTSYS', '
Date: Tue, 30 Apr 2019 08:11:05 -0700
Subject: [PATCH 31/31] one last test of Travis before merging into master

---
 doc/changes.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changes.rst b/doc/changes.rst
index dc1a90aad..fc640056f 100644
--- a/doc/changes.rst
+++ b/doc/changes.rst
@@ -6,7 +6,7 @@ desitarget Change Log
 -------------------
 
 * Further updates and enhancements for DR8 [`PR #490`_]. Includes:
-    * Resolve sky locations and SV targets by North/South regions.
+    * 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.