Tweak edges using a gradient method

doc/releases/1.15.1dev.rst

@@ -19,7 +19,7 @@ that need it; see new ``filter2`` parameter.
   is generated using the A frame, while the existing `skymodel` is generated using the A-B frame.
   This allows to have an actual 1D extracted sky spectrum and to perform the flexure correction for
   `bkg_redux` reduction.
-  
+
 - Various improvements in the flexure correction and added the possibility to use a modeled archive
   sky spectrum generated with `pypeit.wavemodel.nearIR_modelsky()` to perform the flexure correction.
 
@@ -57,8 +57,8 @@ Bug Fixes
   produced error cubes that were not properly propagating the noise. The error cubes of the NGP
   algorithm were unaffected. The error cubes are now regularly inspected with vet tests to ensure
   the error cubes are reliable.
+- Allow `spec_flex_shift()` to take as input either the name of an archive sky spectrum or
+  directly a sky spectrum. This fixes a bug in the flexure correction for SlicerIFU.
 - Fix a bug (introduced in a recent PR) that was generating an error if less than 2 spec1d
   files were used with `pypeit_coadd_1dspec`. Now the script can be run with only one
   file (as it was before).
-- Allow `spec_flex_shift()` to take as input either the name of an archive sky spectrum or
-  directly a sky spectrum. This fixes a bug in the flexure correction for SlicerIFU.

pypeit/alignframe.py

@@ -181,8 +181,8 @@ def build_traces(self, show_peaks=False, debug=False):
             dict:  self.align_dict
         """
         # Generate slits
-        slitid_img_init = self.slits.slit_img(initial=True)
-        left, right, _ = self.slits.select_edges(initial=True)
+        slitid_img_init = self.slits.slit_img()
+        left, right, _ = self.slits.select_edges()
         align_prof = dict({})
 
         # Go through the slits

pypeit/calibrations.py

@@ -607,7 +607,7 @@ def get_scattlight(self):
 
         spatbin = parse.parse_binning(binning)[1]
         pad = self.par['scattlight_pad'] // spatbin
-        offslitmask = self.slits.slit_img(pad=pad, initial=True, flexure=None) == -1
+        offslitmask = self.slits.slit_img(pad=pad, flexure=None) == -1
 
         # Get starting parameters for the scattered light model
         x0, bounds = self.spectrograph.scattered_light_archive(binning, dispname)
@@ -766,6 +766,7 @@ def get_flats(self):
             msgs.info('Creating slit-illumination flat calibration frame using files: ')
             for f in raw_illum_files:
                 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, scattlight=self.msscattlight,

pypeit/coadd3d.py

@@ -838,8 +838,8 @@ def add_grating_corr(self, flatfile, waveimg, slits, spat_flexure=None):
             msgs.info("Calculating relative sensitivity for grating correction")
             # Load the Flat file
             flatimages = flatfield.FlatImages.from_file(flatfile, chk_version=self.chk_version)
-            total_illum = flatimages.fit2illumflat(slits, finecorr=False, frametype='illum', initial=True, spat_flexure=spat_flexure) * \
-                          flatimages.fit2illumflat(slits, finecorr=True, frametype='illum', initial=True, spat_flexure=spat_flexure)
+            total_illum = flatimages.fit2illumflat(slits, finecorr=False, frametype='illum', spat_flexure=spat_flexure) * \
+                          flatimages.fit2illumflat(slits, finecorr=True, frametype='illum', spat_flexure=spat_flexure)
             flatframe = flatimages.pixelflat_raw / total_illum
             if flatimages.pixelflat_spec_illum is None:
                 # Calculate the relative scale
@@ -950,7 +950,7 @@ def get_alignments(self, spec2DObj, slits, spat_flexure=None):
             msgs.info("Using slit edges for astrometric transform")
         # If nothing better was provided, use the slit edges
         if alignments is None:
-            left, right, _ = slits.select_edges(initial=True, flexure=spat_flexure)
+            left, right, _ = slits.select_edges(flexure=spat_flexure)
             locations = [0.0, 1.0]
             traces = np.append(left[:, None, :], right[:, None, :], axis=1)
         else:
@@ -1015,8 +1015,8 @@ def load(self):
             # Initialise the slit edges
             msgs.info("Constructing slit image")
             slits = spec2DObj.slits
-            slitid_img_init = slits.slit_img(pad=0, initial=True, flexure=spat_flexure)
-            slits_left, slits_right, _ = slits.select_edges(initial=True, flexure=spat_flexure)
+            slitid_img = slits.slit_img(pad=0, flexure=spat_flexure)
+            slits_left, slits_right, _ = slits.select_edges(flexure=spat_flexure)
 
             # The order of operations below proceeds as follows:
             #  (1) Get science image
@@ -1041,7 +1041,7 @@ def load(self):
             bpmmask = spec2DObj.bpmmask
 
             # Mask the edges of the spectrum where the sky model is bad
-            sky_is_good = datacube.make_good_skymask(slitid_img_init, spec2DObj.tilts)
+            sky_is_good = datacube.make_good_skymask(slitid_img, spec2DObj.tilts)
 
             # TODO :: Really need to write some detailed information in the docs about all of the various corrections that can optionally be applied
 
@@ -1074,7 +1074,7 @@ def load(self):
             if self.mnmx_wv is None:
                 self.mnmx_wv = np.zeros((len(self.spec2d), slits.nslits, 2))
             for slit_idx, slit_spat in enumerate(slits.spat_id):
-                onslit_init = (slitid_img_init == slit_spat)
+                onslit_init = (slitid_img == slit_spat)
                 self.mnmx_wv[ff, slit_idx, 0] = np.min(waveimg[onslit_init])
                 self.mnmx_wv[ff, slit_idx, 1] = np.max(waveimg[onslit_init])
 
@@ -1097,7 +1097,7 @@ def load(self):
 
             # Construct a good pixel mask
             # TODO: This should use the mask function to figure out which elements are masked.
-            onslit_gpm = (slitid_img_init > 0) & (bpmmask.mask == 0) & sky_is_good
+            onslit_gpm = (slitid_img > 0) & (bpmmask.mask == 0) & sky_is_good
 
             # Generate the alignment splines, and then retrieve images of the RA and Dec of every pixel,
             # and the number of spatial pixels in each slit
@@ -1109,7 +1109,7 @@ def load(self):
             crval_wv = self.cubepar['wave_min'] if self.cubepar['wave_min'] is not None else wave0
             cd_wv = self.cubepar['wave_delta'] if self.cubepar['wave_delta'] is not None else dwv
             self.all_wcs.append(self.spec.get_wcs(spec2DObj.head0, slits, detector.platescale, crval_wv, cd_wv))
-            ra_img, dec_img, minmax = slits.get_radec_image(self.all_wcs[ff], alignSplines, spec2DObj.tilts, initial=True, flexure=spat_flexure)
+            ra_img, dec_img, minmax = slits.get_radec_image(self.all_wcs[ff], alignSplines, spec2DObj.tilts, flexure=spat_flexure)
 
             # Extract wavelength and delta wavelength arrays from the images
             wave_ext = waveimg[onslit_gpm]
@@ -1198,7 +1198,7 @@ def load(self):
             wghts = self.weights[ff] * np.ones(sciImg.shape)
 
             # Get the slit image and then unset pixels in the slit image that are bad
-            slitid_img_gpm = slitid_img_init * onslit_gpm.astype(int)
+            slitid_img_gpm = slitid_img * onslit_gpm.astype(int)
 
             # If individual frames are to be output without aligning them,
             # there's no need to store information, just make the cubes now
@@ -1209,7 +1209,7 @@ def load(self):
                 else:
                     outfile = datacube.get_output_filename(fil, self.cubepar['output_filename'], self.combine, ff + 1)
                 # Get the coordinate bounds
-                slitlength = int(np.round(np.median(slits.get_slitlengths(initial=True, median=True))))
+                slitlength = int(np.round(np.median(slits.get_slitlengths(median=True))))
                 numwav = int((np.max(waveimg) - wave0) / dwv)
                 bins = self.spec.get_datacube_bins(slitlength, minmax, numwav)
                 # Set the wavelength range of the white light image.

pypeit/core/datacube.py

@@ -596,13 +596,16 @@ def get_whitelight_range(wavemin, wavemax, wl_range):
     return wlrng
 
 
-def make_whitelight_fromcube(cube, wave=None, wavemin=None, wavemax=None):
+def make_whitelight_fromcube(cube, bpmcube, wave=None, wavemin=None, wavemax=None):
     """
     Generate a white light image using an input cube.
 
     Args:
         cube (`numpy.ndarray`_):
             3D datacube (the final element contains the wavelength dimension)
+        bpmcube (`numpy.ndarray`_):
+            3D bad pixel mask cube (the final element contains the wavelength dimension).
+            A value of 1 indicates a bad pixel.
         wave (`numpy.ndarray`_, optional):
             1D wavelength array. Only required if wavemin or wavemax are not
             None.
@@ -619,7 +622,6 @@ def make_whitelight_fromcube(cube, wave=None, wavemin=None, wavemax=None):
         A whitelight image of the input cube (of type `numpy.ndarray`_).
     """
     # Make a wavelength cut, if requested
-    cutcube = cube.copy()
     if wavemin is not None or wavemax is not None:
         # Make some checks on the input
         if wave is None:
@@ -635,10 +637,15 @@ def make_whitelight_fromcube(cube, wave=None, wavemin=None, wavemax=None):
         ww = np.where((wave >= wavemin) & (wave <= wavemax))[0]
         wmin, wmax = ww[0], ww[-1]+1
         cutcube = cube[:, :, wmin:wmax]
+        # Cut the bad pixel mask and convert it to a good pixel mask
+        cutgpmcube = np.logical_not(bpmcube[:, :, wmin:wmax])
+    else:
+        cutcube = cube.copy()
+        cutgpmcube = np.logical_not(bpmcube)
     # Now sum along the wavelength axis
-    nrmval = np.sum(cutcube != 0.0, axis=2)
-    nrmval[nrmval == 0.0] = 1.0
-    wl_img = np.sum(cutcube, axis=2) / nrmval
+    nrmval = np.sum(cutgpmcube, axis=2)
+    nrmval[nrmval == 0] = 1.0
+    wl_img = np.sum(cutcube*cutgpmcube, axis=2) / nrmval
     return wl_img
 
 
@@ -1594,7 +1601,7 @@ def generate_cube_subpixel(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_im
             whitelight_range[0], whitelight_range[1]))
         # Get the output filename for the white light image
         out_whitelight = get_output_whitelight_filename(outfile)
-        whitelight_img = make_whitelight_fromcube(flxcube, wave=wave, wavemin=whitelight_range[0], wavemax=whitelight_range[1])
+        whitelight_img = make_whitelight_fromcube(flxcube, bpmcube, wave=wave, wavemin=whitelight_range[0], wavemax=whitelight_range[1])
         msgs.info("Saving white light image as: {0:s}".format(out_whitelight))
         img_hdu = fits.PrimaryHDU(whitelight_img.T, header=whitelight_wcs.to_header())
         img_hdu.writeto(out_whitelight, overwrite=overwrite)
@@ -1698,8 +1705,7 @@ def subpixellate(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_img_gpm, wgh
     Returns:
         :obj:`tuple`: Three or four `numpy.ndarray`_ objects containing (1) the
         datacube generated from the subpixellated inputs, (2) the corresponding
-        variance cube, (3) the corresponding bad pixel mask cube, and (4) the
-        residual cube.  The latter is only returned if debug is True.
+        variance cube, and (3) the corresponding bad pixel mask cube.
     """
     # Check the inputs for combinations of lists or not
     _sciImg, _ivarImg, _waveImg, _gpmImg, _wghtImg, _all_wcs, _tilts, _slits, _astrom_trans, _all_dar, _ra_offset, _dec_offset = \

pypeit/core/flat.py

@@ -12,48 +12,13 @@
 import scipy.interpolate
 import scipy.ndimage
 import matplotlib.pyplot as plt
+from astropy import convolution
 
 from IPython import embed
 
 from pypeit import msgs
-from pypeit.core import parse
-from pypeit.core import pixels
-from pypeit.core import tracewave
 from pypeit.core import coadd
 from pypeit import utils
-from pypeit.core import pydl
-
-# TODO: Put this in utils
-def linear_interpolate(x1, y1, x2, y2, x):
-    r"""
-    Interplate or extrapolate between two points.
-
-    Given a line defined two points, :math:`(x_1,y_1)` and
-    :math:`(x_2,y_2)`, return the :math:`y` value of a new point on
-    the line at coordinate :math:`x`.
-
-    This function is meant for speed. No type checking is performed and
-    the only check is that the two provided ordinate coordinates are not
-    numerically identical. By definition, the function will extrapolate
-    without any warning.
-
-    Args:
-        x1 (:obj:`float`):
-            First abscissa position
-        y1 (:obj:`float`):
-            First ordinate position
-        x2 (:obj:`float`):
-            Second abscissa position
-        y3 (:obj:`float`):
-            Second ordinate position
-        x (:obj:`float`):
-            Abcissa for new value
-
-    Returns:
-        :obj:`float`: Interpolated/extrapolated value of ordinate at
-        :math:`x`.
-    """
-    return y1 if np.isclose(x1,x2) else y1 + (x-x1)*(y2-y1)/(x2-x1)
 
 
 # TODO: Make this function more general and put it in utils.
@@ -505,10 +470,98 @@ def poly_map(rawimg, rawivar, waveimg, slitmask, slitmask_trim, modelimg, deg=3,
     return modelmap, relscale
 
 
+def tweak_slit_edges_gradient(left, right, spat_coo, norm_flat, debug=False):
+    r""" Adjust slit edges based on the gradient of the normalized
+    flat-field illumination profile.
+
+    Args:
+        left (`numpy.ndarray`_):
+            Array with the left slit edge for a single slit. Shape is
+            :math:`(N_{\rm spec},)`.
+        right (`numpy.ndarray`_):
+            Array with the right slit edge for a single slit. Shape
+            is :math:`(N_{\rm spec},)`.
+        spat_coo (`numpy.ndarray`_):
+            Spatial pixel coordinates in fractions of the slit width
+            at each spectral row for the provided normalized flat
+            data. Coordinates are relative to the left edge (with the
+            left edge at 0.). Shape is :math:`(N_{\rm flat},)`.
+            Function assumes the coordinate array is sorted.
+        norm_flat (`numpy.ndarray`_)
+            Normalized flat data that provide the slit illumination
+            profile. Shape is :math:`(N_{\rm flat},)`.
+        debug (:obj:`bool`, optional):
+            If True, the function will output plots to test if the
+            fitting is working correctly.
+
+    Returns:
+        tuple: Returns six objects:
+
+            - The threshold used to set the left edge
+            - The fraction of the slit that the left edge is shifted to
+              the right
+            - The adjusted left edge
+            - The threshold used to set the right edge
+            - The fraction of the slit that the right edge is shifted to
+              the left
+            - The adjusted right edge
+    """
+    # Check input
+    nspec = len(left)
+    if len(right) != nspec:
+        msgs.error('Input left and right traces must have the same length!')
+
+    # Median slit width
+    slitwidth = np.median(right - left)
+
+    # Calculate the gradient of the normalized flat profile
+    grad_norm_flat = np.gradient(norm_flat)
+    # Smooth with a Gaussian kernel
+    sig_res = norm_flat.size / slitwidth
+    gauss_kernel = convolution.Gaussian1DKernel(sig_res)
+    grad_norm_flat_smooth = convolution.convolve(grad_norm_flat, gauss_kernel, boundary='extend')
+
+    # Find the location of the minimum/maximum gradient - this is the amount of shift required
+    left_shift = spat_coo[np.argmax(grad_norm_flat_smooth)]
+    right_shift = spat_coo[np.argmin(grad_norm_flat_smooth)]-1.0
+    msgs.info('Tweaking left slit boundary by {0:.1f}%'.format(100 * left_shift) +
+              ' ({0:.2f} pixels)'.format(left_shift * slitwidth))
+    msgs.info('Tweaking right slit boundary by {0:.1f}%'.format(100 * right_shift) +
+                ' ({0:.2f} pixels)'.format(right_shift * slitwidth))
+
+    # Calculate the tweak for the left edge
+    new_left = np.copy(left) + left_shift * slitwidth
+    new_right = np.copy(right) + right_shift * slitwidth
+
+    # Calculate the value of the threshold at the new slit edges
+    left_thresh = np.interp(left_shift, spat_coo, norm_flat)
+    right_thresh = np.interp(1+right_shift, spat_coo, norm_flat)
+
+    if debug:
+        plt.subplot(211)
+        plt.plot(spat_coo, norm_flat, 'k-')
+        plt.axvline(0.0, color='b', linestyle='-', label='initial')
+        plt.axvline(1.0, color='b', linestyle='-')
+        plt.axvline(left_shift, color='g', linestyle='-', label='tweak (gradient)')
+        plt.axvline(1+right_shift, color='g', linestyle='-')
+        plt.axhline(left_thresh, xmax=0.5, color='lightgreen', linewidth=3.0, zorder=10)
+        plt.axhline(right_thresh, xmin=0.5, color='lightgreen', linewidth=3.0, zorder=10)
+        plt.legend()
+        plt.subplot(212)
+        plt.plot(spat_coo, grad_norm_flat, 'k-')
+        plt.plot(spat_coo, grad_norm_flat_smooth, 'm-')
+        plt.axvline(0.0, color='b', linestyle='-')
+        plt.axvline(1.0, color='b', linestyle='-')
+        plt.axvline(left_shift, color='g', linestyle='-')
+        plt.axvline(1+right_shift, color='g', linestyle='-')
+        plt.show()
+    return left_thresh, left_shift, new_left, right_thresh, right_shift, new_right
+
+
 # TODO: See pypeit/deprecated/flat.py for the previous version. We need
 # to continue to vet this algorithm to make sure there are no
 # unforeseen corner cases that cause errors.
-def tweak_slit_edges(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1, debug=False):
+def tweak_slit_edges_threshold(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1, debug=False):
     r"""
     Adjust slit edges based on the normalized slit illumination profile.
 
@@ -614,10 +667,10 @@ def tweak_slit_edges(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1,
                         100*maxfrac))
             left_shift = maxfrac
         else:
-            left_shift = linear_interpolate(norm_flat[i], spat_coo[i], norm_flat[i+1],
-                                           spat_coo[i+1], left_thresh)
+            left_shift = utils.linear_interpolate(norm_flat[i], spat_coo[i], norm_flat[i+1],
+                                                  spat_coo[i+1], left_thresh)
         msgs.info('Tweaking left slit boundary by {0:.1f}%'.format(100*left_shift) +
-                  ' % ({0:.2f} pixels)'.format(left_shift*slitwidth))
+                  ' ({0:.2f} pixels)'.format(left_shift*slitwidth))
         new_left += left_shift * slitwidth
 
     # ------------------------------------------------------------------
@@ -668,15 +721,15 @@ def tweak_slit_edges(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1,
                         100*maxfrac))
             right_shift = maxfrac
         else:
-            right_shift = 1-linear_interpolate(norm_flat[i-1], spat_coo[i-1], norm_flat[i],
-                                               spat_coo[i], right_thresh)
+            right_shift = 1-utils.linear_interpolate(norm_flat[i-1], spat_coo[i-1], norm_flat[i],
+                                                     spat_coo[i], right_thresh)
         msgs.info('Tweaking right slit boundary by {0:.1f}%'.format(100*right_shift) +
-                  ' % ({0:.2f} pixels)'.format(right_shift*slitwidth))
+                  ' ({0:.2f} pixels)'.format(right_shift*slitwidth))
         new_right -= right_shift * slitwidth
 
     return left_thresh, left_shift, new_left, right_thresh, right_shift, new_right
 
-#def flatfield(sciframe, flatframe, bpm=None, illum_flat=None, snframe=None, varframe=None):
+
 def flatfield(sciframe, flatframe, varframe=None):
     r"""
     Field flatten the input image.