Scattered light model

doc/releases/1.14.1dev.rst

@@ -28,6 +28,8 @@ Script Changes
 
 - When making datacubes, the user can select a separate frame to use for the sky subtraction.
   In this case, it is the processed data that will be used for sky subtraction (akin to nodding).
+- A new script has been written (``chk_scattlight``) to check the generation of the scattered
+  light calibration frame model.
 
 Datamodel Changes
 -----------------

doc/scripts/build_datacontainer_datamodels.py

@@ -67,6 +67,7 @@ def build_datamodel_tbl(obj):
     from pypeit.flatfield import FlatImages
     from pypeit.manual_extract import ManualExtractionObj
     from pypeit.onespec import OneSpec
+    from pypeit.scattlight import ScatteredLight
     from pypeit.sensfunc import SensFunc
     from pypeit.slittrace import SlitTraceSet
     from pypeit.spec2dobj import Spec2DObj
@@ -86,7 +87,7 @@ def build_datamodel_tbl(obj):
 
     from pypeit.images import buildimage
 
-    datacontainers = [Alignments, EdgeTraceSet, FlatImages, ManualExtractionObj, OneSpec,
+    datacontainers = [Alignments, EdgeTraceSet, FlatImages, ManualExtractionObj, OneSpec, ScatteredLight,
                       SensFunc, SlitTraceSet, Spec2DObj, SpecObj, TracePCA, WaveCalib, WaveTilts,
                       bspline, DataCube, PypeItFit, MultiSlitFlexure, Telluric, DetectorContainer,
                       Mosaic, PypeItImage, WaveFit]

doc/spectrographs/keck_kcwi.rst

@@ -112,11 +112,16 @@ Scattered Light Removal
 KCWI suffers from mild scattered light (at the level of ~1 percent),
 and this appears to be worse near regions of the detector where there
 is brighter illumination. We are currently working towards building a
-full model of the scattered light. For the moment, PypeIt uses a robust
-piecewise polynomial to model the scattered light that is detected on
-the left of slice 1, the unilluminated region between slices 12 and 13,
-and the right of slice 24. The model is smooth and continuous, and is
-determined for each spectral pixel. By default, the scattered light is
+full model of the scattered light. For KCWI, the main contributor to
+the scattered light is referred to as the "narcissistic ghost" by
+Morrissey et al. (2018), ApJ, 864, 93. This scattered light is thought
+to be a reflection off the detector that travels back through the optical
+system. Some fraction gets sent back out to space, while the remainder
+comes back through the optical system and a fuzzy version of this is
+re-imaged onto the detector. The current KCWI scattered light model is
+designed to account for these effects. To generate a scattered light model,
+it's a good idea to use a frame that has a lot of counts (e.g. a flatfield
+frame, or a standard star). By default, the scattered light is
 subtracted from the science frame, the pixel flat, and the illumination
 flat. To turn off the scattered light subtraction, you can add the
 following lines to your :ref:`pypeit_file`:

pypeit/calibrations.py

@@ -20,13 +20,15 @@
 from pypeit import alignframe
 from pypeit import flatfield
 from pypeit import edgetrace
+from pypeit import scattlight
 from pypeit import slittrace
 from pypeit import wavecalib
 from pypeit import wavetilts
 from pypeit.calibframe import CalibFrame
 from pypeit.images import buildimage
 from pypeit.metadata import PypeItMetaData
 from pypeit.core import framematch
+from pypeit.core import parse
 from pypeit.par import pypeitpar
 from pypeit.spectrographs.spectrograph import Spectrograph
 from pypeit import io
@@ -181,6 +183,7 @@ def __init__(self, fitstbl, par, spectrograph, caldir, qadir=None,
         self.msbpm = None
         self.wv_calib = None
         self.slits = None
+        self.msscattlight = None
 
         self.wavetilts = None
         self.flatimages = None
@@ -537,6 +540,93 @@ def get_bpm(self, frame=None):
         # Return
         return self.msbpm
 
+    def get_scattlight(self):
+        """
+        Load or generate the scattered light model.
+
+        Returns:
+            :class:`~pypeit.scattlight.ScatteredLight`: The processed calibration image including the model.
+        """
+        # Check for existing data
+        if not self._chk_objs(['msbpm', 'slits']):
+            msgs.warn('Must have the bpm and the slits defined to make a scattered light image!  '
+                      'Skipping and may crash down the line')
+            return self.msscattlight
+
+        # Check internals
+        self._chk_set(['det', 'calib_ID', 'par'])
+
+        # Prep
+        frame = {'type': 'scattlight', 'class': scattlight.ScatteredLight}
+        raw_scattlight_files, cal_file, calib_key, setup, calib_id, detname = \
+            self.find_calibrations(frame['type'], frame['class'])
+        scatt_idx = self.fitstbl.find_frames(frame['type'], calib_ID=self.calib_ID, index=True)
+
+        if len(raw_scattlight_files) == 0 and cal_file is None:
+            msgs.warn(f'No raw {frame["type"]} frames found and unable to identify a relevant '
+                      'processed calibration frame.  Continuing...')
+            return self.msscattlight
+
+        # If a processed calibration frame exists and we want to reuse it, do
+        # so:
+        if cal_file.exists() and self.reuse_calibs:
+            self.msscattlight = frame['class'].from_file(cal_file)
+            return self.msscattlight
+
+        # Scattered light model does not exist or we're not reusing it.
+        # Need to build everything from scratch.  Start with the trace image.
+        msgs.info('Creating edge tracing calibration frame using files: ')
+        for f in raw_scattlight_files:
+            msgs.prindent(f'{Path(f).name}')
+
+        # Reset the BPM
+        self.get_bpm(frame=raw_scattlight_files[0])
+
+        binning = self.fitstbl[scatt_idx[0]]['binning']
+        dispname = self.fitstbl[scatt_idx[0]]['dispname']
+        scattlightImage = buildimage.buildimage_fromlist(self.spectrograph, self.det,
+                                                         self.par['scattlightframe'], raw_scattlight_files,
+                                                         bias=self.msbias, bpm=self.msbpm,
+                                                         dark=self.msdark, calib_dir=self.calib_dir,
+                                                         setup=setup, calib_id=calib_id)
+
+        spatbin = parse.parse_binning(binning)[1]
+        pad = self.par['scattlight']['pad'] // spatbin
+        offslitmask = self.slits.slit_img(pad=pad, initial=True, flexure=None) == -1
+
+        model, modelpar, success = self.spectrograph.scattered_light(scattlightImage.image, offslitmask,
+                                                                     binning=binning, dispname=dispname)
+
+        if not success:
+            # Something went awry
+            msgs.warn('Scattered light modelling failed.  Continuing, but likely to fail soon...')
+            self.success = False
+            return self.msscattlight
+
+        # Now generate the DataModel
+        self.msscattlight = scattlight.ScatteredLight(PYP_SPEC=self.spectrograph.name,
+                                                      pypeline=self.spectrograph.pypeline,
+                                                      detname=scattlightImage.detector.name,
+                                                      nspec=scattlightImage.shape[0], nspat=scattlightImage.shape[1],
+                                                      binning=scattlightImage.detector.binning,
+                                                      pad=self.par['scattlight']['pad'],
+                                                      scattlight_raw=scattlightImage,
+                                                      scattlight_model=model,
+                                                      scattlight_param=modelpar)
+
+        # TODO :: Should we go back and recalculate the slit edges once the scattered light is known?
+
+        if self.msscattlight is not None:
+            # Show the result if requested
+            if self.show:
+                self.msscattlight.show()
+
+            # Save the master scattered light model
+            self.msscattlight.set_paths(self.calib_dir, setup, calib_id, detname)
+            self.msscattlight.to_file()
+
+        return self.msscattlight
+
     def get_flats(self):
         """
         Load or generate the flat-field calibration images.
@@ -622,7 +712,8 @@ def get_flats(self):
             pixel_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                         self.par['pixelflatframe'],
                                                         raw_pixel_files, dark=self.msdark,
-                                                        bias=self.msbias, bpm=self.msbpm)
+                                                        bias=self.msbias, bpm=self.msbpm,
+                                                        scattlight=self.msscattlight)
             if len(raw_lampoff_files) > 0:
                 # Reset the BPM
                 self.get_bpm(frame=raw_lampoff_files[0])
@@ -633,7 +724,7 @@ def get_flats(self):
                                                               self.par['lampoffflatsframe'],
                                                               raw_lampoff_files,
                                                               dark=self.msdark, bias=self.msbias,
-                                                              bpm=self.msbpm)
+                                                              bpm=self.msbpm, scattlight=self.msscattlight)
                 pixel_flat = pixel_flat.sub(lampoff_flat)
 
             # Initialise the pixel flat
@@ -656,7 +747,7 @@ def get_flats(self):
                 msgs.prindent(f'{Path(f).name}')
             illum_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                         self.par['illumflatframe'], raw_illum_files,
-                                                        dark=self.msdark, bias=self.msbias,
+                                                        dark=self.msdark, bias=self.msbias, scattlight=self.msscattlight,
                                                         flatimages=self.flatimages, bpm=self.msbpm)
             if len(raw_lampoff_files) > 0:
                 msgs.info('Subtracting lamp off flats using files: ')
@@ -667,7 +758,9 @@ def get_flats(self):
                                                                   self.par['lampoffflatsframe'],
                                                                   raw_lampoff_files,
                                                                   dark=self.msdark,
-                                                                  bias=self.msbias, bpm=self.msbpm)
+                                                                  bias=self.msbias,
+                                                                  scattlight=self.msscattlight,
+                                                                  bpm=self.msbpm)
                 illum_flat = illum_flat.sub(lampoff_flat)
 
             # Initialise the pixel flat
@@ -776,6 +869,7 @@ def get_slits(self):
         traceImage = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                     self.par['traceframe'], raw_trace_files,
                                                     bias=self.msbias, bpm=self.msbpm,
+                                                    scattlight=self.msscattlight,
                                                     dark=self.msdark, calib_dir=self.calib_dir,
                                                     setup=setup, calib_id=calib_id)
         if len(raw_lampoff_files) > 0:
@@ -789,7 +883,8 @@ def get_slits(self):
             lampoff_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det,
                                                           self.par['lampoffflatsframe'],
                                                           raw_lampoff_files, dark=self.msdark,
-                                                          bias=self.msbias, bpm=self.msbpm)
+                                                          bias=self.msbias, scattlight=self.msscattlight,
+                                                          bpm=self.msbpm)
             traceImage = traceImage.sub(lampoff_flat)
 
         edges = edgetrace.EdgeTraceSet(traceImage, self.spectrograph, self.par['slitedges'],
@@ -1254,7 +1349,7 @@ def default_steps():
         # Order matters!  And the name must match a viable "get_{step}" method
         # in Calibrations.
         # TODO: Does the bpm need to be done after the dark?
-        return ['bias', 'dark', 'bpm', 'slits', 'arc', 'tiltimg', 'wv_calib', 'tilts', 'flats']
+        return ['bias', 'dark', 'bpm', 'slits', 'arc', 'tiltimg', 'wv_calib', 'tilts', 'scattlight', 'flats']
 
 
 class IFUCalibrations(Calibrations):
@@ -1272,7 +1367,7 @@ def default_steps():
         """
         # Order matters!
         return ['bias', 'dark', 'bpm', 'arc', 'tiltimg', 'slits', 'wv_calib', 'tilts', 'align',
-                'flats']
+                'scattlight', 'flats']
 
 
 def check_for_calibs(par, fitstbl, raise_error=True, cut_cfg=None):

pypeit/core/framematch.py

@@ -29,9 +29,10 @@ def __init__(self):
                      ('pinhole', 'Pinhole observation used for tracing slit centers'),
                    ('pixelflat', 'Flat-field exposure used for pixel-to-pixel response'),
                    ('illumflat', 'Flat-field exposure used for illumination flat'),
-                   ('lampoffflats', 'Flat-field exposure with lamps off used to remove '
-                                   'persistence from lamp on flat exposures and/or thermal emission '
-                                   'from the telescope and dome'),
+                ('lampoffflats', 'Flat-field exposure with lamps off used to remove '
+                                 'persistence from lamp on flat exposures and/or thermal emission '
+                                 'from the telescope and dome'),
+                  ('scattlight', 'Frame (ideally with lots of counts) used to determine the scattered light model'),
                      ('science', 'On-sky observation of a primary target'),
                     ('standard', 'On-sky observation of a flux calibrator'),
                        ('trace', 'High-count exposure used to trace slit positions'),

pypeit/display/display.py

@@ -636,6 +636,32 @@ def show_tilts(viewer, ch, tilt_traces, yoff=0., xoff=0., points=True, nspec=Non
     canvas.add('constructedcanvas', canvas_list)
 
 
+def show_scattered_light(image_list, slits=None, wcs_match=True):
+    """
+    Interface to ginga to show the quality of the Scattered Light subtraction
 
-
+    Parameters
+    ----------
+    image_list : zip
+        A zip of the images to show, their names, and the scales
+    slits : :class:`~pypeit.slittrace.SlitTraceSet`, optional
+        The current slit traces
+    wcs_match : :obj:`bool`, optional
+        Use a reference image for the WCS and match all image in other channels to it.
+    """
+    connect_to_ginga(raise_err=True, allow_new=True)
+    if slits is not None:
+        left, right, mask = slits.select_edges()
+        gpm = mask == 0
+    # Loop me
+    clear = True
+    for img, name, cut in image_list:
+        if img is None:
+            continue
+        viewer, ch = show_image(img, chname=name, cuts=cut, wcs_match=wcs_match, clear=clear)
+        if slits is not None:
+            show_slits(viewer, ch, left[:, gpm], right[:, gpm], slit_ids=slits.spat_id[gpm])
+        # Turn off clear
+        if clear:
+            clear = False
 

pypeit/images/buildimage.py

@@ -77,6 +77,18 @@ class TiltImage(pypeitimage.PypeItCalibrationImage):
     calib_type = 'Tiltimg'
 
 
+class ScatteredLightImage(pypeitimage.PypeItCalibrationImage):
+    """
+    Simple DataContainer for the Scattered Light Image
+    """
+    # version is inherited from PypeItImage
+
+    # I/O
+    output_to_disk = ('SCATTLIGHT_IMAGE', 'SCATTLIGHT_FULLMASK', 'SCATTLIGHT_DETECTOR')
+    hdu_prefix = 'SCATTLIGHT_'
+    calib_type = 'ScattLight'
+
+
 class TraceImage(pypeitimage.PypeItCalibrationImage):
     """
     Simple DataContainer for the Trace Image
@@ -133,6 +145,7 @@ def construct_file_name(cls, calib_key, calib_dir=None, basename=None):
     arc=ArcImage,
     tilt=TiltImage,
     trace=TraceImage,
+    scattlight=ScatteredLightImage,
     align=AlignImage)
 """
 The list of classes that :func:`buildimage_fromlist` should use to decorate the
@@ -144,7 +157,7 @@ def construct_file_name(cls, calib_key, calib_dir=None, basename=None):
 
 
 def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=None, dark=None,
-                        flatimages=None, maxiters=5, ignore_saturation=True, slits=None,
+                        scattlight=None, flatimages=None, maxiters=5, ignore_saturation=True, slits=None,
                         mosaic=None, calib_dir=None, setup=None, calib_id=None):
     """
     Perform basic image processing on a list of images and combine the results.
@@ -179,6 +192,8 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
         dark (:class:`~pypeit.images.buildimage.DarkImage`, optional):
             Dark-current image; passed directly to
             :func:`~pypeit.images.rawimage.RawImage.process` for all images.
+        scattlight (:class:`~pypeit.scattlight.ScatteredLight`, optional):
+            Scattered light model to be used to determine scattered light.
         flatimages (:class:`~pypeit.flatfield.FlatImages`, optional):
             Flat-field images for flat fielding; passed directly to
             :func:`~pypeit.images.rawimage.RawImage.process` for all images.
@@ -234,7 +249,7 @@ def buildimage_fromlist(spectrograph, det, frame_par, file_list, bias=None, bpm=
 
     # Do it
     combineImage = combineimage.CombineImage(spectrograph, det, frame_par['process'], file_list)
-    pypeitImage = combineImage.run(bias=bias, bpm=bpm, dark=dark, flatimages=flatimages,
+    pypeitImage = combineImage.run(bias=bias, bpm=bpm, dark=dark, flatimages=flatimages, scattlight=scattlight,
                                    sigma_clip=frame_par['process']['clip'],
                                    sigrej=frame_par['process']['comb_sigrej'],
                                    maxiters=maxiters, ignore_saturation=ignore_saturation,

pypeit/images/combineimage.py

@@ -62,9 +62,8 @@ def __init__(self, spectrograph, det, par, files):
         if self.nfiles == 0:
             msgs.error('CombineImage requires a list of files to instantiate')
 
-    def run(self, bias=None, flatimages=None, ignore_saturation=False, sigma_clip=True,
-            bpm=None, sigrej=None, maxiters=5, slits=None, dark=None, combine_method='mean',
-            mosaic=False):
+    def run(self, bias=None, scattlight=None, flatimages=None, ignore_saturation=False, sigma_clip=True,
+            bpm=None, sigrej=None, maxiters=5, slits=None, dark=None, combine_method='mean', mosaic=False):
         r"""
         Process and combine all images.
 
@@ -137,6 +136,8 @@ def run(self, bias=None, flatimages=None, ignore_saturation=False, sigma_clip=Tr
             bias (:class:`~pypeit.images.buildimage.BiasImage`, optional):
                 Bias image for bias subtraction; passed directly to
                 :func:`~pypeit.images.rawimage.RawImage.process` for all images.
+            scattlight (:class:`~pypeit.scattlight.ScatteredLight`, optional):
+                Scattered light model to be used to determine scattered light.
             flatimages (:class:`~pypeit.flatfield.FlatImages`, optional):
                 Flat-field images for flat fielding; passed directly to
                 :func:`~pypeit.images.rawimage.RawImage.process` for all images.
@@ -194,7 +195,7 @@ def run(self, bias=None, flatimages=None, ignore_saturation=False, sigma_clip=Tr
             # Load raw image
             rawImage = rawimage.RawImage(ifile, self.spectrograph, self.det)
             # Process
-            pypeitImage = rawImage.process(self.par, bias=bias, bpm=bpm, dark=dark,
+            pypeitImage = rawImage.process(self.par, scattlight=scattlight, bias=bias, bpm=bpm, dark=dark,
                                            flatimages=flatimages, slits=slits, mosaic=mosaic)
 
             if self.nfiles == 1: