JP-3593: Saturation step should increase flagging for group 3 saturation

CHANGES.rst

@@ -271,6 +271,11 @@ pixel_replace
 - Added estimated errors and variances for replaced pixels, following the
   interpolation scheme used for the data. [#8504]
 
+saturation
+----------
+
+- Adds a check for saturation bias in group 2 for nframes > 1. [#8499]
+
 ramp_fitting
 ------------
 

jwst/saturation/saturation.py

@@ -49,6 +49,8 @@ def flag_saturation(input_model, ref_model, n_pix_grow_sat):
     """
 
     data = input_model.data
+    ngroups = input_model.meta.exposure.ngroups
+    nframes = input_model.meta.exposure.nframes
 
     # Create the output model as a copy of the input
     output_model = input_model.copy()
@@ -68,9 +70,10 @@ def flag_saturation(input_model, ref_model, n_pix_grow_sat):
         sat_dq = ref_sub_model.dq.copy()
         ref_sub_model.close()
 
+    read_pattern = [[x + 1 + groupstart * nframes for x in range(nframes)] for groupstart in range(ngroups)]
     gdq_new, pdq_new, zframe = flag_saturated_pixels(
         data, gdq, pdq, sat_thresh, sat_dq, ATOD_LIMIT, dqflags.pixel,
-        n_pix_grow_sat=n_pix_grow_sat, zframe=zframe)
+        n_pix_grow_sat=n_pix_grow_sat, read_pattern=read_pattern, zframe=zframe)
 
     # Save the flags in the output GROUPDQ array
     output_model.groupdq = gdq_new
@@ -117,6 +120,8 @@ def irs2_flag_saturation(input_model, ref_model, n_pix_grow_sat):
     nints = data.shape[0]
     ngroups = data.shape[1]
     detector = input_model.meta.instrument.detector
+    nframes = input_model.meta.exposure.nframes
+    read_pattern = [[x + 1 + groupstart * nframes for x in range(nframes)] for groupstart in range(ngroups)]
 
     # create a mask of the appropriate size
     irs2_mask = x_irs2.make_mask(input_model)
@@ -160,6 +165,31 @@ def irs2_flag_saturation(input_model, ref_model, n_pix_grow_sat):
                                         irs2_mask, detector)
             # check for saturation
             flag_temp = np.where(sci_temp >= sat_thresh, SATURATED, 0)
+            if group == 2:
+                # Identify groups which we wouldn't expect to saturate by the third group,
+                # on the basis of the first group
+                scigp1 = x_irs2.from_irs2(data[ints, 0, :, :], irs2_mask, detector)
+                mask = scigp1 / np.mean(read_pattern[0]) * read_pattern[2][-1] < sat_thresh
+
+                # Identify groups with suspiciously large values in the second group
+                scigp2 = x_irs2.from_irs2(data[ints, 1, :, :], irs2_mask, detector)
+                mask &= scigp2 > sat_thresh / len(read_pattern[1])
+
+                # Identify groups that are saturated in the third group
+                gp3mask = np.where(flag_temp & SATURATED, True, False)
+                mask &= gp3mask
+
+                # Flag the 2nd group for the pixels passing that gauntlet
+                dq_temp = x_irs2.from_irs2(groupdq[ints, 1, :, :], irs2_mask, detector)
+                dq_temp[mask] |= SATURATED
+                # flag any pixels that border saturated pixels
+                if n_pix_grow_sat > 0:
+                    dq_temp = adjacency_sat(dq_temp, SATURATED, n_pix_grow_sat)
+                # set the flags in dq array for group 2, i.e. index 1
+                x_irs2.to_irs2(flagarray, dq_temp, irs2_mask, detector)
+                np.bitwise_or(groupdq[ints, 1, ...], flagarray,
+                              groupdq[ints, 1, ...])
+
             # check for A/D floor
             flaglow_temp = np.where(sci_temp <= 0, AD_FLOOR | DONOTUSE, 0)
 

jwst/saturation/tests/test_saturation.py

@@ -28,8 +28,8 @@ def test_basic_saturation_flagging(setup_nrc_cube):
     # Add ramp values up to the saturation limit
     data.data[0, 0, 5, 5] = 0
     data.data[0, 1, 5, 5] = 20000
-    data.data[0, 2, 5, 5] = 40000
-    data.data[0, 3, 5, 5] = 60000   # Signal reaches saturation limit
+    data.data[0, 2, 5, 5] = 60000   # Signal reaches saturation limit
+    data.data[0, 3, 5, 5] = 60000
     data.data[0, 4, 5, 5] = 62000
 
     # Set saturation value in the saturation model
@@ -386,6 +386,7 @@ def _cube(ngroups, nrows, ncols):
         data_model.meta.subarray.xsize = ncols
         data_model.meta.subarray.ysize = nrows
         data_model.meta.exposure.ngroups = ngroups
+        data_model.meta.exposure.nframes = 3
         data_model.meta.instrument.name = 'NIRCAM'
         data_model.meta.instrument.detector = 'NRCA1'
         data_model.meta.observation.date = '2017-10-01'
@@ -432,6 +433,8 @@ def _cube(xstart, ystart, ngroups, nrows, ncols):
         data_model.meta.subarray.xsize = ncols
         data_model.meta.subarray.ystart = ystart
         data_model.meta.subarray.ysize = nrows
+        data_model.meta.exposure.ngroups = ngroups
+        data_model.meta.exposure.nframes = 3
 
         # create a saturation model for the saturation step
         saturation_model = SaturationModel((1032, 1024))
@@ -477,6 +480,8 @@ def _cube():
         data_model.meta.exposure.nrs_normal = 16
         data_model.meta.exposure.nrs_reference = 4
         data_model.meta.exposure.readpatt = 'NRSIRS2RAPID'
+        data_model.meta.exposure.nframes = 3
+        data_model.meta.exposure.ngroups = 5
 
         # create a saturation model for the saturation step
         saturation_model = SaturationModel((2048, 2048))