Merge branch 'master' into fiber-crosstalk

bin/copyprod

@@ -1,39 +1,153 @@
 #!/usr/bin/env python
-# See top-level LICENSE.rst file for Copyright information
 
 """
-Copy a DESI production from one directory to another
-
-Currently this only copies the directory tree and the extracted frame files.
-This script will be expanded to support copying through any level of
-processing, including the metadata of how we got to that processing step.
+Create a nested directory tree with file symlinks to "copy" one prod to another
 
 Stephen Bailey
-April 2015
+March 2021
 """
-from __future__ import absolute_import, division, print_function
 
-import os
+import os, sys, glob
 import shutil
+from collections import Counter
 import optparse
 
+import numpy as np
+from astropy.table import Table
+
+from desiutil.log import get_logger
+from desispec.workflow.tableio import load_table, write_table
+
 parser = optparse.OptionParser(usage = "%prog [options] indir outdir")
-# parser.add_option("-i", "--input", type="string",  help="input data")
-# parser.add_option("-x", "--xxx",   help="some flag", action="store_true")
-opts, args = parser.parse_args()
+parser.add_option("--explist", help="file with NIGHT EXPID to copy")
+parser.add_option("--fullcopy", action="store_true",
+        help="Copy files instead of linking")
+parser.add_option("--abspath", action="store_true",
+        help="Link to absolute path instead of relative path")
+parser.add_option("--filter-exptables", action="store_true",
+        help="Filter exposure_tables to ignore entries not in input explist")
 
+opts, args = parser.parse_args()
 inroot, outroot = args
+log = get_logger()
+
+#- Get list of NIGHT EXPID to copy
+if opts.explist is not None:
+    explist = Table.read(opts.explist, format='ascii')
+else:
+    rows = list()
+    for dirname in sorted(glob.glob(f'{inroot}/exposures/20??????/????????')):
+        nightdir, expid = os.path.split(dirname)
+        try:
+            night = int(os.path.basename(nightdir))
+            expid = int(expid)
+            rows.append((night, expid))
+        except ValueError:
+            pass
+
+    explist = Table(rows=rows, names=('NIGHT', 'EXPID'))
+
+num_nights = len(np.unique(explist['NIGHT']))
+num_expids = len(np.unique(explist['EXPID']))
+assert num_expids == len(explist)
+log.info(f'Linking {num_expids} exposures on {num_nights} nights')
+
+#- What type of copy or link are we making?
+if opts.fullcopy:
+    link = shutil.copy2
+    log.debug('Performing full copy of selected NIGHT/EXPID')
+else:
+    link = os.symlink
+    log.debug('Creating symlinks')
+
+inroot = os.path.abspath(inroot)
+outroot = os.path.abspath(outroot)
+log.info(f'Creating links in {outroot} to files in {inroot}')
 
 if not os.path.exists(outroot):
     os.makedirs(outroot)
 
-#- Copy exposures
-for indir, subdirs, filenames in os.walk(inroot+'/exposures'):
-    outdir = indir.replace(inroot, outroot)
-    print(outdir)
-    if not os.path.exists(outdir):
-        os.makedirs(outdir)
+def link_dirfiles(indir, outdir, abspath=False):
+    """
+    Create relative links in outdir to all files in indir
+    """
+    os.makedirs(outdir, exist_ok=True)
+    if abspath:
+        srcpath = indir
+    else:
+        srcpath = os.path.relpath(indir, outdir)
+
+    for infile in sorted(glob.glob(f'{indir}/*.*')):
+        basename = os.path.basename(infile)
+        src = os.path.join(srcpath, basename)
+        dst = os.path.join(outdir, basename)
+        link(src, dst)
+
+#- calibnight: link all files for requested nights
+for night in np.unique(explist['NIGHT']):
+    log.info(f'calibnight/{night}')
+    indir = f'{inroot}/calibnight/{night}'
+    outdir = f'{outroot}/calibnight/{night}'
+    link_dirfiles(indir, outdir, opts.abspath)
+
+#- preproc, exposures: link all files for requested night/expid
+for subdir in ['preproc', 'exposures']:
+    for night, expid in explist['NIGHT', 'EXPID']:
+        log.info(f'{subdir}/{night}/{expid:08d}')
+        indir = f'{inroot}/{subdir}/{night}/{expid:08d}'
+        outdir = f'{outroot}/{subdir}/{night}/{expid:08d}'
+        link_dirfiles(indir, outdir, opts.abspath)
+
+#- exposure tables: trim to requested exposures
+for night in np.unique(explist['NIGHT']):
+    yearmm = night//100
+    infile = f'{inroot}/exposure_tables/{yearmm}/exposure_table_{night}.csv'
+    outfile = f'{outroot}/exposure_tables/{yearmm}/exposure_table_{night}.csv'
+    os.makedirs(os.path.dirname(outfile), exist_ok=True)
+    if opts.filter_exptables:
+        exptab = load_table(infile, tabletype='exptable')
+        ignore = np.in1d(exptab['EXPID'], explist['EXPID'], invert=True)
+        if np.any(ignore):
+            msg = f'Night {night} ignoring '
+            msg += ', '.join([f'{n} {obstype}' for obstype, n in Counter(exptab['OBSTYPE'][ignore]).items()])
+            msg += ' exposures not included in explist'
+            log.warning(msg)
+
+        exptab['LASTSTEP'][ignore] = 'ignore'
+        write_table(exptab, outfile, tabletype='exptable')
+    else:
+        log.info(f'Night {night} copying exposure tables unchanged')
+        shutil.copy2(infile, outfile)
+
+
+#- processing tables: trim to requested exposures.
+#- Disabled for now while we figure out exactly what we want.
+# for night in np.unique(explist['NIGHT']):
+#     night_expids = explist['EXPID'][explist['NIGHT'] == night]
+#     tmp = glob.glob(f'{inroot}/processing_tables/processing_table_*-{night}.csv')
+#     if len(tmp) == 1:
+#         infile = tmp[0]
+#     elif len(tmp) == 0:
+#         log.error(f'Unable to find processing table for {night}')
+#         continue
+#     elif len(tmp) == 0:
+#         log.error(f'Multiple processing tables for {night} ?!?')
+#         continue
+# 
+#     outfile = '{}/processing_tables/{}'.format(
+#             outroot, os.path.basename(infile))
+#     proctab = load_table(infile, tabletype='proctable')
+# 
+#     #- Convert string "expid1|expid2" entries into keep or not
+#     keep = np.zeros(len(proctab), dtype=bool)
+#     for i, expids in enumerate(proctab['EXPID']):
+#         if np.any(np.in1d(expids, night_expids)):
+#             keep[i] = True
+# 
+#     proctab = proctab[keep]
+# 
+#     os.makedirs(os.path.dirname(outfile), exist_ok=True)
+#     write_table(proctab, outfile, tabletype='proctable')
+
+log.info(f'Production copied into {outroot}')
 
-    for name in filenames:
-        if name.startswith('frame-'):
-            shutil.copy2(indir+'/'+name, outdir+'/'+name)

bin/desi_average_flux_calibration

@@ -68,6 +68,24 @@ if __name__ == '__main__':
             seeing.append(float(header["seeing"]))
 
         cal=read_flux_calibration(filename)
+        if np.any(np.isnan(cal.calib)) :
+            log.error("calib has nan")
+            continue
+
+        mcalib = np.median(cal.calib,axis=0)
+
+        # undo the heliocentric/barycentric correction
+        #- check for heliocentric correction
+        if 'HELIOCOR' in header.keys() :
+            heliocor=header['HELIOCOR']
+            log.info("Undo the heliocentric correction scale factor {} of the calib".format(heliocor))
+            # wavelength are in solar system frame
+            # first divide the multiplicative factor to have wavelength in KPNO frame
+            wave_in_kpno_system = cal.wave/heliocor
+
+            # now we want the wave grid cal.wave to be in the KPNO frame,
+            # we have to resample the calib vector
+            mcalib = np.interp(cal.wave,wave_in_kpno_system,mcalib)
 
         if wave is None :
             wave=cal.wave
@@ -77,10 +95,6 @@ if __name__ == '__main__':
         if exptime <= 2. : # arbitrary cutoff
             print("skip exptime=",exptime)
             continue
-        if np.any(np.isnan(cal.calib)) :
-            print("ERROR calib has nan")
-            continue
-        mcalib = np.median(cal.calib,axis=0)
 
         calibs.append(mcalib/exptime)
 

bin/desi_merge_zbest_tiles

@@ -0,0 +1,143 @@
+#!/usr/bin/env python
+
+"""
+Make summary redshift catalogs from tile-based redshifts, e.g. blanc
+"""
+
+import os, sys, glob, collections
+import argparse
+import numpy as np
+from numpy.lib.recfunctions import rec_append_fields, join_by
+import fitsio
+from astropy.table import Table, vstack
+from desiutil.log import get_logger
+from desiutil import depend
+
+parser = argparse.ArgumentParser(usage = "{prog} [options]")
+parser.add_argument("--group", type=str, required=True,
+        help="exposure grouping, e.g. 'all' or 'deep'")
+parser.add_argument("-o", "--output", type=str, required=True,
+        help="output catalog filename")
+parser.add_argument("--reduxdir", type=str,
+        help="overrides $DESI_SPECTRO_REDUX/$SPECPROD")
+parser.add_argument("--first-tile", type=int,
+        help="First TILEID to include")
+parser.add_argument("--last-tile", type=int,
+        help="Last TILEID to include (inclusive, not python indexing style)")
+parser.add_argument("--tiles", type=int, nargs="*",
+        help="TILEIDs to include (space separated)")
+parser.add_argument("--no-fibermap", action="store_true", help="Do not merge with fibermap")
+parser.add_argument("--no-scores", action="store_true", help="Do not merge with coadd scores")
+# parser.add_argument("-v", "--verbose", action="store_true", help="some flag")
+
+args = parser.parse_args()
+log = get_logger()
+
+if args.reduxdir is None:
+    args.reduxdir = os.path.expandvars('$DESI_SPECTRO_REDUX/$SPECPROD')
+
+tiledir = f'{args.reduxdir}/tiles'
+log.info(f'Looking for per-tile zbest files in {tiledir}/TILEID/{args.group}/')
+assert os.path.isdir(tiledir)
+
+zbestfiles = list()
+for filename in sorted(glob.glob(f'{tiledir}/*/{args.group}/zbest-*.fits')):
+    tileid = int(os.path.basename(os.path.dirname(os.path.dirname(filename))))
+    if args.first_tile is not None and tileid < args.first_tile:
+        continue
+    if args.last_tile is not None and tileid > args.last_tile:
+        continue
+    if args.tiles is not None and tileid not in args.tiles:
+        continue
+
+    zbestfiles.append(filename)
+
+nzbfiles = len(zbestfiles)
+zbx = list()
+
+for i, filename in enumerate(zbestfiles):
+    print(f'{i+1}/{nzbfiles}: {filename}')
+    zb = fitsio.read(filename, 'ZBEST')
+
+    tmp = os.path.basename(filename).replace('zbest-', 'coadd-', 1)
+    coaddfile = os.path.join(os.path.dirname(filename), tmp)
+    if not args.no_fibermap:
+        fm = fitsio.read(coaddfile, 'FIBERMAP')
+        #- Sorted the same
+        assert np.all(zb['TARGETID'] == fm['TARGETID'])
+        #- TARGETID is only column in common
+        assert (set(zb.dtype.names) & set(fm.dtype.names)) == set(['TARGETID',])
+        zb = join_by('TARGETID', zb, fm)
+    else:
+        tileid = np.ones(len(zb), dtype=np.int16)*tileid
+        zb = rec_append_fields(zb, 'TILEID', tileid)
+
+    if not args.no_scores:
+        scores = fitsio.read(coaddfile, 'SCORES')
+        #- Sorted the same
+        assert np.all(zb['TARGETID'] == scores['TARGETID'])
+        #- TARGETID is only column in common
+        assert (set(zb.dtype.names) & set(scores.dtype.names)) == set(['TARGETID',])
+        zb = join_by('TARGETID', zb, scores)
+
+    #- Handle a few dtype special cases
+    zb = Table(zb)
+    if 'NUMOBS_MORE' in zb.colnames and zb['NUMOBS_MORE'].dtype != np.dtype('>i4'):
+        zb['NUMOBS_MORE'] = zb['NUMOBS_MORE'].astype('>i4')
+    if 'RELEASE' in zb.colnames and zb['RELEASE'].dtype != np.dtype('>i2'):
+        zb['RELEASE'] = zb['RELEASE'].astype('>i2')
+
+    zbx.append(zb)
+
+#- Determine union and intersection of columns present in the files
+#- since the fiberassign datamodel evolved for exactly which targeting
+#- columns were provided
+all_columns = set()
+joint_columns = None
+for zb in zbx:
+    all_columns.update(zb.colnames)
+    if joint_columns is None:
+        joint_columns = set(zb.colnames)
+    else:
+        joint_columns &= set(zb.colnames)
+
+#- Add *_TARGET columns if needed, and drop other columns that aren't in common
+add_columns = set()
+drop_columns = set()
+for colname in (all_columns - joint_columns):
+    if colname.endswith('_TARGET'):
+        add_columns.add(colname)
+    else:
+        drop_columns.add(colname)
+
+#- update individual tables to have the same columns
+for zb in zbx:
+    for colname in add_columns:
+        if colname not in zb.colnames:
+            zb[colname] = np.zeros(len(zb), dtype=int)
+    for colname in drop_columns:
+        if colname in list(zb.colnames):
+            zb.remove_column(colname)
+
+#- Standardize column order to match last tile
+columns = zbx[-1].colnames
+for i in range(len(zbx)):
+    zb = zbx[i]
+    if zb.columns != columns:
+        zbx[i] = zb[columns]  #- reorders columns
+
+#- Finally! make the stacks redshift catalog
+zcat = vstack(zbx)
+
+#- Add record of inputs
+zcat.meta['TILEDIR'] = os.path.normpath(tiledir)
+for i, filename in enumerate(zbestfiles):
+    key = f'IN{i:06d}'
+    zcat.meta[key] = filename.replace(tiledir, 'TILEDIR')
+
+depend.add_dependencies(zcat.meta)
+zcat.meta['SPECPROD'] = os.path.basename(os.path.normpath(args.reduxdir))
+zcat.meta['EXTNAME'] = 'ZCATALOG'
+zcat.write(args.output)
+
+

bin/desi_run_night

@@ -21,6 +21,8 @@ def parse_args():  # options=None):
     #                         "Overwrites the environment variable SPECPROD")
     parser.add_argument("-q", "--queue", type=str, required=False, default='realtime',
                         help="The queue to submit jobs to. Default is realtime.")
+    parser.add_argument("-r", "--reservation", type=str, required=False, default=None,
+                        help="The reservation to submit jobs to. If None, it is not submitted to a reservation.")
     parser.add_argument("--exp-table-path", type=str, required=False, default=None,
                         help="Directory name where the output exposure table should be saved.")
     parser.add_argument("--proc-table-path", type=str, required=False, default=None,
@@ -33,6 +35,8 @@ def parse_args():  # options=None):
     parser.add_argument("--error-if-not-available", action="store_true",
                         help="Raise an error instead of reporting and moving on if an exposure "+\
                              "table doesn't exist.")
+    parser.add_argument("--overwrite-existing", action="store_true",
+                        help="Give this flag if you want to submit jobs even if scripts already exist.")
 
     args = parser.parse_args()