KCWI scattered light

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

@@ -124,6 +124,19 @@ def parse_binning(binning:str):
 
     Args:
         binning (str, `numpy.ndarray`_, 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

pypeit/core/skysub.py

@@ -1450,7 +1450,7 @@ def convolve_skymodel(input_img, fwhm_map, thismask, subpixel=5, nsample=10):
         xkern = np.arange(kernsize, dtype=int) - midp
         kern = np.exp(-0.5 * (xkern / kernwids[kk]) ** 2)
         kern = kern / np.sum(kern)
-        conv_allkern[:, :, kk] = utils.rebinND(utils.convolve_fft(_input_img, _input_msk, kern), input_img.shape)
+        conv_allkern[:, :, kk] = utils.rebinND(utils.convolve_fft(_input_img, kern, _input_msk), input_img.shape)
 
     # Collect all of the images
     msgs.info(f"Collating all convolution steps")

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.spectrograph.procimg.subtract_pattern`.
+
+        """
+        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,12 @@ 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['calibrations']['standardframe']['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
@@ -1017,6 +1023,39 @@ 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 `~pypeit.core.parse.parse_binning`.
+
+        Returns
+        -------
+        scatt_img : `numpy.ndarray`_
+            A 2D image of the scattered light determined from the input frame
+        """
+        spatbin = parse.parse_binning(binning)[1]
+        ym = 4096 // (2 * spatbin)
+        y0 = ym - 180 // spatbin
+        y1 = ym + 180 // spatbin
+        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]
+        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

@@ -1879,6 +1879,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 `~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."""