JP-4196: Update NIS_WFSS exptype to expect per-wavelength units in photom ref file

changes/10199.photom.rst

@@ -0,0 +1 @@
+Update NIS_WFSS exptype to expect per-wavelength units in photom ref file

changes/10199.wfss_contam.rst

@@ -0,0 +1 @@
+Update NIS_WFSS exptype to expect per-wavelength units in photom ref file

jwst/photom/photom.py

@@ -884,7 +884,7 @@ def photom_io(self, tabdata, order=None, time_correction=None):
                     )
                     slit.photom_point = conversion  # store the result
 
-                elif self.exptype in ["NRC_WFSS", "NRC_TSGRISM"]:
+                elif self.exptype in ["NRC_WFSS", "NRC_TSGRISM", "NIS_WFSS"]:
                     log.info("Including spectral dispersion in 2-d flux calibration")
                     conversion, no_cal = self.create_2d_conversion(
                         slit,
@@ -909,7 +909,7 @@ def photom_io(self, tabdata, order=None, time_correction=None):
                 )
             else:
                 # NRC_TSGRISM data produces a SpecModel, which is handled here
-                if self.exptype in ["NRC_WFSS", "NRC_TSGRISM"]:
+                if self.exptype in ["NRC_WFSS", "NRC_TSGRISM", "NIS_WFSS"]:
                     log.info("Including spectral dispersion in 2-d flux calibration")
                     conversion, no_cal = self.create_2d_conversion(
                         self.input,
@@ -1101,14 +1101,15 @@ def create_2d_conversion(
             dispaxis = get_dispersion_direction(self.exptype, self.grating, self.filter, self.pupil)
             if dispaxis is not None:
                 dispersion_array = self.get_dispersion_array(wl_array, dispaxis)
-                # Convert dispersion from micron/pixel to angstrom/pixel
-                dispersion_array *= 1.0e4
+                if self.exptype != "NIS_WFSS":
+                    # Convert dispersion from micron/pixel to angstrom/pixel
+                    # except for NIRISS WFSS, which has conv2d in micron/pixel
+                    dispersion_array *= 1.0e4
                 conv_2d /= np.abs(dispersion_array)
             else:
                 log.warning(
                     "Unable to get dispersion direction, so cannot calculate dispersion array"
                 )
-        # Combine the scalar and 2D conversion factors
         conversion = conversion * conv_2d
         no_cal = np.isnan(conv_2d)
         conversion[no_cal] = 0.0

jwst/photom/tests/test_photom.py

@@ -1499,6 +1499,8 @@ def test_niriss_wfss():
         input_data = slit.data  # this is from save_input
         output = ds.input.slits[k].data  # ds.input is the output
         sp_order = slit.meta.wcsinfo.spectral_order
+
+        # retrieve relevant photom table data
         rownum = find_row_in_ftab(
             save_input, ftab, ["filter", "pupil"], slitname=None, order=sp_order
         )
@@ -1510,8 +1512,15 @@ def test_niriss_wfss():
         ix = shape[1] // 2
         iy = shape[0] // 2
         wl = slit.wavelength[iy, ix]
+
+        # compute dispersion
+        dispersion_array = np.gradient(slit.wavelength[:, ix])
+        dispersion = dispersion_array[iy]
+
+        # compute expected value
         rel_resp = np.interp(wl, wavelength, relresponse, left=np.nan, right=np.nan)
-        compare = photmjsr * rel_resp
+        compare = photmjsr * rel_resp / np.abs(dispersion)
+
         # Compare the values at the center pixel.
         ratio = output[iy, ix] / input_data[iy, ix]
         result.append(np.allclose(ratio, compare, rtol=1.0e-7))

jwst/wfss_contam/disperse.py

@@ -199,7 +199,6 @@ def disperse(
     grism_wcs,
     naxis,
     oversample_factor=2,
-    phot_per_lam=True,
 ):
     """
     Compute the dispersed image pixel values from the direct image.
@@ -211,7 +210,8 @@ def disperse(
     ys : ndarray
         Flat array of Y coordinates of pixels in the direct image
     fluxes : ndarray
-        Fluxes of the pixels in the direct image corresponding to xs, ys
+        Fluxes of the pixels in the direct image corresponding to xs, ys.
+        These should have units of MJy/sr.
     source_ids_per_pixel : int array
         Source IDs of the input pixels in the segmentation map
     order : int
@@ -221,9 +221,10 @@ def disperse(
     wmax : float
         Maximum wavelength for dispersed spectra
     sens_waves : float array
-        Wavelength array from photom reference file
+        Wavelength array from photom reference file. Expected unit is micron.
     sens_resp : float array
-        Response (flux calibration) array from photom reference file
+        Response (flux calibration) array from photom reference file.
+        Expected units are (micron) * (MJy / sr) / (ADU/s).
     direct_image_wcs : WCS object
         WCS object for the direct image and segmentation map
     grism_wcs : WCS object
@@ -232,12 +233,6 @@ def disperse(
         Dimensions of the grism image (naxis[0], naxis[1])
     oversample_factor : int, optional
         Factor by which to oversample the wavelength grid
-    phot_per_lam : bool
-        If True, it is assumed that the photometric response is calibrated per wavelength,
-        as is done for NIRCam, and therefore we need to multiply by dlam to un-do the calibration
-        in this step.
-        If False, it is assumed that the response is calibrated per pixel, and there is no need
-        to multiply by dlam. This is what's done for NIRISS.
 
     Returns
     -------
@@ -323,16 +318,11 @@ def disperse(
     # The input direct image data is already photometrically calibrated,
     # so we need to basically apply a reverse flux calibration here.
     # Divide out the response values to convert from Mjy/sr to DN/s.
-    # Note that the photom reference files are constructed differently for NIRCam and NIRISS
-    # in the way they handle the dispersion: NIRCam uses per/wavelength, NIRISS uses per-pixel.
+    # Note that the photom reference files are constructed with per-wavelength units,
+    # so oversampling is accounted for by the spacing of dlam.
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", category=RuntimeWarning, message="divide by zero")
-        if phot_per_lam:
-            # NIRCam case. Oversampling is accounted for by the spacing of dlam.
-            counts = fluxes * areas * dlam * 1e4 / sens
-        else:
-            # NIRISS case. Oversampling must be handled directly to avoid double-counting.
-            counts = fluxes * areas / (sens * oversample_factor)
+        counts = fluxes * areas * dlam / sens
     counts[no_cal] = 0.0  # set to zero where no flux cal info available
     del fluxes, areas, sens, dlam, no_cal
 

jwst/wfss_contam/observations.py

@@ -124,7 +124,6 @@ def __init__(
         max_cpu=1,
         max_pixels_per_chunk=5e4,
         oversample_factor=2,
-        phot_per_lam=True,
     ):
         """
         Initialize all data and metadata for a given observation.
@@ -147,10 +146,6 @@ def __init__(
             Maximum number of pixels per chunk when dispersing sources
         oversample_factor : int, optional
             Factor by which to oversample the wavelength grid
-        phot_per_lam : bool, optional
-            Whether to compute photometry per wavelength bin (True) or per pixel (False).
-            This depends on how the photom reference file has been delivered.
-            True should be used for NIRCam, and False should be used for NIRISS.
         """
         if boundaries is None:
             boundaries = []
@@ -164,7 +159,6 @@ def __init__(
         self.max_cpu = max_cpu
         self.max_pixels_per_chunk = max_pixels_per_chunk
         self.oversample_factor = oversample_factor
-        self.phot_per_lam = phot_per_lam
 
         # ensure the direct image has background subtracted
         self.dimage = background_subtract(direct_image)
@@ -274,7 +268,6 @@ def chunk_sources(
                     self.grism_wcs,
                     self.naxis,
                     self.oversample_factor,
-                    self.phot_per_lam,
                 ]
             )
 

jwst/wfss_contam/sens1d.py

@@ -1,6 +1,7 @@
 import logging
 
 import numpy as np
+import stdatamodels.jwst.datamodels as dm
 
 from jwst.photom.photom import find_row
 
@@ -14,6 +15,10 @@ def get_photom_data(phot_model, filter_name, pupil, order):
     Retrieve wavelength and response data from photom ref file.
 
     Wavelengths from the reference file are expected to be in units of microns.
+    Units of the relative response values depend on the instrument:
+    NIRCam: (Angstrom) * (MJy/sr) / (ADU/s)
+    NIRISS: (micron) * (MJy / sr) / (ADU/s)
+    Output units are converted to (micron) * (MJy / sr) / (ADU/s).
 
     Parameters
     ----------
@@ -32,7 +37,7 @@ def get_photom_data(phot_model, filter_name, pupil, order):
         Wavelengths from the ref file.
     relresps : float array
         Wavelength-dependent response (flux calibration) values from the ref file,
-        same shape as `ref_waves`.
+        same shape as `ref_waves`, units of (micron) * (MJy / sr) / (ADU/s).
     """
     # Get the appropriate row of data from the reference table
     phot_table = phot_model.phot_table
@@ -41,7 +46,15 @@ def get_photom_data(phot_model, filter_name, pupil, order):
     tabdata = phot_table[row]
 
     # Scalar conversion factor
-    scalar_conversion = tabdata["photmjsr"]  # unit is MJy / sr
+    if isinstance(phot_model, dm.NrcWfssPhotomModel):
+        # NIRCam: convert from Angstrom to micron
+        scalar_conversion = tabdata["photmjsr"] / 1e4
+    elif isinstance(phot_model, dm.NisWfssPhotomModel):
+        # NIRISS: no conversion needed
+        scalar_conversion = tabdata["photmjsr"]
+    else:
+        log.error("Unsupported photom model type: %s", type(phot_model))
+        raise TypeError("Unsupported photom model type")
 
     # Get the length of the relative response arrays in this row
     nelem = tabdata["nelem"]

jwst/wfss_contam/tests/conftest.py

@@ -109,7 +109,7 @@ def grism_wcs():
 
 
 @pytest.fixture(scope="module")
-def photom_ref_model():
+def phot_table():
     """
     Make a mock photom reference model for NIRISS WFSS.
 
@@ -158,4 +158,56 @@ def photom_ref_model():
     phot_table["wavelength"] = wavelength
     phot_table["relresponse"] = relresponse
     phot_table["reluncertainty"] = reluncertainty
+    return phot_table
+
+
+@pytest.fixture(scope="module")
+def photom_ref_model_niriss(phot_table):
+    """
+    Make a mock photom reference model for NIRISS WFSS.
+
+    Parameters
+    ----------
+    phot_table : np.recarray
+        Photometry table.
+
+    Returns
+    -------
+    `~stdatamodels.jwst.datamodels.NisWfssPhotomModel`
+        Photom ref file model.
+    """
     return dm.NisWfssPhotomModel(phot_table=phot_table)
+
+
+@pytest.fixture(scope="module")
+def photom_ref_model_nircam(phot_table):
+    """
+    Make a mock photom reference model for NIRCam WFSS.
+
+    Parameters
+    ----------
+    phot_table : np.recarray
+        Photometry table.
+
+    Returns
+    -------
+    `~stdatamodels.jwst.datamodels.NrcWfssPhotomModel`
+        Photom ref file model.
+    """
+    return dm.NrcWfssPhotomModel(phot_table=phot_table)
+
+
+@pytest.fixture
+def photom_ref_model(request):
+    """
+    Make a mock photom reference model that can be parameterized for any WFSS mode.
+
+    If request.param is not set, defaults to NIRISS.
+
+    Returns
+    -------
+    `~stdatamodels.jwst.datamodels.JwstDataModel`
+        Photom ref file model of the appropriate type based on the parameterization.
+    """
+    fixture_name = getattr(request, "param", "photom_ref_model_niriss")
+    return request.getfixturevalue(fixture_name)

jwst/wfss_contam/tests/test_disperse.py

@@ -1,12 +1,10 @@
 import numpy as np
-import pytest
 from numpy.testing import assert_allclose
 
 from jwst.wfss_contam.disperse import disperse
 
 
-@pytest.mark.parametrize("phot_per_lam", [True, False])
-def test_disperse_oversample_same_result(grism_wcs, direct_image_with_gradient, phot_per_lam):
+def test_disperse_oversample_same_result(grism_wcs, direct_image_with_gradient):
     """Coverage for bug where wavelength oversampling led to double-counted fluxes."""
     x0 = np.array([200.5])
     y0 = np.array([200.5])
@@ -35,7 +33,6 @@ def test_disperse_oversample_same_result(grism_wcs, direct_image_with_gradient,
             grism_wcs,
             naxis,
             oversample_factor=os,
-            phot_per_lam=phot_per_lam,
         )
         output_images.append(src[source_id[0]]["image"])
 

jwst/wfss_contam/tests/test_observations.py

@@ -24,7 +24,6 @@ def observation(direct_image_with_gradient, segmentation_map, grism_wcs):
         segmentation_map.data,
         grism_wcs,
         direct_image_with_gradient.meta.wcs,
-        phot_per_lam=False,
     )
 
 
@@ -115,4 +114,4 @@ def test_disperse_order(observation, segmentation_map):
     assert slit.data.shape == (slit.ysize, slit.xsize)
 
     # check for regression by hard-coding one value of slit.data
-    assert np.isclose(slit.data[5, 60], 20.996877670288086)
+    assert np.isclose(slit.data[5, 60], 0.09994397)

jwst/wfss_contam/tests/test_sens_waves.py → jwst/wfss_contam/tests/test_sens1d.py

@@ -1,8 +1,13 @@
 import numpy as np
+import pytest
+import stdatamodels.jwst.datamodels as dm
 
 from jwst.wfss_contam.sens1d import create_1d_sens, get_photom_data
 
 
+@pytest.mark.parametrize(
+    "photom_ref_model", ["photom_ref_model_niriss", "photom_ref_model_nircam"], indirect=True
+)
 def test_get_photom_data(photom_ref_model):
     filter_name = "GR150C"
     pupil = "F200W"
@@ -28,8 +33,13 @@ def test_get_photom_data(photom_ref_model):
 
     # get the flux scaling out of the table: F200W is the 3rd row
     scalar_conversion = refmodel.phot_table["photmjsr"][2]
-    assert np.isclose(ref_relresp[0], 10 * scalar_conversion)
-    assert np.isclose(ref_relresp[-1], 99 * scalar_conversion)
+    if isinstance(photom_ref_model, dm.NrcWfssPhotomModel):
+        # NIRCam gets scaled by 1e-4 to account for per-inverse-angstrom units
+        assert np.isclose(ref_relresp[0], 10 * 1e-4 * scalar_conversion)
+        assert np.isclose(ref_relresp[-1], 99 * 1e-4 * scalar_conversion)
+    else:
+        assert np.isclose(ref_relresp[0], 10 * scalar_conversion)
+        assert np.isclose(ref_relresp[-1], 99 * scalar_conversion)
 
 
 def test_create_1d_sens(photom_ref_model):

jwst/wfss_contam/wfss_contam.py

@@ -435,10 +435,8 @@ def contam_corr(
     # filter needs to come from the PUPIL keyword value.
     if input_model.meta.instrument.name == "NIRISS":
         filter_name = pupil_kwd
-        phot_per_lam = False
     else:
         filter_name = filter_kwd
-        phot_per_lam = True
 
     # Read the source catalog to perform magnitude-based source selection later
     # mag limit will be scaled according to order 1 sensitivity
@@ -459,7 +457,6 @@ def contam_corr(
         max_cpu=ncpus,
         max_pixels_per_chunk=max_pixels_per_chunk,
         oversample_factor=oversample_factor,
-        phot_per_lam=phot_per_lam,
     )
 
     no_sources = True