Enable resample_spec for NIRSpec line lamp exposures

CHANGES.rst

@@ -12,6 +12,7 @@ assign_wcs
 ----------
 
 - Add nrs_verify to the NIRSpec exposure list [#5403]
+- Enable resample_spec for NIRSpec line lamp exposures [#5484]
 
 associations
 ------------
@@ -113,6 +114,7 @@ pipeline
 
 - Update ``Image3Pipeline`` to allow sky subtraction when input contains
   only one image (group). [#5423]
+- Enable resample_spec for NIRSpec line lamp exposures in Spec2Pipeline [#5484]
 
 ramp_fitting
 ------------
@@ -129,6 +131,7 @@ resample
 - Implement memory check in resample to prevent huge arrays [#5354]
 
 - Add ``pixel_scale_ratio`` parameter to allow finer output grid. [#5389]
+- Enable resample_spec for NIRSpec line lamp exposures [#5484]
 
 saturation
 ----------

jwst/assign_wcs/nirspec.py

@@ -983,7 +983,8 @@ def gwa_to_slit(open_slits, input_model, disperser,
     # The wavelength units up to this point are
     # meters as required by the pipeline but the desired output wavelength units is microns.
     # So we are going to Scale the spectral units by 1e6 (meters -> microns)
-    if input_model.meta.instrument.filter == 'OPAQUE':
+    is_lamp_exposure = input_model.meta.exposure.type in ['NRS_LAMP', 'NRS_AUTOWAVE', 'NRS_AUTOFLAT']
+    if input_model.meta.instrument.filter == 'OPAQUE' or is_lamp_exposure:
         lgreq = lgreq | Scale(1e6)
 
     msa = MSAModel(reference_files['msa'])

jwst/lib/exposure_types.py

@@ -41,6 +41,12 @@ def is_nrs_flatlamp(datamodel):
     lamp_state = datamodel.meta.instrment.lamp_state.lower()
     return lamp_state[0:4] == 'flat'
 
+def is_nrs_slit_linelamp(datamodel):
+    lamp_mode = datamodel.meta.instrument.lamp_mode.lower()
+    exp_type = datamodel.meta.exposure.type.lower()
+    return lamp_mode in ['msaspec', 'fixedslit', 'brightobj'] and is_nrs_linelamp(datamodel) and \
+           exp_type in ['nrs_autowave', 'nrs_lamp']
+
 def is_nrs_ifu_linelamp(datamodel):
     lamp_mode = datamodel.meta.instrument.lamp_mode.lower()
     return lamp_mode == 'ifu' and is_nrs_linelamp(datamodel)

jwst/pipeline/calwebb_spec2.py

@@ -5,7 +5,7 @@
 
 from .. import datamodels
 from ..assign_wcs.util import NoDataOnDetectorError
-from ..lib.exposure_types import is_nrs_ifu_flatlamp, is_nrs_ifu_linelamp
+from ..lib.exposure_types import is_nrs_ifu_flatlamp, is_nrs_ifu_linelamp, is_nrs_slit_linelamp
 from ..stpipe import Pipeline
 
 # step imports
@@ -254,6 +254,11 @@ def process_exposure_product(
             # Call the resample_spec step for 2D slit data
             resampled = self.resample_spec(calibrated)
 
+        elif is_nrs_slit_linelamp(calibrated):
+
+            # Call resample_spec for NRS 2D line lamp slit data
+            resampled = self.resample_spec(calibrated)
+
         elif (exp_type in ['MIR_MRS', 'NRS_IFU']) or is_nrs_ifu_linelamp(calibrated):
 
             # Call the cube_build step for IFU data;

jwst/resample/gwcs_drizzle.py

@@ -9,7 +9,6 @@
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
 
-
 class GWCSDrizzle:
     """
     Combine images using the drizzle algorithm
@@ -420,9 +419,10 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon,
     # print("x slice: ", pixmap[y_mid,:,0])
     # print("y slice: ", pixmap[:,x_mid,1])
     # print("insci: ", insci)
-
     # Call 'drizzle' to perform image combination
+
     log.info(f"Drizzling {insci.shape} --> {outsci.shape}")
+
     _vers, nmiss, nskip = cdrizzle.tdriz(
         insci, inwht, pixmap,
         outsci, outwht, outcon,
@@ -436,5 +436,4 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon,
         wtscale=wt_scl,
         fillstr=fillval
         )
-
     return _vers, nmiss, nskip

jwst/resample/resample_spec.py

@@ -157,15 +157,20 @@ def build_interpolated_output_wcs(self, refmodel=None):
                 ra_center_pt = np.nanmean(ra_center)
                 dec_center_pt = np.nanmean(dec_center)
 
-                # ra and dec this converted to tangent projection
-                tan = Pix2Sky_TAN()
-                native2celestial = RotateNative2Celestial(ra_center_pt, dec_center_pt, 180)
-                undist2sky1 = tan | native2celestial
-                # Filter out RuntimeWarnings due to computed NaNs in the WCS
-                warnings.simplefilter("ignore")
-                # at this center of slit find x,y tangent proction - x_tan, y_tan
-                x_tan, y_tan = undist2sky1.inverse(ra, dec)
-                warnings.resetwarnings()
+                if resample_utils.is_sky_like(model.meta.wcs.output_frame):
+                    # convert ra and dec to tangent projection
+                    tan = Pix2Sky_TAN()
+                    native2celestial = RotateNative2Celestial(ra_center_pt, dec_center_pt, 180)
+                    undist2sky1 = tan | native2celestial
+                    # Filter out RuntimeWarnings due to computed NaNs in the WCS
+                    warnings.simplefilter("ignore")
+                    # at this center of slit find x,y tangent projection - x_tan, y_tan
+                    x_tan, y_tan = undist2sky1.inverse(ra, dec)
+                    warnings.resetwarnings()
+                else:
+                    # for non sky-like output frames, no need to do tangent plane projections
+                    # but we still use the same variables
+                    x_tan, y_tan = ra, dec
 
                 # pull out data from center
                 if spectral_axis == 0:
@@ -286,11 +291,13 @@ def build_interpolated_output_wcs(self, refmodel=None):
             ra_center_final = ra_center_pt
             dec_center_final = dec_center_pt
 
-        native2celestial = RotateNative2Celestial(ra_center_final, dec_center_final, 180)
-        undist2sky = tan | native2celestial
-
-        # find the spatial size of the output - same in x,y
-        x_tan_all, y_tan_all = undist2sky.inverse(all_ra, all_dec)
+        if resample_utils.is_sky_like(model.meta.wcs.output_frame):
+            native2celestial = RotateNative2Celestial(ra_center_final, dec_center_final, 180)
+            undist2sky = tan | native2celestial
+            # find the spatial size of the output - same in x,y
+            x_tan_all, _ = undist2sky.inverse(all_ra, all_dec)
+        else:
+            x_tan_all, _ = all_ra, all_dec
         x_min = np.amin(x_tan_all)
         x_max = np.amax(x_tan_all)
         x_size = int(np.ceil((x_max - x_min)/np.absolute(pix_to_xtan.slope)))
@@ -301,11 +308,17 @@ def build_interpolated_output_wcs(self, refmodel=None):
             x_size = len(x_tan_array)
 
         # define the output wcs
-        transform = mapping | (pix_to_xtan & pix_to_ytan | undist2sky) & pix_to_wavelength
+        if resample_utils.is_sky_like(model.meta.wcs.output_frame):
+            transform = mapping | (pix_to_xtan & pix_to_ytan | undist2sky) & pix_to_wavelength
+        else:
+            transform = mapping | (pix_to_xtan & pix_to_ytan) & pix_to_wavelength
 
         det = cf.Frame2D(name='detector', axes_order=(0, 1))
-        sky = cf.CelestialFrame(name='sky', axes_order=(0, 1),
-                                reference_frame=coord.ICRS())
+        if resample_utils.is_sky_like(model.meta.wcs.output_frame):
+            sky = cf.CelestialFrame(name='sky', axes_order=(0, 1),
+                                    reference_frame=coord.ICRS())
+        else:
+            sky = cf.Frame2D(name=f'resampled_{model.meta.wcs.output_frame.name}', axes_order=(0, 1))
         spec = cf.SpectralFrame(name='spectral', axes_order=(2,),
                                 unit=(u.micron,), axes_names=('wavelength',))
         world = cf.CompositeFrame([sky, spec], name='world')

jwst/resample/resample_utils.py

@@ -4,6 +4,8 @@
 import numpy as np
 from astropy import wcs as fitswcs
 from astropy.modeling import Model
+from astropy.modeling.models import Mapping
+from astropy import units as u
 from gwcs import WCS, wcstools
 
 from jwst.assign_wcs.util import wcs_from_footprints, wcs_bbox_from_shape
@@ -108,7 +110,13 @@ def reproject(wcs1, wcs2):
     if isinstance(wcs2, fitswcs.WCS):
         backward_transform = wcs2.all_world2pix
     elif isinstance(wcs2, WCS):
-        backward_transform = wcs2.backward_transform
+        if not is_sky_like(wcs1.output_frame):
+            # nirspec lamps: simplify backward transformation by omitting the msa_x (it's constant)
+            # and just using the wavelength lookup table [1] and linear msa_y transformation [2]
+            log.info("Custom transform for NRS Lamp exposure")
+            backward_transform = Mapping((2, 1)) | wcs2.backward_transform[2] & wcs2.backward_transform[1]
+        else:
+            backward_transform = wcs2.backward_transform
     elif issubclass(wcs2, Model):
         backward_transform = wcs2.inverse
     else:
@@ -158,3 +166,7 @@ def build_mask(dqarr, bitvalue):
     if bitvalue is None:
         return (np.ones(dqarr.shape, dtype=np.uint8))
     return np.logical_not(np.bitwise_and(dqarr, ~bitvalue)).astype(np.uint8)
+
+def is_sky_like(frame):
+    # Differentiate between sky-like and cartesian frames
+    return u.Unit("deg") in frame.unit or u.Unit("arcsec") in frame.unit