Merge branch 'develop' into spatflex_improve

doc/releases/1.17.1dev.rst

@@ -11,10 +11,12 @@ Dependency Changes
 ------------------
 
 
-
 Functionality/Performance Improvements and Additions
 ----------------------------------------------------
 
+- The WCS for datacubes now adopts the convention of North
+  is up and East is left. In previous version of PypeIt,
+  East was right.
 - Improved the code to compute the spatial flexure correction and added a QA
   plot to inspect the results.  Two new parset parameters, `spat_flexure_sigdetect`
   and `spat_flexure_vrange`, were added to control the detection of the slit edges
@@ -25,7 +27,6 @@ Instrument-specific Updates
 ---------------------------
 
 
-
 Script Changes
 --------------
 
@@ -40,7 +41,11 @@ Under-the-hood Improvements
 ---------------------------
 
 
-
 Bug Fixes
 ---------
 
+- Fix the WCS for the datacube generation. There was an offset
+  in both spatial dimensions equal to half the field of view.
+- Fix the code for the extraction of the 1D flat spectrum, so that
+  the spectrum is extracted even when `pixelflat_model` does not exist.
+

doc/whatsnew.rst

@@ -11,6 +11,10 @@ What's New in PypeIt
 
 ----
 
+.. include:: releases/1.17.1dev.rst
+
+----
+
 .. include:: releases/1.17.0.rst
 
 ----

presentations/py/users.py

@@ -21,9 +21,9 @@ def set_fontsize(ax, fsz):
                  ax.get_xticklabels() + ax.get_yticklabels()):
         item.set_fontsize(fsz)
 
-user_dates = ["2021-03-11", "2022-04-29", "2022-11-07", "2022-12-06", "2023-06-08", "2023-06-29", "2023-07-11", "2023-09-03", "2023-10-13", "2023-12-01", "2023-12-15", "2024-02-22", "2024-03-21", "2024-04-09", "2024-05-02", "2024-05-19", "2024-06-06", "2024-06-10", "2024-08-20"]
+user_dates = ["2021-03-11", "2022-04-29", "2022-11-07", "2022-12-06", "2023-06-08", "2023-06-29", "2023-07-11", "2023-09-03", "2023-10-13", "2023-12-01", "2023-12-15", "2024-02-22", "2024-03-21", "2024-04-09", "2024-05-02", "2024-05-19", "2024-06-06", "2024-06-10", "2024-08-20", "2024-11-15"]
 user_dates = numpy.array([numpy.datetime64(date) for date in user_dates])
-user_number = numpy.array([125, 293, 390, 394, 477, 487, 506, 518, 531, 544, 551, 568, 579, 588, 596, 603, 616, 620, 643])
+user_number = numpy.array([125, 293, 390, 394, 477, 487, 506, 518, 531, 544, 551, 568, 579, 588, 596, 603, 616, 620, 643, 688])
 
 user_pred_dates = numpy.array([numpy.datetime64(date)
                             for date in ["2024-06-10", "2024-12-31", "2025-12-31", "2026-12-31",

pypeit/alignframe.py

@@ -397,7 +397,7 @@ def build_splines(self):
             xcoord = np.arange(np.floor(np.min(xlr)), np.ceil(np.max(xlr))+1, 1.0)
             out_transform = np.zeros((self.nspec, xcoord.size))
             for sp in range(self.nspec):
-                out_transform[sp,:] = (self.spl_loc[sl][sp](xcoord) - 0.5) * self.spl_slen[sl](sp)
+                out_transform[sp,:] = self.spl_loc[sl][sp](xcoord) * self.spl_slen[sl](sp)
             self.spl_transform[sl] = RegularGridInterpolator((ycoord, xcoord), out_transform, method='linear',
                                                              bounds_error=False, fill_value=None) # This will extrapolate
             # TODO :: Remove these notes...

pypeit/core/datacube.py

@@ -979,7 +979,7 @@ def create_wcs(raImg, decImg, waveImg, slitid_img_gpm, dspat, dwave,
 
     # Generate a master WCS to register all frames
     coord_min = [_ra_min, _dec_min, _wave_min]
-    coord_dlt = [dspat, dspat, dwave]
+    coord_dlt = [-dspat, dspat, dwave]
 
     # If a reference image is being used and a white light image is requested (collapse=True) update the celestial parts
     reference_image = None
@@ -991,7 +991,7 @@ def create_wcs(raImg, decImg, waveImg, slitid_img_gpm, dspat, dwave,
         coord_dlt[:2] = imgwcs.wcs.cdelt
         numra, numdec = reference_image.shape
 
-    cubewcs = generate_WCS(coord_min, coord_dlt, equinox=equinox, name=specname)
+    cubewcs = generate_WCS(coord_min, coord_dlt, numra, equinox=equinox, name=specname)
     msgs.info(msgs.newline() + "-" * 40 +
               msgs.newline() + "Parameters of the WCS:" +
               msgs.newline() + "RA   min = {0:f}".format(coord_min[0]) +
@@ -1009,7 +1009,7 @@ def create_wcs(raImg, decImg, waveImg, slitid_img_gpm, dspat, dwave,
     return cubewcs, voxedges, reference_image
 
 
-def generate_WCS(crval, cdelt, equinox=2000.0, name="PYP_SPEC"):
+def generate_WCS(crval, cdelt, numra, equinox=2000.0, name="PYP_SPEC"):
     """
     Generate a WCS that will cover all input spec2D files
 
@@ -1020,6 +1020,10 @@ def generate_WCS(crval, cdelt, equinox=2000.0, name="PYP_SPEC"):
         cdelt (list):
             3 element list containing the delta values of the [RA,
             DEC, WAVELENGTH]
+        numra (int):
+            Number of RA values in the WCS. This is used to ensure
+            that the convention of the WCS is so that North is up
+            and East is to the left.
         equinox (float, optional):
             Equinox of the WCS
 
@@ -1037,7 +1041,7 @@ def generate_WCS(crval, cdelt, equinox=2000.0, name="PYP_SPEC"):
     w.wcs.cunit = [units.degree, units.degree, units.Angstrom]
     w.wcs.ctype = ["RA---TAN", "DEC--TAN", "WAVE"]
     w.wcs.crval = crval  # RA, DEC, and wavelength zeropoints
-    w.wcs.crpix = [0, 0, 0]  # RA, DEC, and wavelength reference pixels
+    w.wcs.crpix = [numra, 0, 0]  # RA, DEC, and wavelength reference pixels
     #w.wcs.cd = np.array([[cdval[0], 0.0, 0.0], [0.0, cdval[1], 0.0], [0.0, 0.0, cdval[2]]])
     w.wcs.cdelt = cdelt
     w.wcs.lonpole = 180.0  # Native longitude of the Celestial pole
@@ -1309,18 +1313,20 @@ def compute_weights(raImg, decImg, waveImg, sciImg, ivarImg, slitidImg,
     #idx_max = np.unravel_index(np.argmax(whitelight_img), whitelight_img.shape)
     msgs.info("Highest S/N object located at spaxel (x, y) = {0:d}, {1:d}".format(idx_max[0], idx_max[1]))
 
-    # Generate a 2D WCS to register all frames
-    coord_min = [_ra_min, _dec_min, _wave_min]
-    coord_dlt = [dspat, dspat, dwv]
-    whitelightWCS = generate_WCS(coord_min, coord_dlt)
-    wcs_scale = (1.0 * whitelightWCS.spectral.wcs.cunit[0]).to(units.Angstrom).value  # Ensures the WCS is in Angstroms
     # Make the bin edges to be at +/- 1 pixels around the maximum (i.e. summing 9 pixels total)
     numwav = int((_wave_max - _wave_min) / dwv)
     xbins = np.array([idx_max[0]-1, idx_max[0]+2]) - 0.5
     ybins = np.array([idx_max[1]-1, idx_max[1]+2]) - 0.5
     spec_bins = np.arange(1 + numwav) - 0.5
     bins = (xbins, ybins, spec_bins)
 
+    # Generate a 2D WCS to register all frames
+    numra = xbins.size - 1
+    coord_min = [_ra_min, _dec_min, _wave_min]
+    coord_dlt = [-dspat, dspat, dwv]
+    whitelightWCS = generate_WCS(coord_min, coord_dlt, numra)
+    wcs_scale = (1.0 * whitelightWCS.spectral.wcs.cunit[0]).to(units.Angstrom).value  # Ensures the WCS is in Angstroms
+
     # Extract the spectrum of the highest S/N object
     flux_stack = np.zeros((numwav, numframes))
     ivar_stack = np.zeros((numwav, numframes))

pypeit/extraction.py

@@ -15,8 +15,7 @@
 
 from pypeit import msgs, utils
 from pypeit.display import display
-from pypeit.core import skysub, extract, flexure
-
+from pypeit.core import skysub, extract, flexure, flat
 from IPython import embed
 
 
@@ -64,7 +63,7 @@ class Extract:
     @classmethod
     def get_instance(cls, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_sky=None,
                      bkg_redux_global_sky=None, waveTilts=None, tilts=None, wv_calib=None, waveimg=None,
-                     flatimg=None, bkg_redux=False, return_negative=False, std_redux=False, show=False, basename=None):
+                     flatimages=None, bkg_redux=False, return_negative=False, std_redux=False, show=False, basename=None):
         """
         Instantiate the Extract subclass appropriate for the provided
         spectrograph.
@@ -101,9 +100,9 @@ def get_instance(cls, sciImg, slits, sobjs_obj, spectrograph, par, objtype, glob
                 This is the waveCalib object which is optional, but either wv_calib or waveimg must be provided.
             waveimg (`numpy.ndarray`_, optional):
                 Wave image. Either a wave image or wv_calib object (above) must be provided
-            flatimg (`numpy.ndarray`_, optional):
-                Flat image. This is optional, but if provided, it is used to extract the
-                normalized blaze profile. Same shape as ``sciImg``.
+            flatimages (:class:`~pypeit.flatfield.FlatImages`, optional):
+                FlatImages class. This is optional, but if provided, it is used to extract the
+                normalized blaze profile.
             bkg_redux (:obj:`bool`, optional):
                 If True, the sciImg has been subtracted by
                 a background image (e.g. standard treatment in the IR)
@@ -131,12 +130,12 @@ def get_instance(cls, sciImg, slits, sobjs_obj, spectrograph, par, objtype, glob
                     if c.__name__ == (spectrograph.pypeline + 'Extract'))(
             sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_sky=global_sky,
             bkg_redux_global_sky=bkg_redux_global_sky, waveTilts=waveTilts, tilts=tilts,
-            wv_calib=wv_calib, waveimg=waveimg, flatimg=flatimg, bkg_redux=bkg_redux, return_negative=return_negative,
-            std_redux=std_redux, show=show, basename=basename)
+            wv_calib=wv_calib, waveimg=waveimg, flatimages=flatimages, bkg_redux=bkg_redux,
+            return_negative=return_negative, std_redux=std_redux, show=show, basename=basename)
 
     def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_sky=None,
                  bkg_redux_global_sky=None, waveTilts=None, tilts=None, wv_calib=None, waveimg=None,
-                 flatimg=None, bkg_redux=False, return_negative=False, std_redux=False, show=False,
+                 flatimages=None, bkg_redux=False, return_negative=False, std_redux=False, show=False,
                  basename=None):
 
         # Setup the parameters sets for this object. NOTE: This uses objtype, not frametype!
@@ -149,7 +148,6 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_
         self.par = par
         self.global_sky = global_sky if global_sky is not None else np.zeros_like(sciImg.image)
         self.bkg_redux_global_sky = bkg_redux_global_sky
-        self.flatimg = flatimg
 
         self.basename = basename
         # Parse
@@ -247,6 +245,18 @@ def __init__(self, sciImg, slits, sobjs_obj, spectrograph, par, objtype, global_
         else:
             msgs.warn("Spectral FWHM image could not be generated")
 
+        # get flatfield image for blaze function
+        self.flatimg = None
+        if flatimages is not None:
+            # This way of getting the flat image (instead of just reading flatimages.pixelflat_model) ensures that the
+            # flat image is available also when an archival pixel flat is used.
+            flat_raw = flatimages.pixelflat_raw if flatimages.pixelflat_raw is not None else flatimages.illumflat_raw
+            if flat_raw is not None and flatimages.pixelflat_norm is not None:
+                # TODO: Can we just use flat_raw if flatimages.pixelflat_norm is None?
+                self.flatimg, _ = flat.flatfield(flat_raw, flatimages.pixelflat_norm)
+        if self.flatimg is None:
+            msgs.warn("No flat image was found. A spectrum of the flatfield will not be extracted!")
+
         # Now apply a global flexure correction to each slit provided it's not a standard star
         if self.par['flexure']['spec_method'] != 'skip' and not self.std_redux:
             # Update slitshift values

pypeit/pypeit.py

@@ -1029,18 +1029,14 @@ def extract_one(self, frames, det, sciImg, bkg_redux_sciimg, objFind, initial_sk
 
         if not self.par['reduce']['extraction']['skip_extraction']:
             msgs.info(f"Extraction begins for {self.basename} on det={det}")
-            # set the flatimg, if it exists
-            flatimg = None if self.caliBrate.flatimages is None else self.caliBrate.flatimages.pixelflat_model
-            if flatimg is None:
-                msgs.warn("No flat image was found. A spectrum of the flatfield will not be extracted!")
             # Instantiate Reduce object
             # Required for pipeline specific object
             # At instantiation, the fullmask in self.sciImg is modified
             # TODO Are we repeating steps in the init for FindObjects and Extract??
             self.exTract = extraction.Extract.get_instance(
                 sciImg, slits, sobjs_obj, self.spectrograph,
                 self.par, self.objtype, global_sky=final_global_sky, bkg_redux_global_sky=bkg_redux_global_sky,
-                waveTilts=self.caliBrate.wavetilts, wv_calib=self.caliBrate.wv_calib, flatimg=flatimg,
+                waveTilts=self.caliBrate.wavetilts, wv_calib=self.caliBrate.wv_calib, flatimages=self.caliBrate.flatimages,
                 bkg_redux=self.bkg_redux, return_negative=self.par['reduce']['extraction']['return_negative'],
                 std_redux=self.std_redux, basename=self.basename, show=self.show)
             # Perform the extraction

pypeit/slittrace.py

@@ -556,7 +556,7 @@ def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_off
             minmax[slit_idx, 0] = np.min(evalpos)
             minmax[slit_idx, 1] = np.max(evalpos)
             # Calculate the WCS from the pixel positions
-            slitID = np.ones(evalpos.size) * slit_idx + slice_offset - wcs.wcs.crpix[0]
+            slitID = np.ones(evalpos.size) * slit_idx + slice_offset
             world_ra, world_dec, _ = wcs.wcs_pix2world(slitID, evalpos, tilts[onslit_init]*(self.nspec-1), 0)
             # Set the RA first and DEC next
             raimg[onslit] = world_ra.copy()

pypeit/specobjs.py

@@ -237,8 +237,8 @@ def unpack_object(self, ret_flam=False, log10blaze=False, min_blaze_value=1e-3,
 
         # Check for missing data
         none_flux = [f is None for f in getattr(self, flux_key)]
+        other = 'OPT' if extract_type == 'BOX' else 'BOX'
         if np.any(none_flux):
-            other = 'OPT' if extract_type == 'BOX' else 'BOX'
             msg = f"{extract_type} extracted flux is not available for all slits/orders. " \
                   f"Consider trying the {other} extraction."
             if not remove_missing:
@@ -249,6 +249,12 @@ def unpack_object(self, ret_flam=False, log10blaze=False, min_blaze_value=1e-3,
                 # Remove missing data
                 r_indx = np.where(none_flux)[0]
                 self.remove_sobj(r_indx)
+        # check for missing blaze
+        if extract_blaze:
+            none_blaze = [f is None for f in getattr(self, blaze_key)]
+            if np.any(none_blaze):
+                msgs.error(f"{extract_type} extracted blaze is not available for all slits/orders. "
+                           f"Consider trying the {other} extraction, or NOT using the flat.")
 
         #
         norddet = self.nobj
@@ -299,7 +305,7 @@ def unpack_object(self, ret_flam=False, log10blaze=False, min_blaze_value=1e-3,
         # Return
         if self[0].PYPELINE in ['MultiSlit', 'SlicerIFU'] and self.nobj == 1:
             meta_spec['ECH_ORDERS'] = None
-            blaze_ret = blaze_function.reshape(nspec) if extract_blaze else None
+            blaze_ret = blaze_function.reshape(nspec) if blaze_function is not None else None
             return wave.reshape(nspec), flux.reshape(nspec), flux_ivar.reshape(nspec), \
                    flux_gpm.reshape(nspec), blaze_ret, meta_spec, self.header
         else:

pypeit/spectrographs/keck_kcwi.py

@@ -720,7 +720,7 @@ def get_datacube_bins(self, slitlength, minmax, num_wave):
             when constructing a histogram of the spec2d files. The elements
             are :math:`(x,y,\lambda)`.
         """
-        xbins = np.arange(1 + 24) - 24/2 - 0.5
+        xbins = np.arange(1 + 24) - 0.5
         ybins = np.linspace(np.min(minmax[:, 0]), np.max(minmax[:, 1]), 1+slitlength) - 0.5
         spec_bins = np.arange(1+num_wave) - 0.5
         return xbins, ybins, spec_bins

pypeit/wavetilts.py

@@ -797,6 +797,13 @@ def run(self, doqa=True, debug=False, show=False):
                                        self.spat_order[slit_idx], self.spec_order[slit_idx],
                                        slit_idx,
                                        doqa=doqa, show_QA=show)
+            # Flag slits with high number of rejected pixels (>95%)
+            # TODO: Is 95% the right threshold?
+            _gpm = self.all_fit_dict[slit_idx]['pypeitFit'].bool_gpm
+            if np.sum(np.logical_not(_gpm)) > 0.95 * _gpm.size:
+                msgs.warn(f'Large number of pixels rejected in the fit. This slit/order will not be reduced!')
+                self.slits.mask[slit_idx] = self.slits.bitmask.turn_on(self.slits.mask[slit_idx], 'BADTILTCALIB')
+                continue
             self.coeffs[:self.spec_order[slit_idx]+1,:self.spat_order[slit_idx]+1,slit_idx] = coeff_out
 
             # TODO: Need a way to assess the success of fit_tilts and
@@ -807,6 +814,12 @@ def run(self, doqa=True, debug=False, show=False):
             # images, trace images, and pixelflats etc.
             self.tilts = tracewave.fit2tilts(self.slitmask_science.shape, coeff_out,
                                              self.par['func2d'])
+            # Check that the tilts image has values that span a reasonable range
+            # TODO: Is this the right threshold?
+            if np.nanmax(self.tilts) - np.nanmin(self.tilts) < 0.8:
+                msgs.warn('Tilts image fit not good. This slit/order will not be reduced!')
+                self.slits.mask[slit_idx] = self.slits.bitmask.turn_on(self.slits.mask[slit_idx], 'BADTILTCALIB')
+                continue
             # Save to final image
             thismask_science = self.slitmask_science == slit_spat
             self.final_tilts[thismask_science] = self.tilts[thismask_science]