Trace shifts in pipeline

py/desispec/pipeline/db.py

@@ -190,6 +190,11 @@ def all_tasks(night, nside):
                     props["spec"] = spec
                     props["expid"] = int(ex)
                     props["state"] = "waiting" # see defs.task_states
+
+                    # Add traceshift
+                    full["traceshift"].append(props)
+
+                    # Add extractions
                     full["extract"].append(props)
 
                 if flavor == "flat" :

py/desispec/pipeline/tasks/__init__.py

@@ -61,6 +61,6 @@
                 .format(taskname))
         # add the class to the dictionary.
         base.task_classes[taskname] = taskclass()
-    base.default_task_chain = ["pix", "psf", "psfnight", "extract",
+    base.default_task_chain = ["pix", "psf", "psfnight", "traceshift", "extract",
         "fiberflat", "fiberflatnight", "sky", "starfit", "fluxcalib",
         "cframe", "spectra", "redshift"]

py/desispec/pipeline/tasks/extract.py

@@ -64,7 +64,7 @@ def _deps(self, name, db, inputs):
         deptasks = {
             "input" : task_classes["pix"].name_join(props),
             "fibermap" : task_classes["fibermap"].name_join(props),
-            "psf" : task_classes["psfnight"].name_join(props)
+            "psf" : task_classes["traceshift"].name_join(props)
             }
         return deptasks
 

py/desispec/pipeline/tasks/psfnight.py

@@ -149,7 +149,7 @@ def postprocessing(self, db, name, cur):
         # run getready for all extraction with same night,band,spec
         props = self.name_split(name)
         log  = get_logger()
-        tt  = "extract"
+        tt  = "traceshift"
         cmd = "select name from {} where night={} and band='{}' and spec={} and state=0".format(tt,props["night"],props["band"],props["spec"])
         cur.execute(cmd)
         tasks = [ x for (x,) in cur.fetchall() ]

py/desispec/pipeline/tasks/traceshift.py

@@ -0,0 +1,142 @@
+#
+# See top-level LICENSE.rst file for Copyright information
+#
+# -*- coding: utf-8 -*-
+
+from __future__ import absolute_import, division, print_function
+
+from collections import OrderedDict
+
+from ..defs import (task_name_sep, task_state_to_int, task_int_to_state)
+
+from ...util import option_list
+
+from ...io import findfile
+
+from .base import (BaseTask, task_classes)
+
+from desiutil.log import get_logger
+
+import sys,re,os,copy
+
+# NOTE: only one class in this file should have a name that starts with "Task".
+
+class TaskTraceShift(BaseTask):
+    """Class containing the properties of one trace shift task.
+    """
+    def __init__(self):
+        super(TaskTraceShift, self).__init__()
+        # then put int the specifics of this class
+        # _cols must have a state
+        self._type = "traceshift"
+        self._cols = [
+            "night",
+            "band",
+            "spec",
+            "expid",
+            "state"
+        ]
+        self._coltypes = [
+            "integer",
+            "text",
+            "integer",
+            "integer",
+            "integer"
+        ]
+        # _name_fields must also be in _cols
+        self._name_fields  = ["night","band","spec","expid"]
+        self._name_formats = ["08d","s","d","08d"]
+
+    def _paths(self, name):
+        """See BaseTask.paths.
+        """
+        props = self.name_split(name)
+        camera = "{}{}".format(props["band"], props["spec"])
+        return [ findfile("psf", night=props["night"], expid=props["expid"],
+            camera=camera, groupname=None, nside=None, band=props["band"],
+            spectrograph=props["spec"]) ]
+
+    def _deps(self, name, db, inputs):
+        """See BaseTask.deps.
+        """
+        from .base import task_classes
+        props = self.name_split(name)
+        deptasks = {
+            "image" : task_classes["pix"].name_join(props),
+            "psf" : task_classes["psfnight"].name_join(props)
+            }
+        return deptasks
+
+    def _run_max_procs(self, procs_per_node):
+        """See BaseTask.run_max_procs.
+        """
+        return 1 
+
+
+    def _run_time(self, name, procs_per_node, db=None):
+        """See BaseTask.run_time.
+        """
+        return 15 # 8 minute per bundle of 25 fibers on edison, but can be slower
+
+
+    def _run_defaults(self):
+        """See BaseTask.run_defaults.
+        """
+        opts = {}
+        opts["degxx"] = 2
+        opts["degxy"] = 2
+        opts["degyx"] = 2
+        opts["degyy"] = 2
+        opts["auto"]  = True
+        return opts
+
+
+    def _option_list(self, name, opts):
+        """Build the full list of options.
+
+        This includes appending the filenames and incorporating runtime
+        options.
+        """
+        from .base import task_classes, task_type
+
+        deps = self.deps(name)
+        options = {}
+        options["image"]    = task_classes["pix"].paths(deps["image"])[0]
+        options["psf"]      = task_classes["psfnight"].paths(deps["psf"])[0]
+        options["outpsf"]   = self.paths(name)[0]
+
+        options.update(opts)
+        return option_list(options)
+
+    def _run_cli(self, name, opts, procs, db):
+        """See BaseTask.run_cli.
+        """
+        entry = "desi_compute_trace_shifts"
+        optlist = self._option_list(name, opts)
+        com = "{} {}".format(entry, " ".join(optlist))
+        return com
+
+    def _run(self, name, opts, comm, db):
+        """See BaseTask.run.
+        """
+        from ...scripts import trace_shifts
+        optlist = self._option_list(name, opts)
+        args = trace_shifts.parse(optlist)
+        if comm is None :
+            trace_shifts.main(args)
+        else :
+            trace_shifts.main_mpi(args, comm=comm)
+        return
+
+    def postprocessing(self, db, name, cur):
+        """For successful runs, postprocessing on DB"""
+        # run getready for all extraction with same night,band,spec
+        props = self.name_split(name)
+        log  = get_logger()
+        for tt in ["extract"] :
+            cmd = "select name from {} where night={} and expid={} and band='{}' and spec={} and state=0".format(tt,props["night"],props["expid"],props["band"],props["spec"])
+            cur.execute(cmd)
+            tasks = [ x for (x,) in cur.fetchall() ]
+            log.debug("checking {} {}".format(tt,tasks))
+            for task in tasks :
+                task_classes[tt].getready( db=db,name=task,cur=cur)

py/desispec/scripts/trace_shifts.py

@@ -7,6 +7,7 @@
 from numpy.linalg.linalg import LinAlgError
 import astropy.io.fits as pyfits
 from numpy.polynomial.legendre import legval,legfit
+from pkg_resources import resource_exists, resource_filename
 
 from desispec.io import read_image
 from desiutil.log import get_logger
@@ -31,7 +32,10 @@ def parse(options=None):
     parser.add_argument('--lines', type = str, default = None, required=False,
                         help = 'path of lines ASCII file. Using this option changes the fit method.')
     parser.add_argument('--spectrum', type = str, default = None, required=False,
-                        help = 'path to a spectrum ASCII file (e.g. the DESIMODEL sky spectrum)')
+                        help = 'path to an spectrum ASCII file for external wavelength calibration')
+    parser.add_argument('--sky', action = 'store_true',
+                        help = 'use sky spectrum desispec/data/spec-sky.dat for external wavelength calibration')
+
     parser.add_argument('--outpsf', type = str, default = None, required=True,
                         help = 'path of output PSF with shifted traces')
     parser.add_argument('--outoffsets', type = str, default = None, required=False,
@@ -46,6 +50,9 @@ def parse(options=None):
                         help = 'polynomial degree for y shifts along y')
     parser.add_argument('--continuum', action='store_true',
                         help = 'only fit shifts along x for continuum input image')
+    parser.add_argument('--auto', action='store_true',
+                        help = 'choose best method (sky,continuum or just internal calib) from the FLAVOR keyword in the input image header')
+
     parser.add_argument('--nfibers', type = int, default = None, required=False,
                         help = 'limit the number of fibers for debugging')
     parser.add_argument('--max-error', type = float, default = 0.05 , required=False,
@@ -98,8 +105,40 @@ def main(args) :
     else :
         log.info("We will do an internal calibration of trace coordinates without using the psf shape in a first step")
 
-
-
+
+    internal_wavelength_calib    = (not args.continuum)
+    external_wavelength_calib   = args.sky | ( args.spectrum is not None )
+
+    if args.auto :
+        log.debug("read flavor of input image {}".format(args.image))
+        hdus = pyfits.open(args.image)
+        if "FLAVOR" not in hdus[0].header :
+            log.error("no FLAVOR keyword in image header, cannot run with --auto option")
+            raise KeyError("no FLAVOR keyword in image header, cannot run with --auto option")
+        flavor = hdus[0].header["FLAVOR"].strip().lower()
+        hdus.close()
+        log.info("Input is a '{}' image".format(flavor))
+        if flavor == "flat" : 
+            internal_wavelength_calib = False
+            external_wavelength_calib = False
+        elif flavor == "arc" :
+            internal_wavelength_calib = True
+            external_wavelength_calib = False
+        else :
+            internal_wavelength_calib = True
+            external_wavelength_calib = True
+        log.info("wavelength calib, internal={}, external={}".format(internal_wavelength_calib,external_wavelength_calib))
+
+    spectrum_filename = args.spectrum
+    if external_wavelength_calib and spectrum_filename is None :
+        srch_file = "data/spec-sky.dat"
+        if not resource_exists('desispec', srch_file):
+            log.error("Cannot find sky spectrum file {:s}".format(srch_file))
+            raise RuntimeError("Cannot find sky spectrum file {:s}".format(srch_file))
+        else :
+            spectrum_filename=resource_filename('desispec', srch_file)
+            log.info("Use external calibration from cross-correlation with {}".format(spectrum_filename))
+
     if args.nfibers is not None :
         nfibers = args.nfibers # FOR DEBUGGING
 
@@ -129,7 +168,7 @@ def main(args) :
 
         # measure x shifts
         x_for_dx,y_for_dx,dx,ex,fiber_for_dx,wave_for_dx = compute_dx_from_cross_dispersion_profiles(xcoef,ycoef,wavemin,wavemax, image=image, fibers=fibers, width=7, deg=args.degxy)
-        if not args.continuum :
+        if internal_wavelength_calib :
             # measure y shifts
             x_for_dy,y_for_dy,dy,ey,fiber_for_dy,wave_for_dy = compute_dy_using_boxcar_extraction(xcoef,ycoef,wavemin,wavemax, image=image, fibers=fibers, width=7)
         else :
@@ -188,7 +227,7 @@ def main(args) :
 
             if degxy>0 and degyy>0 and degxy>degxx and degyy>degyx : # first along wavelength
                 if degxy>0 : degxy-=1
-                if degxy>0 : degyy-=1
+                if degyy>0 : degyy-=1
             else : # then along fiber
                 if degxx>0 : degxx-=1
                 if degyx>0 : degyx-=1
@@ -227,15 +266,15 @@ def main(args) :
 
 
     # use an input spectrum as an external calibration of wavelength
-    if args.spectrum :
-         
+    if spectrum_filename  :
+
 
         log.info("write and reread PSF to be sure predetermined shifts were propagated")
         write_traces_in_psf(args.psf,args.outpsf,xcoef,ycoef,wavemin,wavemax)
         psf,xcoef,ycoef,wavemin,wavemax = read_psf_and_traces(args.outpsf)
 
         ycoef=shift_ycoef_using_external_spectrum(psf=psf,xcoef=xcoef,ycoef=ycoef,wavemin=wavemin,wavemax=wavemax,
-                                                  image=image,fibers=fibers,spectrum_filename=args.spectrum,degyy=args.degyy,width=7)
+                                                  image=image,fibers=fibers,spectrum_filename=spectrum_filename,degyy=args.degyy,width=7)
 
         write_traces_in_psf(args.psf,args.outpsf,xcoef,ycoef,wavemin,wavemax)
         log.info("wrote modified PSF in %s"%args.outpsf)

py/desispec/trace_shifts.py

@@ -360,10 +360,9 @@ def compute_dy_from_spectral_cross_correlation(flux,wave,refflux,ivar=None,hw=3.
         if c[0]>0 :
             delta=-c[1]/(2.*c[0])
             sigma=np.sqrt(1./c[0] + error_floor**2)
-            if ndata>1 :
-                chi2pdf=(c[0]*delta**2+c[1]*delta+c[2])/(ndata+1)
-                if chi2pdf>1 : sigma *= np.sqrt(chi2pdf) 
-
+            # do not scale error bars using chi2pdf because the
+            # two spectra do not necessarily match
+            # (for instance dark vs bright time sky spectrum)                
         else :
             # something else went wrong
             delta=0.
@@ -767,7 +766,6 @@ def shift_ycoef_using_external_spectrum(psf,xcoef,ycoef,wavemin,wavemax,image,fi
         eps=0.1
         x,yp=psf.xy(fiber_for_psf_evaluation,bin_wave+eps)
         dydw=(yp-y)/eps
-
         if err*dydw<1 :
             dy=np.append(dy,-dwave*dydw)
             ey=np.append(ey,err*dydw)