KCWI scattered light

doc/help/run_pypeit.rst

@@ -4,7 +4,7 @@
     usage: run_pypeit [-h] [-v VERBOSITY] [-r REDUX_PATH] [-m] [-s] [-o] [-c]
                       pypeit_file
 
-    ##  PypeIt : The Python Spectroscopic Data Reduction Pipeline v1.14.1.dev9+gb73d74939
+    ##  PypeIt : The Python Spectroscopic Data Reduction Pipeline v1.13.1.dev308+g83230725a
     ##  
     ##  Available spectrographs include:
     ##   bok_bc, gemini_flamingos1, gemini_flamingos2, gemini_gmos_north_e2v,

doc/releases/1.14.1dev.rst

@@ -8,12 +8,18 @@ Dependency Changes
 Functionality/Performance Improvements and Additions
 ----------------------------------------------------
 
+- Started the development of instrument-specific scattered light removal. In this
+  release, we only model KCWI/KCRM scattered light.
+
 Instrument-specific Updates
 ---------------------------
 
 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).
+
 Datamodel Changes
 -----------------
 

doc/scripts/make_example_files.py

@@ -204,6 +204,7 @@ def make_meta_examples():
     if otmp.exists():
         otmp.unlink()
 
+
 if __name__ == '__main__':
     t = time.perf_counter()
     print('Making shane_kast_blue_A.pypeit.rst')

doc/spectrographs/keck_kcwi.rst

@@ -106,6 +106,32 @@ overscan regions. Also note that this pattern noise is different
 from the detector structure mentioned above for pixelflats. The
 pattern noise is additive, the detector structure is multiplicative.
 
+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
+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`:
+
+.. code-block:: ini
+
+    [scienceframe]
+        [[process]]
+            subtract_scattlight = False
+    [calibrations]
+        [[pixelflatframe]]
+            [[[process]]]
+                subtract_scattlight = False
+
+
 Relative spectral illumination correction
 -----------------------------------------
 

pypeit/bspline/bspline.py

@@ -74,7 +74,7 @@ class bspline(datamodel.DataContainer):
     Attributes
     ----------
     breakpoints
-        Breakpoints for bspline, spacing for these breakpoints are determinated by keywords inputs;
+        Breakpoints for bspline, spacing for these breakpoints are determined by keywords inputs;
     nord
         Order of bspline; [default=4]
     npoly

pypeit/core/datacube.py

@@ -1881,7 +1881,7 @@ def coadd_cube(files, opts, spectrograph=None, parset=None, overwrite=False):
         except:
             msgs.error("Could not load skysub image from spec2d file:" + msgs.newline() + cubepar['skysub_frame'])
         skysub_default = cubepar['skysub_frame']
-        skyImgDef = spec2DObj.skymodel/skysub_exptime  # Sky counts/second
+        skyImgDef = spec2DObj.sciimg/skysub_exptime  # Sky counts/second
         skySclDef = spec2DObj.scaleimg
 
     # Load all spec2d files and prepare the data for making a datacube
@@ -1894,8 +1894,8 @@ def coadd_cube(files, opts, spectrograph=None, parset=None, overwrite=False):
 
         # Load the header
         hdr = spec2DObj.head0
-        ifu_ra = np.append(ifu_ra, spec.compound_meta([hdr], 'ra'))
-        ifu_dec = np.append(ifu_dec, spec.compound_meta([hdr], 'dec'))
+        ifu_ra = np.append(ifu_ra, spec.get_meta_value([hdr], 'ra'))
+        ifu_dec = np.append(ifu_dec, spec.get_meta_value([hdr], 'dec'))
 
         # Get the exposure time
         exptime = hdr['EXPTIME']
@@ -1933,7 +1933,7 @@ def coadd_cube(files, opts, spectrograph=None, parset=None, overwrite=False):
                     this_skysub = "image"  # Use the current spec2d for sky subtraction
                 else:
                     skyImg = skyImgDef.copy() * exptime
-                    skyScl = skySclDef.copy() * exptime
+                    skyScl = skySclDef.copy()
                     this_skysub = skysub_default  # Use the global value for sky subtraction
             elif opts['skysub_frame'][ff].lower() == 'image':
                 skyImg = spec2DObj.skymodel
@@ -1953,7 +1953,8 @@ def coadd_cube(files, opts, spectrograph=None, parset=None, overwrite=False):
                     msgs.error("Could not load skysub image from spec2d file:" + msgs.newline() + opts['skysub_frame'][ff])
                 # TODO :: Consider allowing the actual frame (instead of the skymodel) to be used as the skysub image - make sure the BPM is carried over.
                 #      :: Allow sky data fitting (i.e. scale the flux of a skysub frame to the science frame data)
-                skyImg = spec2DObj_sky.skymodel * exptime / skysub_exptime  # Sky counts
+                skyImg = spec2DObj_sky.sciimg * exptime / skysub_exptime  # Sky counts
+                # skyImg = spec2DObj_sky.skymodel * exptime / skysub_exptime  # Sky counts
                 skyScl = spec2DObj_sky.scaleimg
                 this_skysub = opts['skysub_frame'][ff]  # User specified spec2d for sky subtraction
         if this_skysub == "none":

pypeit/core/flat.py

@@ -325,7 +325,7 @@ def illum_profile_spectral_poly(rawimg, waveimg, slitmask, slitmask_trim, model,
     gpm = gpmask if (gpmask is not None) else np.ones_like(rawimg, dtype=bool)
     # Extract the list of  spatial IDs from the slitmask
     slitmask_spatid = np.unique(slitmask)
-    slitmask_spatid = np.sort(slitmask_spatid[slitmask_spatid != 0])
+    slitmask_spatid = np.sort(slitmask_spatid[slitmask_spatid > 0])
     # Initialise the scale image that will be returned
     scaleImg = np.ones_like(rawimg)
     # Divide the slit into several bins and calculate the median of each bin

pypeit/core/parse.py

@@ -101,16 +101,14 @@ def binning2string(binspectral, binspatial):
     convention order, spectral then spatial.
 
     Args:
-        binspectral (:obj:`int`):
-            Number of on-detector pixels binned in the spectral
+        binspectral (:obj:`int`): Number of on-detector pixels binned in the spectral
             direction (along the first axis in the PypeIt convention).
-        binspatial (int):
-            Number of on-detector pixels binned in the spatial direction
+        binspatial (:obj:`int`): Number of on-detector pixels binned in the spatial direction
             (along the second axis in the PypeIt convention).
 
     Returns:
         str: Comma-separated binning along the spectral and spatial
-        directions; e.g., '2,1'
+        directions; e.g., ``2,1``
     """
     return '{0},{1}'.format(binspectral, binspatial)
 
@@ -123,7 +121,28 @@ def parse_binning(binning:str):
     parsed directly from the Header.  The developer needs to react accordingly..
 
     Args:
-        binning (str, `numpy.ndarray`_, tuple):
+        binning (:obj:`str`, `numpy.ndarray`_, :obj:`tuple`): The spectral and spatial binning.
+            Several formats are supported, including the following examples. Note that in all
+            examples, the binning in the spectral direction is 2, and the binning in the
+            spatial direction is 1:
+
+            - string format
+
+                *  comma delimited string (e.g. ``2,1``)
+
+                *  x delimited string (e.g. ``2x1``)
+
+                *  space delimited string (e.g. ``2 1``)
+
+                *  ``'None'`` will always assume 1x1 binning
+
+            - tuple format
+
+                * this must be of the form of tuple, for example: ``(2,1)``
+
+            - numpy array
+
+                * this must be of the form of tuple, for example: ``np.array([2,1])``
 
     Returns:
         tuple: binspectral, binspatial as integers
@@ -176,7 +195,7 @@ def sec2slice(subarray, one_indexed=False, include_end=False, require_dim=None,
     a list of slice objects.
 
     Args:
-        subarray (str):
+        subarray (:obj:`str`):
             The string to convert.  Should have the form of normal slice
             operation, 'start:stop:step'.  The parser ignores whether or
             not the string has the brackets '[]', but the string must

pypeit/find_objects.py

@@ -1244,6 +1244,10 @@ def global_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0),
         correction using the sky spectrum, if requested. See Reduce.global_skysub()
         for parameter definitions.
         """
+        if self.par['reduce']['findobj']['skip_skysub']:
+            msgs.info("Skipping global sky sub as per user request")
+            return np.zeros_like(self.sciImg.image)
+
         # Generate a global sky sub for all slits separately
         global_sky_sep = super().global_skysub(skymask=skymask, update_crmask=update_crmask,
                                                trim_edg=trim_edg, show_fit=show_fit, show=show,

pypeit/images/rawimage.py

@@ -177,6 +177,7 @@ def __init__(self, ifile, spectrograph, det):
                           subtract_pattern=False,
                           subtract_overscan=False,
                           subtract_continuum=False,
+                          subtract_scattlight=False,
                           trim=False,
                           orient=False,
                           subtract_bias=False,
@@ -633,6 +634,10 @@ def process(self, par, bpm=None, flatimages=None, bias=None, slits=None, dark=No
         self.spat_flexure_shift = self.spatial_flexure_shift(slits) \
                                     if self.par['spat_flexure_correct'] else None
 
+        #   - Subtract scattered light... this needs to be done before flatfielding.
+        if self.par['subtract_scattlight']:
+            self.subtract_scattlight()
+
         # Flat-field the data.  This propagates the flat-fielding corrections to
         # the variance.  The returned bpm is propagated to the PypeItImage
         # bitmask below.
@@ -1106,6 +1111,27 @@ def subtract_continuum(self, force=False):
         #cont = ndimage.median_filter(self.image, size=(1,101,3), mode='reflect')
         self.steps[step] = True
 
+    def subtract_scattlight(self):
+        """
+        Analyze and subtract the scattered light from the image.
+
+        This is primarily a wrapper for
+        :func:`~pypeit.spectrographs.spectrograph.scattered_light`.
+
+        """
+        step = inspect.stack()[0][3]
+        if self.steps[step]:
+            # Already pattern subtracted
+            msgs.warn("The scattered light has already been subtracted from the image!")
+            return
+
+        # Loop over the images
+        for ii in range(self.nimg):
+            binning = self.detector[0].binning
+            scatt_img = self.spectrograph.scattered_light(self.image[ii, ...], binning)
+            self.image[ii, ...] -= scatt_img
+        self.steps[step] = True
+
     def trim(self, force=False):
         """
         Trim image attributes to include only the science data.
@@ -1244,5 +1270,3 @@ def build_mosaic(self):
     def __repr__(self):
         return f'<{self.__class__.__name__}: file={self.filename}, nimg={self.nimg}, ' \
                f'steps={self.steps}>'
-
-

pypeit/par/pypeitpar.py

@@ -219,7 +219,8 @@ def __init__(self, trim=None, apply_gain=None, orient=None,
                  dark_expscale=None,
                  empirical_rn=None, shot_noise=None, noise_floor=None,
                  use_pixelflat=None, use_illumflat=None, use_specillum=None,
-                 use_pattern=None, subtract_continuum=None, spat_flexure_correct=None):
+                 use_pattern=None, subtract_scattlight=None, subtract_continuum=None,
+                 spat_flexure_correct=None):
 
         # Grab the parameter names and values from the function
         # arguments
@@ -298,6 +299,12 @@ def __init__(self, trim=None, apply_gain=None, orient=None,
                                       'be set to True to combine arcs with multiple different lamps. ' \
                                       'For all other cases, this parameter should probably be False.'
 
+        defaults['subtract_scattlight'] = False
+        dtypes['subtract_scattlight'] = bool
+        descr['subtract_scattlight'] = 'Subtract off the scattered light from an image. This parameter should only ' \
+                                       'be set to True for spectrographs that have dedicated methods to subtract ' \
+                                       'scattered light. For all other cases, this parameter should be False.'
+
         defaults['empirical_rn'] = False
         dtypes['empirical_rn'] = bool
         descr['empirical_rn'] = 'If True, use the standard deviation in the overscan region to ' \
@@ -429,8 +436,8 @@ def __init__(self, trim=None, apply_gain=None, orient=None,
     @classmethod
     def from_dict(cls, cfg):
         k = np.array([*cfg.keys()])
-        parkeys = ['trim', 'apply_gain', 'orient', 'use_biasimage', 'subtract_continuum', 'use_pattern',
-                   'use_overscan', 'overscan_method', 'overscan_par', 'use_darkimage',
+        parkeys = ['trim', 'apply_gain', 'orient', 'use_biasimage', 'subtract_continuum', 'subtract_scattlight',
+                   'use_pattern', 'use_overscan', 'overscan_method', 'overscan_par', 'use_darkimage',
                    'dark_expscale', 'spat_flexure_correct', 'use_illumflat', 'use_specillum',
                    'empirical_rn', 'shot_noise', 'noise_floor', 'use_pixelflat', 'combine',
                    'satpix', #'calib_setup_and_bit',

pypeit/spectrographs/keck_kcwi.py

@@ -886,6 +886,11 @@ def default_pypeit_par(cls):
         par['calibrations']['illumflatframe']['process']['use_pattern'] = True
         par['calibrations']['standardframe']['process']['use_pattern'] = True
         par['scienceframe']['process']['use_pattern'] = True
+        # Subtract scattered light, but only for the pixel and illum flats,
+        # as well as and science/standard star data.
+        par['calibrations']['pixelflatframe']['process']['subtract_scattlight'] = True
+        par['calibrations']['illumflatframe']['process']['subtract_scattlight'] = True
+        par['scienceframe']['process']['subtract_scattlight'] = True
 
         # Correct the illumflat for pixel-to-pixel sensitivity variations
         par['calibrations']['illumflatframe']['process']['use_pixelflat'] = True
@@ -1016,6 +1021,47 @@ def get_rawimage(self, raw_file, det):
         # Return
         return detpar, raw_img, hdu, exptime, rawdatasec_img, oscansec_img
 
+    def scattered_light(self, frame, binning):
+        """
+        Calculate a model of the scattered light of the input frame.
+
+        Parameters
+        ----------
+        frame : `numpy.ndarray`_
+            Raw 2D data frame to be used to compute the scattered light.
+        binning : str, `numpy.ndarray`_, tuple
+            Binning of the frame (e.g. '2x1' refers to a binning of 2 in the spectral
+            direction, and a binning of 1 in the spatial direction). For the supported
+            formats, refer to :func:`~pypeit.core.parse.parse_binning`.
+
+        Returns
+        -------
+        scatt_img : `numpy.ndarray`_
+            A 2D image of the scattered light determined from the input frame
+        """
+        # Determine the spatial binning. Currently, we hard-code the regions on the KCWI detector
+        # to be used for the determination of the scattered light. This is crude, and a better model
+        # is currently under development (RJC)
+        spatbin = parse.parse_binning(binning)[1]
+        ym = 4096 // (2 * spatbin)
+        y0 = ym - 180 // spatbin
+        y1 = ym + 180 // spatbin
+        # Obtain a robust (median) "spectrum" of the scattered light at the left, right, and middle of the detector
+        scattlightl = np.nanmedian(frame[:, :20//spatbin], axis=1)[:, None]
+        scattlightr = np.nanmedian(frame[:, -40//spatbin:], axis=1)[:, None]
+        scattlight = np.nanmedian(frame[:, y0:y1], axis=1)[:, None]
+        # The following algorithm is a polynomial fit to:
+        # (1) The left half of the detector (using the left and middle scattered light spectrum)
+        # (2) The right half of the detector (using the right and middle scattered light spectrum)
+        # We then ensure that the model is continuous and smooth at the middle of the detector.
+        medl = np.median(scattlightl / scattlight)
+        medr = np.median(scattlightr / scattlight)
+        spatimg = np.meshgrid(np.arange(frame.shape[1]), np.arange(frame.shape[0]))[0]
+        scatt_img = np.outer(scattlight, np.ones(frame.shape[1]))
+        scatt_img[spatimg < ym] *= 1 + (medl - 1) * (spatimg[spatimg < ym] / ym - 1) ** 2
+        scatt_img[spatimg >= ym] *= 1 + (medr - 1) * (spatimg[spatimg >= ym] / (4095//spatbin - ym) - ym / (4095//spatbin - ym)) ** 2
+        return scatt_img
+
     def fit_2d_det_response(self, det_resp, gpmask):
         r"""
         Perform a 2D model fit to the KCWI detector response.

pypeit/spectrographs/spectrograph.py

@@ -1894,6 +1894,28 @@ def calc_pattern_freq(self, frame, rawdatasec_img, oscansec_img, hdu):
         msgs.info("Pattern noise removal is not implemented for spectrograph {0:s}".format(self.name))
         return []
 
+    def scattered_light(self, frame, binning):
+        """
+        Calculate a model of the scattered light of the input frame.
+
+        Parameters
+        ----------
+        frame : `numpy.ndarray`_
+            Raw 2D data frame to be used to compute the scattered light.
+        binning : str, `numpy.ndarray`_, tuple
+            Binning of the frame (e.g. '2x1' refers to a binning of 2 in the spectral
+            direction, and a binning of 1 in the spatial direction). For the supported
+            formats, refer to :func:`~pypeit.core.parse.parse_binning`.
+
+        Returns
+        -------
+        scatt_img : `numpy.ndarray`_, float
+            A 2D image of the scattered light determined from the input frame.
+            Alternatively, if a constant value is used, a constant floating point
+            value can be returned as well.
+        """
+        msgs.info("Scattered light removal is not implemented for spectrograph {0:s}".format(self.name))
+        return 0.0
 
     def __repr__(self):
         """Return a string representation of the instance."""

sphinx.readme

@@ -1,5 +1,8 @@
 # Buiding the docs:
 
+# First, make sure that you have the PypeIt Development Suite installed
+# on your system. Then execute the following commands in a terminal:
+#
 # ./update_docs
 # cd doc
 # make html
@@ -8,7 +11,7 @@
 We use Sphinx to build our documentation:
 https://www.sphinx-doc.org/en/master/index.html
 
-Dostrings should be in rst format:
+Docstrings should be in rst format:
 http://docutils.sourceforge.net/rst.html
 
 Google format example (preferred):