Merge branch 'main' into lastframe_codestyle

.pre-commit-config.yaml

@@ -40,14 +40,12 @@ repos:
             jwst/assign_wcs/.* |
             jwst/associations/.* |
             jwst/background/.* |
-            jwst/badpix_selfcal/.* |
             jwst/charge_migration/.* |
             jwst/combine_1d/.* |
             jwst/coron/.* |
             jwst/cube_build/.* |
             jwst/cube_skymatch/.* |
             jwst/dark_current/.* |
-            jwst/datamodels/.* |
             jwst/dq_init/.* |
             jwst/emicorr/.* |
             jwst/engdblog/.* |
@@ -84,12 +82,10 @@ repos:
             jwst/skymatch/.* |
             jwst/spectral_leak/.* |
             jwst/srctype/.* |
-            jwst/stpipe/.* |
             jwst/straylight/.* |
             jwst/superbias/.* |
             jwst/tests/.* |
             jwst/tso_photometry/.* |
-            jwst/tweakreg/.* |
             jwst/wavecorr/.* |
             jwst/wfs_combine/.* |
             jwst/wfss_contam/.* |

.ruff.toml

@@ -27,16 +27,16 @@ exclude = [
     "jwst/assign_wcs/**.py",
     "jwst/associations/**.py",
     "jwst/background/**.py",
-    "jwst/badpix_selfcal/**.py",
-    #"jwst/barshadow/**.py",
+    # "jwst/badpix_selfcal/**.py",
+    # "jwst/barshadow/**.py",
     "jwst/charge_migration/**.py",
     # "jwst/clean_flicker_noise/**.py",
     "jwst/combine_1d/**.py",
     "jwst/coron/**.py",
     "jwst/cube_build/**.py",
     "jwst/cube_skymatch/**.py",
     "jwst/dark_current/**.py",
-    "jwst/datamodels/**.py",
+    # "jwst/datamodels/**.py",
     "jwst/dq_init/**.py",
     "jwst/emicorr/**.py",
     "jwst/engdblog/**.py",
@@ -79,12 +79,12 @@ exclude = [
     # "jwst/source_catalog/**.py",
     "jwst/spectral_leak/**.py",
     "jwst/srctype/**.py",
-    "jwst/stpipe/**.py",
+    # "jwst/stpipe/**.py",
     "jwst/straylight/**.py",
     "jwst/superbias/**.py",
     "jwst/tests/**.py",
     "jwst/tso_photometry/**.py",
-    "jwst/tweakreg/**.py",
+    # "jwst/tweakreg/**.py",
     "jwst/wavecorr/**.py",
     "jwst/wfs_combine/**.py",
     "jwst/wfss_contam/**.py",
@@ -157,16 +157,16 @@ ignore-fully-untyped = true  # Turn of annotation checking for fully untyped cod
 "jwst/assign_wcs/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/associations/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/background/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
-"jwst/badpix_selfcal/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
-#"jwst/barshadow/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
+# "jwst/badpix_selfcal/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
+# "jwst/barshadow/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/charge_migration/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 # "jwst/clean_flicker_noise/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/combine_1d/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/coron/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/cube_build/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/cube_skymatch/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/dark_current/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
-"jwst/datamodels/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
+# "jwst/datamodels/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/dq_init/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/emicorr/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/engdblog/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
@@ -209,12 +209,12 @@ ignore-fully-untyped = true  # Turn of annotation checking for fully untyped cod
 # "jwst/source_catalog/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/spectral_leak/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/srctype/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
-"jwst/stpipe/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
+# "jwst/stpipe/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/straylight/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/superbias/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/tests/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/tso_photometry/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
-"jwst/tweakreg/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
+# "jwst/tweakreg/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/wavecorr/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/wfs_combine/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]
 "jwst/wfss_contam/**.py" = ["D", "N", "A", "ARG", "B", "C4", "ICN", "INP", "ISC", "LOG", "NPY", "PGH", "PTH", "S", "SLF", "SLOT", "T20", "TRY", "UP", "YTT", "E501"]

changes/9126.straylight.rst

@@ -0,0 +1 @@
+Add an additional optional step to the MIRI MRS straylight correction code to measure and remove residual cosmic ray showers using information from the inter-slice pixels.

docs/jwst/references_general/references_general.rst

@@ -169,6 +169,8 @@ documentation on each reference file.
 |                                                       | :ref:`ABVEGAOFFSET <abvegaoffset_reffile>`       |
 +-------------------------------------------------------+--------------------------------------------------+
 | :ref:`straylight <straylight_step>`                   | :ref:`MRSXARTCORR <mrsxartcorr_reffile>`         |
++                                                       +--------------------------------------------------+
+|                                                       | :ref:`REGIONS <regions_reffile>`                 |
 +-------------------------------------------------------+--------------------------------------------------+
 | :ref:`spectral_leak <spectral_leak_step>`             | :ref:`MRSPTCORR <mrsptcorr_reffile>`             |
 +-------------------------------------------------------+--------------------------------------------------+
@@ -273,6 +275,8 @@ documentation on each reference file.
 | :ref:`REGIONS <regions_reffile>`                 | :ref:`assign_wcs <assign_wcs_step>`                   |
 +                                                  +-------------------------------------------------------+
 |                                                  | :ref:`residual_fringe <residual_fringe_step>`         |
++                                                  +-------------------------------------------------------+
+|                                                  | :ref:`straylight <straylight_step>`                   |
 +--------------------------------------------------+-------------------------------------------------------+
 | :ref:`RESET <reset_reffile>`                     | :ref:`reset <reset_step>`                             |
 +--------------------------------------------------+-------------------------------------------------------+

docs/jwst/references_general/regions_reffile.rst

@@ -28,6 +28,8 @@ The IFU takes a region reference file that defines the pixels on the detector th
 IFU slice they correspond to.  The reference file is a three-dimensional array
 measuring 1032 x 1024 x 9.  Each of the nine planes within the array represents
 a two-dimensional detector image for which each pixel has a value indicating the slice to which it belongs.
+Negative values represent slice pixels that receive light from the sky but for data quality reasons
+are not mapped to the sky using the world coordinate solution.
 The nine planes each correspond to different throughput thresholds ranging from 10% - 90%; these differ slightly
 in the effective along-slice size.
 

docs/jwst/straylight/arguments.rst

@@ -1,3 +1,32 @@
 Step Arguments
 ==============
-The ``straylight`` step has no step-specific arguments.
+The ``straylight`` step has the following optional arguments.
+
+``--clean_showers`` (boolean, default=False)
+  A flag to indicate whether to remove straylight due to residual cosmic
+  ray showers using pixels between the IFU slices.
+
+``--shower_plane`` (int, default=3)
+  Identifies the throughput plane to use from the MRS regions reference
+  files to identify between-slice pixels.  Lower values identify pixels
+  with less science signal, but result in fewer (or no) pixels between
+  slices in some bands.
+
+``--shower_x_stddev`` (float, default=18)
+   Standard deviation to use for gaussian kernel convolution in the X
+   (between-slice) direction.  Must be large enough to reach the middle
+   of the IFU slices.
+
+``--shower_y_stddev`` (float, default=5)
+   Standard deviation to use for gaussian kernel convolution in the Y
+   (along-slice) direction.
+
+``--shower_low_reject`` (float, default=0.1)
+   Percentile of low-valued pixels to reject when performing gaussian
+   kernel convolution.  Important to ensure unmasked bad pixels do
+   not bias the kernel convolution.
+
+``--shower_high_reject`` (float, default=99.9)
+   Percentile of high-valued pixels to reject when performing gaussian
+   kernel convolution.  Important to ensure unmasked bad pixels do
+   not bias the kernel convolution.

docs/jwst/straylight/main.rst

@@ -10,18 +10,19 @@ The ``straylight`` correction is only valid for MIRI MRS data.
 
 Overview
 --------
-This routine removes contamination of MIRI MRS spectral data
-by the MIRI cross artifact feature produced by internal reflections
-within the detector arrays.  As discussed in depth for the MIRI Imager
-by A. Gáspár et al. 2021 (PASP, 133, 4504), the cross artifact manifests
+This routine removes contamination of MIRI MRS spectral data due to a variety of stray light
+sources including the MIRI cross artifact, variable detector dark current, and residual cosmic ray showers.
+
+As discussed in depth for the MIRI Imager
+by A. Gáspár et al. 2021 (PASP, 133, 4504), the MIRI cross artifact feature is produced by internal reflections
+within the detector arrays.  For the MIRI MRS, the cross artifact manifests
 as a signal extending up to hundreds of pixels along the detector column and row directions from
 bright sources.  This signal has both smooth and structured components whose
 profiles vary as a function of wavelength.
 Although the peak intensity of the cross artifact is at
 most 1% of the source intensity in Channel 1 (decreasing toward longer wavelengths),
 the total integrated light in this feature can be of order 20% of the total light from a given source.
 
-
 In the MIRI MRS, such a signal extending along detector rows is more disruptive
 than for the MIRI imager.
 Since the individual IFU slices are interleaved on the detector
@@ -32,11 +33,18 @@ The purpose of this routine is thus to model the cross artifact feature in a giv
 and subtract it at the detector level prior to reformatting
 the data into three-dimensional cubes.
 
-At the same time, this step also ensures that the median count rate (in DN/s) in regions of the detector that
-see no direct light from the sky is zero for consistency with the applied flux calibration vectors.
+This step also ensures that the median count rate (in DN/s) in regions of the detector that
+see no direct light from the sky is zero for consistency with the applied flux calibration vectors.  This
+correction is necessary because the MRS darks are strongly time-variable, and significant signal can remain
+after subtraction of the reference dark in the ``dark_current`` step.
 
-Algorithm
----------
+Finally, this step optionally applies a correction for residual cosmic ray showers that can produce
+blotchy straylight features across the detector.  While the ``jump`` detection step masks the cores of the largest
+such artifacts, the extended halos of showers can cover hundreds of pixels and are best corrected as
+coherent extended stray light artifacts.
+
+Algorithm: Cross-Artifact Correction
+------------------------------------
 The basic idea of the cross artifact correction is to convolve a given science detector image with a
 kernel function that has been pre-calibrated based on observations
 of isolated sources and subtract the corresponding convolved image.
@@ -85,3 +93,25 @@ as the magnitude of the correction is typically too small to be visible from poi
 channels use only a Lorentzian kernel with the Gaussian amplitudes set to zero as such structured components are less
 obvious at these longer wavelengths.  In Channel 4 no correction appears to be necessary,
 and the amplitudes of all model components are set equal to zero.
+
+Algorithm: Residual Shower Correction
+-------------------------------------
+The residual cosmic ray shower correction relies upon the fact that cosmic rays shower artifacts tend to
+have extremely large, blobby halos that extend across one or more IFU slices and the inter-slice pixels between
+them.  These inter-slice pixels receive little to no direct light from the sky, and can be used to construct
+a model of residual showers which are coherent in detector space.
+
+Algorithmically, this routine first uses the :ref:`REGIONS <regions_reffile>` reference file to locate and
+mask all pixels on the detector that receive direct light from the sky.  Pixels that are flagged as
+DO_NOT_USE or REFERENCE_PIXEL are similarly masked, as are some highest and lowest percentile of pixel values
+to help ensure robustness against unmasked bad pixels.
+
+The unmasked remaining pixels are then used to construct a model of the detector rates in DN/s across
+both the masked and unmasked regions using a 2D Gaussian kernel.  This model of the residual shower artifacts
+is then subtracted from the data.
+
+Since the inter-slice pixels never reach truly zero throughput from the sky (particularly at short wavelengths)
+this correction is most reliable in cases for which the observations contain only faint sources and
+have had a pixel-based dedicated background subtraction applied.  While the correction can be applied to brighter
+sources, this may result in artifacts as the subtracted model may be dominated more by spilled source flux
+than by residual cosmic ray showers.

docs/jwst/straylight/reference_files.rst

@@ -6,4 +6,9 @@ for each MRS band.  In Channel 1 these vectors include power in a broad Lorentzi
 plus a pair of double-Gaussian profiles.  In Channels 2 and 3 these vectors include only
 power in the broad Lorentzian, while in Channel 4 there is no correction.
 
+Optionally, the ``straylight`` step can also use the :ref:`REGIONS <regions_reffile>`
+reference files during correction of residual cosmic ray showers.  These reference files
+describe the footprint of the dispersed IFU traces on the detector and are used to
+identify inter-slice pixels that do not see direct light from the sky.
+
 .. include:: ../references_general/mrsxartcorr_reffile.inc

jwst/badpix_selfcal/__init__.py

@@ -1,4 +1,6 @@
+"""Self-calibration of bad pixels in JWST data."""
+
 from .badpix_selfcal_step import BadpixSelfcalStep
 
 
-__all__ = ['BadpixSelfcalStep']
+__all__ = ["BadpixSelfcalStep"]

jwst/badpix_selfcal/badpix_selfcal.py

@@ -8,13 +8,15 @@
 from stdatamodels.jwst.datamodels.dqflags import pixel
 
 
-def badpix_selfcal(minimg: np.ndarray,
-                   flagfrac_lower: float = 0.001,
-                   flagfrac_upper: float = 0.001,
-                   kernel_size: int = 15,
-                   dispaxis=None) -> np.ndarray:
+def badpix_selfcal(
+    minimg: np.ndarray,
+    flagfrac_lower: float = 0.001,
+    flagfrac_upper: float = 0.001,
+    kernel_size: int = 15,
+    dispaxis=None,
+) -> np.ndarray:
     """
-    Flag residual artifacts as bad pixels in the DQ array of a JWST exposure
+    Flag residual artifacts as bad pixels in the DQ array of a JWST exposure.
 
     Parameters
     ----------
@@ -56,15 +58,19 @@ def badpix_selfcal(minimg: np.ndarray,
     minimg_hpf = minimg - smoothed
 
     # Flag outliers using percentile cutoff
-    flag_low, flag_high = np.nanpercentile(minimg_hpf, [flagfrac_lower * 100, (1 - flagfrac_upper) * 100])
+    flag_low, flag_high = np.nanpercentile(
+        minimg_hpf, [flagfrac_lower * 100, (1 - flagfrac_upper) * 100]
+    )
     bad = (minimg_hpf > flag_high) | (minimg_hpf < flag_low)
     flagged_indices = np.where(bad)
     return flagged_indices
 
 
 def apply_flags(input_model: IFUImageModel, flagged_indices: np.ndarray) -> IFUImageModel:
     """
-    Apply the flagged indices to the input model. Sets the flagged pixels to NaN
+    Apply the flagged indices to the input model.
+
+    Sets the flagged pixels to NaN
     and the DQ flag to DO_NOT_USE + OTHER_BAD_PIXEL
 
     Parameters
@@ -80,7 +86,6 @@ def apply_flags(input_model: IFUImageModel, flagged_indices: np.ndarray) -> IFUI
     output_model : IFUImageModel
         Flagged data model
     """
-
     input_model.dq[flagged_indices] |= pixel["DO_NOT_USE"] + pixel["OTHER_BAD_PIXEL"]
 
     input_model.data[flagged_indices] = np.nan