(JP-3337) Renaming of undersampling_correction to charge_migration. (s…

…pacetelescope#7825)

Co-authored-by: Howard Bushouse <[email protected]>

CHANGES.rst

@@ -28,6 +28,13 @@ calwebb_spec2
   that is returned by the pipeline to ensure a file is created with the
   expected ``_cal`` suffix. [#7772]
 
+charge_migration
+----------------
+
+- Step was renamed from undersampling_migration. Changed default signal threshold,
+  added efficient routine to flag neighborhood pixels, added new unit test,
+  improved earlier unit tests, updated docs. [#7825]
+
 cube_build
 ----------
 
@@ -110,7 +117,7 @@ residual_fringe
 
 - Use scipy.interpolate.BSpline instead of astropy.modeling.Spline1D in
   residual_fringe fitting utils [#7764]
-
+  
 undersampling_correction
 ------------------------
 

docs/jwst/charge_migration/arguments.rst

@@ -0,0 +1,7 @@
+Arguments
+=========
+
+The ``charge migration`` step has one optional argument that can be set by the user:
+
+* ``--signal_threshold``: A floating-point value in units of ADU for the science value above which
+  a group's DQ will be flagged as CHARGELOSS and DO_NOT_USE.

docs/jwst/charge_migration/description.rst

@@ -0,0 +1,51 @@
+Description
+===========
+
+:Class: `jwst.charge_migration.ChargeMigrationStep`
+:Alias: charge_migration
+
+
+Overview
+--------
+This step corrects for an artifact seen in undersampled NIRISS images that may depress flux 
+in resampled images. The artifact is seen in dithered images where the star is centered in 
+a pixel. When the peak pixels of such stars approach the saturation level, they suffer from 
+significant :ref:`charge migration <charge_migration>`:
+the spilling of charge from a saturated pixel into its neighboring pixels. This charge migration 
+causes group-to-group differences to decrease significantly once the signal level is greater than 
+~25,000 ADU. As a result, the last several groups of these ramps get flagged by the :ref:`jump <jump_step>`
+step.  The smaller number of groups used for these pixels in the :ref:`ramp_fitting <ramp_fitting_step>`
+step results in them having  larger read noise variances, which in turn leads to lower weights used
+during resampling. This ultimately leads to a lower than normal flux for the star in resampled images.
+
+Once a group in a ramp has been flagged as affected by charge migration, all subsequent 
+groups in the ramp are also flagged. By flagging these groups, they are not used in the 
+computation of slopes in the :ref:`ramp_fitting <ramp_fitting_step>` step. However, as described 
+in the algorithm section below, they _are_ used in the calculation of the variance of the slope 
+due to readnoise.
+
+Input details
+-------------
+The input must be in the form of a `~jwst.datamodels.RampModel`.
+
+
+Algorithm
+--------- 
+The first group, and all subsequent groups, exceeding the value of the 
+``signal_threshold`` parameter is flagged as CHARGELOSS. ``signal_threshold`` is in units 
+of ADUs. These groups will also be flagged as DO_NOT_USE, and will not 
+be included in the slope calculation during the ``ramp_fitting`` step. Despite being flagged 
+as DO_NOT_USE, these CHARGELOSS groups are still included in the calculation of the
+variance due to readnoise. 
+This results in a readnoise variance for undersampled pixels that is similar to that of 
+pixels unaffected by charge migration. For the Poisson noise variance calculation in 
+:ref:`ramp_fitting <ramp_fitting_step>`, the CHARGELOSS/DO_NOT_USE groups are not included.
+
+For integrations having only 1 or 2 groups, no flagging will be performed.
+
+
+Output product
+--------------
+The output is a new copy of the input `~jwst.datamodels.RampModel`, with the updated DQ flags
+added to the GROUPDQ array.
+

docs/jwst/charge_migration/index.rst

@@ -0,0 +1,14 @@
+.. _charge_migration_step:
+
+================
+Charge Migration
+================
+
+.. toctree::
+   :maxdepth: 2
+
+   description.rst
+   arguments.rst
+   reference_files.rst
+
+.. automodapi:: jwst.charge_migration

docs/jwst/charge_migration/reference_files.rst

@@ -0,0 +1,3 @@
+Reference Files
+===============
+This step does not use any reference files.

docs/jwst/package_index.rst

@@ -14,6 +14,7 @@ Package Index
    background_step/index.rst
    background_subtraction/index.rst
    barshadow/index.rst
+   charge_migration/index.rst
    combine_1d/index.rst
    cube_build/index.rst
    dark_current/index.rst
@@ -63,7 +64,6 @@ Package Index
    superbias/index.rst
    tso_photometry/index.rst
    tweakreg/index.rst
-   undersampling_correction/index.rst
    wavecorr/index.rst
    wfs_combine/index.rst
    wfss_contam/index.rst

docs/jwst/pipeline/calwebb_detector1.rst

@@ -72,8 +72,8 @@ on either the first group or frame zero pixel values.
 +----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
 |                                                                      |         |         | :ref:`refpix <refpix_step>`             | |check| | |check| |
 +----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
-| :ref:`undersampling_correction <undersampling_correction_step>` [3]_ | |check| |         |                                         |         |         |
-+----------------------------------------------------------------------+---------+---------+----------------------+------------------+---------+---------+
+| :ref:`charge_migration <charge_migration_step>` [3]_                 | |check| |         |                                         |         |         |
++----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
 | :ref:`jump <jump_step>`                                              | |check| | |check| | :ref:`jump <jump_step>`                 | |check| | |check| |
 +----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
 | :ref:`ramp_fitting <ramp_fitting_step>`                              | |check| | |check| | :ref:`ramp_fitting <ramp_fitting_step>` | |check| | |check| |
@@ -86,7 +86,7 @@ on either the first group or frame zero pixel values.
    retrieved from CRDS will skip the :ref:`ipc <ipc_step>` step for all instruments.
 .. [2] The :ref:`persistence <persistence_step>` step is currently hardwired to be skipped in
    the `Detector1Pipeline` module for all NIRSpec exposures.
-.. [3] By default, the :ref:`undersampling_correction <undersampling_correction_step>` step is skipped in
+.. [3] By default, the :ref:`charge_migration <charge_migration_step>` step is skipped in
    the `Detector1Pipeline` module for all instruments.
    
 

docs/jwst/ramp_fitting/description.rst

@@ -322,9 +322,9 @@ extension.
 
 Weighted Readnoise Variance
 +++++++++++++++++++++++++++
-If the :ref:`undersampling correction <undersampling_correction_step>`
+If the :ref:`charge migration <charge_migration_step>`
 step has been performed prior to ramp fitting, any group whose value exceeds the
-`signal_threshold` parameter will have been flagged with the UNDERSAMP and DO_NOT_USE
+`signal_threshold` parameter will have been flagged with the CHARGELOSS and DO_NOT_USE
 data quality flags. Due to the DO_NOT_USE flags, such groups will be excluded
 from the slope calculations.
 

docs/jwst/references_general/references_general.rst

@@ -602,7 +602,7 @@ Bit  Value         Name              Description
 4    16            OUTLIER           Flagged by outlier detection
 5    32            PERSISTENCE       High persistence
 6    64            AD_FLOOR          Below A/D floor
-7    128           UNDERSAMP         Undersampling correction
+7    128           CHARGELOSS        Charge Migration
 8    256           UNRELIABLE_ERROR  Uncertainty exceeds quoted error
 9    512           NON_SCIENCE       Pixel not on science portion of detector
 10   1024          DEAD              Dead pixel

jwst/charge_migration/__init__.py

@@ -0,0 +1,3 @@
+from .charge_migration_step import ChargeMigrationStep
+
+__all__ = ['ChargeMigrationStep']

jwst/charge_migration/charge_migration.py

@@ -0,0 +1,143 @@
+#  Module for charge migration
+#
+import logging
+import numpy as np
+
+from stdatamodels.jwst.datamodels import dqflags
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+GOOD = dqflags.group["GOOD"]
+DNU = dqflags.group["DO_NOT_USE"]
+CHLO = dqflags.group["CHARGELOSS"]
+
+CHLO_DNU = CHLO + DNU
+
+
+def charge_migration(input_model, signal_threshold):
+    """
+    Correct for chargemigration
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodels.RampModel`
+        The input science data to be corrected
+
+    signal_threshold : float
+        Science value above which a group will be flagged as CHARGELOSS
+
+    Returns
+    -------
+    output_model : `~jwst.datamodels.RampModel`
+        Data model with charge_migration applied; add CHARGELOSS and
+
+        DO_NOT_USE flags to groups exceeding signal_threshold
+    """
+    data = input_model.data
+    gdq = input_model.groupdq
+
+    # Create the output model as a copy of the input
+    output_model = input_model.copy()
+
+    log.info('Using signal_threshold: %.2f', signal_threshold)
+
+    gdq_new = flag_pixels(data, gdq, signal_threshold)
+
+    # Save the flags in the output GROUPDQ array
+    output_model.groupdq = gdq_new
+
+    return output_model
+
+
+def flag_pixels(data, gdq, signal_threshold):
+    """
+    Flag each group in each ramp that exceeds signal_threshold as CHARGELOSS and DO_NOT_USE,
+    skipping groups already flagged as DO_NOT_USE.
+
+    Parameters
+    ----------
+    data : float, 4D array
+        science array
+
+    gdq : int, 4D array
+        group dq array
+
+    signal_threshold : float
+        Science value above which a group will be flagged as CHARGELOSS and DO_NOT_USE
+
+    Returns
+    -------
+    new_gdq : int, 4D array
+        updated group dq array
+    """
+    n_ints, n_grps, n_rows, n_cols = gdq.shape
+
+    new_gdq = gdq.copy()   # Updated gdq
+
+    # Flag all exceedances with CHARGELOSS and NO_NOT_USE
+    chargeloss_pix = (data > signal_threshold) & (gdq != DNU)
+    new_gdq[chargeloss_pix] = np.bitwise_or(new_gdq[chargeloss_pix], CHLO | DNU)
+
+    # Reset groups previously flagged as DNU
+    gdq_orig = gdq.copy()  # For resetting to previously flagged DNU
+    wh_gdq_DNU = np.bitwise_and(gdq_orig, DNU)
+
+    # Get indices for exceedances
+    arg_where = np.argwhere(new_gdq == CHLO_DNU)
+
+    a_int = arg_where[:, 0]  # array of integrations
+    a_grp = arg_where[:, 1]  # array of groups
+    a_row = arg_where[:, 2]  # array of rows
+    a_col = arg_where[:, 3]  # array of columns
+
+    # Process the 4 nearest neighbors of each exceedance
+    # Pixel to the east
+    xx_max_p1 = a_col[a_col < (n_cols-1)] + 1
+    i_int = a_int[a_col < (n_cols-1)]
+    i_grp = a_grp[a_col < (n_cols-1)]
+    i_row = a_row[a_col < (n_cols-1)]
+
+    if len(xx_max_p1) > 0:
+        new_gdq[i_int, i_grp, i_row, xx_max_p1] = \
+            np.bitwise_or(new_gdq[i_int, i_grp, i_row, xx_max_p1], CHLO | DNU)
+
+    new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1]  # reset for earlier DNUs
+
+    # Pixel to the west
+    xx_m1 = a_col[a_col > 0] - 1
+    i_int = a_int[a_col > 0]
+    i_grp = a_grp[a_col > 0]
+    i_row = a_row[a_col > 0]
+
+    if len(xx_m1) > 0:
+        new_gdq[i_int, i_grp, i_row, xx_m1] = \
+            np.bitwise_or(new_gdq[i_int, i_grp, i_row, xx_m1], CHLO | DNU)
+
+    new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1]  # reset for earlier DNUs
+
+    # Pixel to the north
+    yy_m1 = a_row[a_row > 0] - 1
+    i_int = a_int[a_row > 0]
+    i_grp = a_grp[a_row > 0]
+    i_col = a_col[a_row > 0]
+
+    if len(yy_m1) > 0:
+        new_gdq[i_int, i_grp, yy_m1, i_col] = \
+            np.bitwise_or(new_gdq[i_int, i_grp, yy_m1, i_col], CHLO | DNU)
+
+    new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1]  # reset for earlier DNUs
+
+    # Pixel to the south
+    yy_max_p1 = a_row[a_row < (n_rows-1)] + 1
+    i_int = a_int[a_row < (n_rows-1)]
+    i_grp = a_grp[a_row < (n_rows-1)]
+    i_col = a_col[a_row < (n_rows-1)]
+
+    if len(yy_max_p1) > 0:
+        new_gdq[i_int, i_grp, yy_max_p1, i_col] = \
+            np.bitwise_or(new_gdq[i_int, i_grp, yy_max_p1, i_col], CHLO | DNU)
+
+    new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1]  # reset for earlier DNUs
+
+    return new_gdq

jwst/charge_migration/charge_migration_step.py

@@ -0,0 +1,46 @@
+#! /usr/bin/env python
+import logging
+
+from ..stpipe import Step
+
+from stdatamodels.jwst import datamodels
+
+from . import charge_migration
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+__all__ = ["ChargeMigrationStep"]
+
+
+class ChargeMigrationStep(Step):
+    """
+    This Step sets the CHARGELOSS flag for groups exhibiting significant
+    charge migration.
+    """
+    class_alias = "charge_migration"
+
+    spec = """
+        signal_threshold = float(default=30000)
+        skip = boolean(default=True)
+    """
+
+    def process(self, input):
+
+        # Open the input data model
+        with datamodels.RampModel(input) as input_model:
+            if (input_model.data.shape[1] < 3):  # skip step if only 1 or 2 groups/integration
+                log.info('Too few groups per integration; skipping charge_migration')
+
+                result = input_model
+                result.meta.cal_step.charge_migration = 'SKIPPED'
+
+                return result
+
+            # Retrieve the parameter value(s)
+            signal_threshold = self.signal_threshold
+
+            result = charge_migration.charge_migration(input_model, signal_threshold)
+            result.meta.cal_step.charge_migration = 'COMPLETE'
+
+        return result