JP-3547: Flux conservation for spectral resampling

CHANGES.rst

@@ -26,6 +26,23 @@ outlier_detection
 - Fixed failures due to a missing ``wcs.array_shape`` attribute when the
   ``outlier_detection`` step was run standalone using e.g. ``strun`` [#8645]
 
+resample_spec
+-------------
+
+- Modified the output NIRSpec spectral WCS to sample the input data linearly in sky
+  coordinates, rather than slit coordinates, in order to conserve spectral
+  flux in default reductions. [#8596]
+
+- Updated handling for the ``pixel_scale_ratio`` parameter to apply only to the
+  spatial dimension, to match the sense of the parameter application to the
+  documented intent, and to conserve spectral flux when applied. [#8596]
+
+- Implemented handling for the ``pixel_scale`` parameter, which was previously
+  ignored for spectral resampling. [#8596]
+
+- Fixed a bug resulting in incorrect output slit coordinates for NIRSpec moving
+  targets in the ``calwebb_spec3`` pipeline. [#8596]
+
 stpipe
 ------
 

docs/jwst/resample/arguments.rst

@@ -15,35 +15,73 @@ image.
     image.  Available kernels are `square`, `gaussian`, `point`,
     `turbo`, `lanczos2`, and `lanczos3`.
 
+    For spectral data, only the `square` and `point` kernels should be used.
+    The other kernels do not conserve spectral flux.
+
 ``--pixel_scale_ratio`` (float, default=1.0)
-    Ratio of input to output pixel scale.  A value of 0.5 means the output
+    Ratio of input to output pixel scale.
+
+    For imaging data, a value of 0.5 means the output
     image would have 4 pixels sampling each input pixel.
+
+    For spectral data, values greater than 1 indicate that the input
+    pixels have a larger spatial scale, so more output pixels will
+    sample the same input pixel.  For example, a value of 2.0
+    means the output image would have 2 pixels sampling each input
+    spatial pixel. If the input data has units of flux density (MJy/pixel),
+    the output flux per pixel will be half the input flux per pixel.
+    If the input data has units of surface brightness (MJy/sr), the output
+    flux per pixel is not scaled.
+
+    Note that this parameter is only applied in the cross-dispersion
+    direction for spectral data: sampling wavelengths are not affected.
+
     Ignored when ``pixel_scale`` or ``output_wcs`` are provided.
 
+    .. note::
+        If this parameter is modified for spectral data, the extraction
+        aperture for the :ref:`extract_1d <extract_1d_step>` step must
+        also be modified, since it is specified in pixels.
+
 ``--pixel_scale`` (float, default=None)
     Absolute pixel scale in ``arcsec``. When provided, overrides
     ``pixel_scale_ratio``. Ignored when ``output_wcs`` is provided.
 
+    For spectral data, if the input data has units of flux density
+    (MJy/pixel), the output flux per pixel will be scaled by the ratio
+    of the selected output pixel scale to an average input pixel scale.
+    If the input data has units of surface brightness (MJy/sr),
+    the output flux per pixel is not scaled.
+
+    Note that this parameter is only applied in the cross-dispersion
+    direction for spectral data: sampling wavelengths are not affected.
+
+    .. note::
+        If this parameter is modified for spectral data, the extraction
+        aperture for the :ref:`extract_1d <extract_1d_step>` step must
+        also be modified, since it is specified in pixels.
+
 ``--rotation`` (float, default=None)
     Position angle of output image’s Y-axis relative to North.
     A value of 0.0 would orient the final output image to be North up.
     The default of `None` specifies that the images will not be rotated,
     but will instead be resampled in the default orientation for the camera
     with the x and y axes of the resampled image corresponding
     approximately to the detector axes. Ignored when ``pixel_scale``
-    or ``output_wcs`` are provided.
+    or ``output_wcs`` are provided.  Also ignored for all spectral data.
 
 ``--crpix`` (tuple of float, default=None)
     Position of the reference pixel in the image array in the ``x, y`` order.
     If ``crpix`` is not specified, it will be set to the center of the bounding
     box of the returned WCS object. When supplied from command line, it should
     be a comma-separated list of floats. Ignored when ``output_wcs``
-    is provided.
+    is provided. Also ignored for all spectral data.
 
 ``--crval`` (tuple of float, default=None)
     Right ascension and declination of the reference pixel. Automatically
     computed if not provided. When supplied from command line, it should be a
     comma-separated list of floats. Ignored when ``output_wcs`` is provided.
+    Also ignored for all spectral data.
 
 ``--output_shape`` (tuple of int, default=None)
     Shape of the image (data array) using "standard" ``nx`` first and ``ny``
@@ -117,6 +155,8 @@ image.
     For example, if set to ``0.5``, only resampled images that use less than
     half the available memory can be created.
 
+    This parameter is ignored for spectral data.
+
 ``--in_memory`` (boolean, default=True)
   Specifies whether or not to load and create all images that are used during
   processing into memory. If ``False``, input files are loaded from disk when

docs/jwst/resample/resample.rst

@@ -1,6 +1,7 @@
 .. resample_:
 
-Python Interface to Drizzle: ResampleData()
-===========================================
+Python Interface to Drizzle: ResampleData() and ResampleSpecData()
+==================================================================
 
 .. automodapi:: jwst.resample.resample
+.. automodapi:: jwst.resample.resample_spec

docs/jwst/resample/resample_step.rst

@@ -1,6 +1,7 @@
 .. resample_step_:
 
-Python Step Interface: ResampleStep()
-=====================================
+Python Step Interface: ResampleStep() and ResampleSpecStep()
+============================================================
 
 .. automodapi:: jwst.resample.resample_step
+.. automodapi:: jwst.resample.resample_spec_step

jwst/outlier_detection/outlier_detection.py

@@ -572,10 +572,16 @@ def gwcs_blot(median_model, blot_img, interp='poly5', sinscl=1.0):
     # Set array shape, needed to compute image pixel area
     blot_img.meta.wcs.array_shape = blot_img.shape
     if 'SPECTRAL' not in blot_img.meta.wcs.output_frame.axes_type:
+        # Account for intensity scaling, needed if there is a difference
+        # between nominal pixel area and average pixel area,
+        # computed from the WCS
         input_pixflux_area = blot_img.meta.photometry.pixelarea_steradians
         input_pixel_area = compute_image_pixel_area(blot_img.meta.wcs)
         pix_ratio = np.sqrt(input_pixflux_area / input_pixel_area)
     else:
+        # Note: spectral scaling is only needed if the pixel ratio
+        # is not set to 1.0, which is not supported for
+        # outlier detection.
         pix_ratio = 1.0
     log.info('Blotting {} <-- {}'.format(blot_img.data.shape, median_model.data.shape))
 

jwst/regtest/test_nirspec_fs_spec3.py

@@ -46,24 +46,13 @@ def test_nirspec_fs_spec3(run_pipeline, rtdata_module, fitsdiff_default_kwargs,
     diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
     assert diff.identical, diff.report()
 
-    if output == "s2d":
-        # Compare the calculated wavelengths
+    # Check output wavelength array against its own wcs
+    if suffix == "s2d":
         tolerance = 1e-03
-        dmt = datamodels.open(rtdata.truth)
         dmr = datamodels.open(rtdata.output)
-        if isinstance(dmt, datamodels.MultiSlitModel):
-            names = [s.name for s in dmt.slits]
-            for name in names:
-                st_idx = [(s.wcs, s.wavelength) for s in dmt.slits if s.name==name]
-                w = dmt.slits[st_idx].meta.wcs
-                x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
-                _, _, wave = w(x, y)
-                sr_idx = [(s.wcs, s.wavelength) for s in dmr.slits if s.name==name]
-                wlr = dmr.slits[sr_idx].wavelength
-                assert np.all(np.isclose(wave, wlr, atol=tolerance))
-        else:
-            w = dmt.meta.wcs
-            x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
-            _, _, wave = w(x, y)
-            wlr = dmr.wavelength
-            assert np.all(np.isclose(wave, wlr, atol=tolerance))
+
+        w = dmr.meta.wcs
+        x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
+        _, _, wave = w(x, y)
+        wlr = dmr.wavelength
+        assert np.all(np.isclose(wave, wlr, atol=tolerance))

jwst/regtest/test_nirspec_fs_spec3_moving_target.py

@@ -0,0 +1,57 @@
+from astropy.io.fits.diff import FITSDiff
+import pytest
+import numpy as np
+from gwcs import wcstools
+
+from jwst.stpipe import Step
+from stdatamodels.jwst import datamodels
+
+
+@pytest.fixture(scope="module")
+def run_pipeline(rtdata_module):
+    """
+    Run the calwebb_spec3 pipeline on a NIRSpec FS moving target.
+    """
+    rtdata = rtdata_module
+
+    # Get the ASN file and input exposures
+    rtdata.get_asn('nirspec/fs/jw01245-o002_20240701t053319_spec3_00001_asn.json')
+
+    # Run the calwebb_spec3 pipeline; save results from intermediate steps
+    args = ["calwebb_spec3", rtdata.input,
+            "--steps.outlier_detection.save_results=true",
+            "--steps.resample_spec.save_results=true",
+            "--steps.extract_1d.save_results=true"]
+    Step.from_cmdline(args)
+
+
+@pytest.mark.bigdata
+@pytest.mark.parametrize("suffix", ["cal", "crf", "s2d", "x1d"])
+def test_nirspec_fs_spec3_moving_target(
+        run_pipeline, rtdata_module, fitsdiff_default_kwargs, suffix):
+    """Test spec3 pipeline on a NIRSpec FS moving target."""
+    rtdata = rtdata_module
+
+    output = f"jw01245-o002_s000000001_nirspec_clear-prism-s200a1-subs200a1_{suffix}.fits"
+    rtdata.output = output
+    rtdata.get_truth(f"truth/test_nirspec_fs_spec3_moving_target/{output}")
+
+    # Adjust tolerance for machine precision with float32 drizzle code
+    if suffix == "s2d":
+        fitsdiff_default_kwargs["rtol"] = 1e-2
+        fitsdiff_default_kwargs["atol"] = 2e-4
+
+    # Compare the results
+    diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
+    assert diff.identical, diff.report()
+
+    # Check output wavelength array against its own wcs
+    if suffix == "s2d":
+        tolerance = 1e-03
+        dmr = datamodels.open(rtdata.output)
+
+        w = dmr.meta.wcs
+        x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
+        _, _, wave = w(x, y)
+        wlr = dmr.wavelength
+        assert np.all(np.isclose(wave, wlr, atol=tolerance))

jwst/regtest/test_nirspec_mos_spec3.py

@@ -42,24 +42,13 @@ def test_nirspec_mos_spec3(run_pipeline, suffix, source_id, fitsdiff_default_kwa
     diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs)
     assert diff.identical, diff.report()
 
+    # Check output wavelength array against its own wcs
     if suffix == "s2d":
-        # Compare the calculated wavelengths
         tolerance = 1e-03
-        dmt = datamodels.open(rtdata.truth)
         dmr = datamodels.open(rtdata.output)
-        if isinstance(dmt, datamodels.MultiSlitModel):
-            names = [s.name for s in dmt.slits]
-            for name in names:
-                st_idx = [(s.wcs, s.wavelength) for s in dmt.slits if s.name==name]
-                w = dmt.slits[st_idx].meta.wcs
-                x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
-                _, _, wave = w(x, y)
-                sr_idx = [(s.wcs, s.wavelength) for s in dmr.slits if s.name==name]
-                wlr = dmr.slits[sr_idx].wavelength
-                assert np.all(np.isclose(wave, wlr, atol=tolerance))
-        else:
-            w = dmt.meta.wcs
-            x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
-            _, _, wave = w(x, y)
-            wlr = dmr.wavelength
-            assert np.all(np.isclose(wave, wlr, atol=tolerance))
+
+        w = dmr.meta.wcs
+        x, y = wcstools.grid_from_bounding_box(w.bounding_box, step=(1, 1), center=True)
+        _, _, wave = w(x, y)
+        wlr = dmr.wavelength
+        assert np.all(np.isclose(wave, wlr, atol=tolerance))