Hotfix: Datacube WCS fix + adopt East-West convention

doc/releases/1.17.1dev.rst

@@ -11,16 +11,17 @@ 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.
 
 Instrument-specific Updates
 ---------------------------
 
 
-
 Script Changes
 --------------
 
@@ -35,6 +36,9 @@ 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.
+

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/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/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