Merge pull request #424 from desihub/qa_frame_script

QA frame script

bin/desi_qa_frame

@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+#
+# See top-level LICENSE.rst file for Copyright information
+#
+# -*- coding: utf-8 -*-
+
+"""
+This script generates QA related to a production
+"""
+
+import desispec.scripts.qa_frame as qa_frame
+
+
+if __name__ == '__main__':
+    args = qa_frame.parse()
+    qa_frame.main(args)

doc/changes.rst

@@ -8,6 +8,11 @@ desispec Change Log
 * Small fixes to desi_qa_prod and qa_prod
 * Fixes integration tests for desisim newexp refactor
 * Removes spectra grouping by brick; nside=64 healpix grouping default
+* Add get_nights method to io.meta
+* Add search_for_framefile method to io.frame
+* Add desi_qa_frame script to generate frame QA
+
+.. _`#422`: https://github.com/desihub/desispec/pull/422
 
 0.15.2 (2017-07-12)
 -------------------

doc/index.rst

@@ -13,6 +13,7 @@ Contents
     pipeline.rst
     dev.rst
     coadd.rst
+    qa.rst
     changes.rst
     api.rst
 

doc/qa.rst

@@ -0,0 +1,108 @@
+.. _qa:
+
+*****************
+Quality Assurance
+*****************
+
+Overview
+========
+
+The DESI spectroscopic pipeline includes a series of
+routines that monitor the quality of the pipeline products
+and may be used to inspect outputs across exposures, nights,
+or a full production.
+
+
+Scripts
+=======
+
+desi_qa_frame
++++++++++++++
+
+Generate the QA for an input frame file.
+The code can be written anywhere and the
+output is written to its "proper" location.
+
+usage
+-----
+
+Here is the usage::
+
+    usage: desi_qa_frame [-h] --frame_file FRAME_FILE [--reduxdir PATH]
+                         [--make_plots]
+
+    Generate Frame Level QA [v1.0]
+
+    optional arguments:
+      -h, --help            show this help message and exit
+      --frame_file FRAME_FILE
+                            Frame filename. Full path is not required nor desired.
+      --reduxdir PATH       Override default path ($DESI_SPECTRO_REDUX/$SPECPROD)
+                            to processed data.
+      --make_plots          Generate QA figs too?
+
+
+examples
+--------
+
+Generate the QA YAML file::
+
+    desi_qa_frame --frame_file=frame-r7-00000077.fits
+
+Generate the QA YAML file and figures::
+
+    desi_qa_frame --frame_file=frame-r7-00000077.fits --make_plots
+
+desi_qa_prod
+++++++++++++
+
+This script is used to both generate and analyze the
+QA outputs for a complete production.
+
+usage
+-----
+
+Here is the usage::
+
+    usage: desi_qa_prod [-h] --specprod_dir SPECPROD_DIR
+                        [--make_frameqa MAKE_FRAMEQA] [--slurp] [--remove]
+                        [--clobber] [--channel_hist CHANNEL_HIST]
+
+    Generate/Analyze Production Level QA [v1.2]
+
+    optional arguments:
+      -h, --help            show this help message and exit
+      --specprod_dir SPECPROD_DIR
+                            Path containing the exposures/directory to use
+      --make_frameqa MAKE_FRAMEQA
+                            Bitwise flag to control remaking the QA files (1) and
+                            figures (2) for each frame in the production
+      --slurp               slurp production QA files into one?
+      --remove              remove frame QA files?
+      --clobber             clobber existing QA files?
+      --channel_hist CHANNEL_HIST
+                            Generate channel histogram(s)
+
+
+
+frameqa
+-------
+
+One generates the frame QA, the YAML and/or figure files
+with the --make_frameqa flag::
+
+    desi_qa_prod --make_frameqa=1  # Generate all the QA YAML files
+    desi_qa_prod --make_frameqa=2  # Generate all the QA figure files
+    desi_qa_prod --make_frameqa=3  # Generate YAML and figures
+
+The optional --remove and --clobber flags can be used to remove/clobber
+the QA files.
+
+slurp
+-----
+
+By using the --slurp flag, one generates a full
+YAML file of all the QA outputs::
+
+    desi_qa_prod --slurp   # Collate all the QA YAML files
+    desi_qa_prod --slurp --remove  # Collate and remove the individual files

py/desispec/io/__init__.py

@@ -23,7 +23,7 @@
 from .meta import (findfile, get_exposures, get_files, get_raw_files,
                    rawdata_root, specprod_root, validate_night,
                    get_pipe_plandir, get_pipe_rundir, get_pipe_scriptdir,
-                   get_pipe_logdir, get_pipe_faildir)
+                   get_pipe_logdir, get_pipe_faildir, get_nights)
 from .params import read_params
 from .qa import (read_qa_frame, read_qa_data, write_qa_frame, write_qa_brick,
                  load_qa_frame, write_qa_exposure, write_qa_prod)

py/desispec/io/frame.py

@@ -15,11 +15,10 @@
 from desiutil.io import encode_table
 
 from ..frame import Frame
-from .meta import findfile
+from .meta import findfile, get_nights, get_exposures
 from .util import fitsheader, native_endian, makepath
 from desiutil.log import get_logger
 
-
 def write_frame(outfile, frame, header=None, fibermap=None, units=None):
     """Write a frame fits file and returns path to file written.
 
@@ -172,3 +171,32 @@ def read_frame(filename, nspec=None):
         log.error("Frame did not pass simple vetting test. diagnosis={:d}".format(diagnosis))
     # Return
     return frame
+
+
+def search_for_framefile(frame_file):
+    """ Search for an input frame_file in the desispec redux hierarchy
+    Args:
+        frame_file:  str
+
+    Returns:
+        mfile: str,  full path to frame_file if found else raise error
+
+    """
+    log=get_logger()
+    # Parse frame file
+    path, ifile = os.path.split(frame_file)
+    splits = ifile.split('-')
+    root = splits[0]
+    camera = splits[1]
+    fexposure = int(splits[2].split('.')[0])
+
+    # Loop on nights
+    nights = get_nights()
+    for night in nights:
+        for exposure in get_exposures(night):
+            if exposure == fexposure:
+                mfile = findfile(root, camera=camera, night=night, expid=exposure)
+                if os.path.isfile(mfile):
+                    return mfile
+                else:
+                    log.error("Expected file {:s} not found..".format(mfile))

py/desispec/io/meta.py

@@ -300,6 +300,29 @@ def get_exposures(night, raw=False, rawdata_dir=None, specprod_dir=None):
     return sorted(exposures)
 
 
+def get_nights(strip_path=True, rawdata_dir=None, specprod_dir=None):
+    """
+    Args:
+        strip_path:  bool, optional; Strip the path to the nights folders
+        rawdata_dir:
+        specprod_dir:
+
+    Returns:
+        nights: list of nights (without or with paths)
+    """
+    # Init
+    if specprod_dir is None:
+        specprod_dir = specprod_root()
+    # Glob for nights
+    exp_path = os.path.join(specprod_dir,'exposures')
+    nights_path = glob.glob(exp_path+'/*')
+    # Strip off path?
+    if strip_path:
+        return [night_path.split('/')[-1] for night_path in nights_path]
+    else:
+        return nights_path
+
+
 def rawdata_root():
     """Returns directory root for raw data, i.e. ``$DESI_SPECTRO_DATA``
 

py/desispec/io/qa.py

@@ -125,7 +125,7 @@ def write_qa_brick(outfile, qabrick):
     return outfile
 
 
-def write_qa_frame(outfile, qaframe):
+def write_qa_frame(outfile, qaframe, verbose=False):
     """Write QA for a given frame
 
     Args:
@@ -134,6 +134,7 @@ def write_qa_frame(outfile, qaframe):
         qa_exp : QA_Frame object, with the following attributes
             qa_data: dict of QA info
     """
+    log=get_logger()
     outfile = makepath(outfile, 'qa')
 
     # Generate the dict
@@ -143,6 +144,8 @@ def write_qa_frame(outfile, qaframe):
     # Simple yaml
     with open(outfile, 'w') as yamlf:
         yamlf.write( yaml.dump(ydict))#, default_flow_style=True) )
+    if verbose:
+        log.info("Wrote QA frame file: {:s}".format(outfile))
 
     return outfile
 

py/desispec/qa/qa_frame.py

@@ -4,8 +4,7 @@
 
 from __future__ import print_function, absolute_import, division
 
-import numpy as np
-import os
+import warnings
 
 from desiutil.log import get_logger
 from desispec.io import read_params
@@ -194,3 +193,91 @@ def __repr__(self):
         """
         return ('{:s}: night={:s}, expid={:d}, camera={:s}, flavor={:s}'.format(
                 self.__class__.__name__, self.night, self.expid, self.camera, self.flavor))
+
+
+def qaframe_from_frame(frame_file, specprod_dir=None, make_plots=False):
+    """  Generate a qaframe object from an input frame_file name (and night)
+    Write QA to disk
+    Will also make plots if directed
+    Args:
+        frame_file: str
+        specprod_dir: str, optional
+        make_plots: bool, optional
+
+    Returns:
+
+    """
+    from desispec.io import read_frame
+    from desispec.io import meta
+    from desispec.io.qa import load_qa_frame, write_qa_frame
+    from desispec.io.frame import search_for_framefile
+    from desispec.io.fiberflat import read_fiberflat
+    from desispec.qa import qa_plots
+    from desispec.io.sky import read_sky
+    from desispec.io.fluxcalibration import read_flux_calibration
+
+    if '/' in frame_file:  # If present, assume full path is used here
+        pass
+    else: # Find the frame file in the desispec hierarchy?
+        frame_file = search_for_framefile(frame_file)
+    # Load frame
+    frame = read_frame(frame_file)
+    frame_meta = frame.meta
+    night = frame_meta['NIGHT'].strip()
+    camera = frame_meta['CAMERA'].strip()
+    expid = frame_meta['EXPID']
+    spectro = int(frame_meta['CAMERA'][-1])
+    if frame_meta['FLAVOR'] in ['flat', 'arc']:
+        qatype = 'qa_calib'
+    else:
+        qatype = 'qa_data'
+    # Load
+    qafile = meta.findfile(qatype, night=night, camera=camera, expid=expid, specprod_dir=specprod_dir)
+    qaframe = load_qa_frame(qafile, frame, flavor=frame.meta['FLAVOR'])
+    # Flat QA
+    if frame.meta['FLAVOR'] in ['flat']:
+        fiberflat_fil = meta.findfile('fiberflat', night=night, camera=camera, expid=expid,
+                                      specprod_dir=specprod_dir)
+        fiberflat = read_fiberflat(fiberflat_fil)
+        qaframe.run_qa('FIBERFLAT', (frame, fiberflat), clobber=True)
+        if make_plots:
+            # Do it
+            qafig = meta.findfile('qa_flat_fig', night=night, camera=camera, expid=expid, specprod_dir=specprod_dir)
+            qa_plots.frame_fiberflat(qafig, qaframe, frame, fiberflat)
+    # SkySub QA
+    if qatype == 'qa_data':
+        sky_fil = meta.findfile('sky', night=night, camera=camera, expid=expid, specprod_dir=specprod_dir)
+        try:
+            skymodel = read_sky(sky_fil)
+        except FileNotFoundError:
+            warnings.warn("Sky file {:s} not found.  Skipping..".format(sky_fil))
+        else:
+            qaframe.run_qa('SKYSUB', (frame, skymodel))
+            if make_plots:
+                qafig = meta.findfile('qa_sky_fig', night=night, camera=camera, expid=expid,
+                                      specprod_dir=specprod_dir)
+                qa_plots.frame_skyres(qafig, frame, skymodel, qaframe)
+    # FluxCalib QA
+    if qatype == 'qa_data':
+        # Standard stars
+        stdstar_fil = meta.findfile('stdstars', night=night, camera=camera, expid=expid, specprod_dir=specprod_dir,
+                                    spectrograph=spectro)
+        # try:
+        #    model_tuple=read_stdstar_models(stdstar_fil)
+        # except FileNotFoundError:
+        #    warnings.warn("Standard star file {:s} not found.  Skipping..".format(stdstar_fil))
+        # else:
+        flux_fil = meta.findfile('calib', night=night, camera=camera, expid=expid, specprod_dir=specprod_dir)
+        try:
+            fluxcalib = read_flux_calibration(flux_fil)
+        except FileNotFoundError:
+            warnings.warn("Flux file {:s} not found.  Skipping..".format(flux_fil))
+        else:
+            qaframe.run_qa('FLUXCALIB', (frame, fluxcalib))  # , model_tuple))#, indiv_stars))
+            if make_plots:
+                qafig = meta.findfile('qa_flux_fig', night=night, camera=camera, expid=expid,
+                                      specprod_dir=specprod_dir)
+                qa_plots.frame_fluxcalib(qafig, qaframe, frame, fluxcalib)  # , model_tuple)
+    # Write
+    write_qa_frame(qafile, qaframe, verbose=True)
+    return qaframe

py/desispec/qa/qa_plots.py

@@ -97,6 +97,7 @@ def frame_skyres(outfil, frame, skymodel, qaframe, quick_look=False):
     skymodel: SkyModel object
     qaframe: QAFrame object
     """
+    from desispec.sky import subtract_sky
 
     # Access metrics
     '''
@@ -108,6 +109,7 @@ def frame_skyres(outfil, frame, skymodel, qaframe, quick_look=False):
     pchi2_med = scipy.stats.chisqprob(chi2_med, dof_wavg)
     '''
     skyfibers = np.array(qaframe.qa_data['SKYSUB']["METRICS"]["SKY_FIBERID"])
+    subtract_sky(frame, skymodel)
     res=frame.flux[skyfibers]
     res_ivar=frame.ivar[skyfibers]
     if quick_look:
@@ -119,7 +121,10 @@ def frame_skyres(outfil, frame, skymodel, qaframe, quick_look=False):
 
     # Plot
     fig = plt.figure(figsize=(8, 10.0))
-    gs = gridspec.GridSpec(4,2)
+    if quick_look:
+        gs = gridspec.GridSpec(4,2)
+    else:
+        gs = gridspec.GridSpec(2,2)
     xmin,xmax = np.min(frame.wave), np.max(frame.wave)
 
     # Simple residual plot