JP-3593: 2nd group saturation part 2

CHANGES.rst

@@ -236,6 +236,12 @@ resample_spec
 - Ensure that NaNs and DO_NOT_USE flags match up in all input data before
   resampling. [#8557]
 
+saturation
+----------
+
+- Add option for using the readout pattern information to improve saturation flagging
+  in grouped data. [#8731]
+
 scripts
 -------
 

docs/jwst/saturation/arguments.rst

@@ -6,3 +6,8 @@ The ``saturation`` step has one optional argument:
   The distance to use when growing saturation flag values to neighboring pixels,
   in order to account for charge migration (spilling). The total region size is
   2*n_pix_grow_sat+1 pixels, centered on the primary pixel.
+
+``--use_readpatt`` (boolean, default=True)
+  Use information about the readout pattern for grouped data to perform additional
+  checks to find data that saturates mid-group.  This can be important for
+  finding extremely bright cosmic rays that saturate in group 2.

docs/jwst/saturation/description.rst

@@ -34,12 +34,13 @@ The "NO_SAT_CHECK" flag is propagated to the
 PIXELDQ array in the output science data to indicate which pixels fall into
 this category.
 
-For grouped data, this step also looks for and flags data that saturated during
+If the "use_readpatt" keyword is set, this step will use information about the
+read pattern to find pixels that saturated in the middle of grouped data.  This
+can be particularly important for flagging data that saturated during
 the second group but did not trigger the normal saturation threshold due to the
 grouped data averaging.  This requires that the third group be saturated, and
 the first group sufficiently low that the third group would not have been expected
 to saturate (i.e., flagging due to cosmic rays but not sources).
-At present, this is only implemented for NIRSpec IRS2 data.
 
 .. _charge_migration:
 

jwst/saturation/saturation.py

@@ -22,7 +22,7 @@
 ATOD_LIMIT = 65535.  # Hard DN limit of 16-bit A-to-D converter
 
 
-def flag_saturation(input_model, ref_model, n_pix_grow_sat):
+def flag_saturation(input_model, ref_model, n_pix_grow_sat, use_readpatt):
     """
     Short Summary
     -------------
@@ -41,6 +41,9 @@ def flag_saturation(input_model, ref_model, n_pix_grow_sat):
         as saturated (i.e '1' will flag the surrouding 8 pixels) to account for
         charge spilling.
 
+    use_readpatt : bool
+        Use grouped read pattern information to assist with flagging
+
     Returns
     -------
     output_model : `~jwst.datamodels.RampModel`
@@ -49,6 +52,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 +73,16 @@ def flag_saturation(input_model, ref_model, n_pix_grow_sat):
         sat_dq = ref_sub_model.dq.copy()
         ref_sub_model.close()
 
+    # Enable use of read_pattern specific treatment if selected
+    if use_readpatt:
+        read_pattern = [[x + 1 + groupstart * nframes for x in range(nframes)] for groupstart in range(ngroups)]
+        log.info(f"Using read_pattern with nframes {nframes}")
+    else:
+        read_pattern=None
+
     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
@@ -84,7 +96,7 @@ def flag_saturation(input_model, ref_model, n_pix_grow_sat):
     return output_model
 
 
-def irs2_flag_saturation(input_model, ref_model, n_pix_grow_sat):
+def irs2_flag_saturation(input_model, ref_model, n_pix_grow_sat, use_readpatt):
     """
     Short Summary
     -------------
@@ -106,6 +118,9 @@ def irs2_flag_saturation(input_model, ref_model, n_pix_grow_sat):
         as saturated (i.e '1' will flag the surrouding 8 pixels) to account for
         charge spilling.
 
+    use_readpatt : bool
+        Use grouped read pattern information to assist with flagging
+
     Returns
     -------
     output_model : `~jwst.datamodels.RampModel`
@@ -118,7 +133,12 @@ def irs2_flag_saturation(input_model, ref_model, n_pix_grow_sat):
     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)]
+
+    if use_readpatt:
+        read_pattern = [[x + 1 + groupstart * nframes for x in range(nframes)] for groupstart in range(ngroups)]
+        log.info(f"Using read_pattern with nframes {nframes}")
+    else:
+        read_pattern=None
 
     # create a mask of the appropriate size
     irs2_mask = x_irs2.make_mask(input_model)
@@ -162,7 +182,8 @@ 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:
+            # Additional checks for group 2 saturation in grouped data
+            if ((group == 2) & (read_pattern is not None)):
                 # 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)

jwst/saturation/saturation_step.py

@@ -18,6 +18,7 @@ class SaturationStep(Step):
 
     spec = """
         n_pix_grow_sat = integer(default=1) # number of layers adjacent pixels to flag
+        use_readpatt = boolean(default=True) # Use grouped read pattern information to assist with flagging
     """
 
     reference_file_types = ['saturation']
@@ -44,9 +45,9 @@ def process(self, input):
 
             # Do the saturation check
             if pipe_utils.is_irs2(input_model):
-                sat = saturation.irs2_flag_saturation(input_model, ref_model, self.n_pix_grow_sat)
+                sat = saturation.irs2_flag_saturation(input_model, ref_model, self.n_pix_grow_sat, self.use_readpatt)
             else:
-                sat = saturation.flag_saturation(input_model, ref_model, self.n_pix_grow_sat)
+                sat = saturation.flag_saturation(input_model, ref_model, self.n_pix_grow_sat, self.use_readpatt)
 
             # Close the reference file and update the step status
             ref_model.close()

jwst/saturation/tests/test_saturation.py

@@ -36,7 +36,7 @@ def test_basic_saturation_flagging(setup_nrc_cube):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Make sure that groups with signal > saturation limit get flagged
     satindex = np.argmax(output.data[0, :, 4:6, 4:6] == satvalue)
@@ -51,12 +51,53 @@ def test_nirspec_irs2_saturation_flagging(setup_nrs_irs2_cube):
     data.data[0, 4, pixx, pixy] = 67000
 
     # Run saturation detection
-    output = irs2_flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = irs2_flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Make sure that groups with signal > saturation limit get flagged
     assert np.all(output.groupdq[0, 3:, pixx - 1: pixx + 1, pixy - 1: pixy + 1] == dqflags.group['SATURATED'])
 
 
+def test_nirspec_irs2_readpatt(setup_nrs_irs2_cube):
+    """
+    Tests that the readpatt framework (saturation in grouped data) is working for IRS2 processing.
+    """
+    data, satmap = setup_nrs_irs2_cube()
+
+    pixx, pixy = 1000, 1000
+    data.data[0, 0, pixx, pixy] = 500  # Low signal
+    data.data[0, 1, pixx, pixy] = 50000  # Suspiciously high signal
+    data.data[0, 2, pixx, pixy] = 65000  # Signal exceeds saturation limit of 60000
+    data.data[0, 3, pixx, pixy] = 66000
+    data.data[0, 4, pixx, pixy] = 67000
+
+    # Run saturation detection
+    output = irs2_flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=True)
+
+    # Make sure that group 2 gets flagged
+    assert np.all(output.groupdq[0, 1, pixx, pixy] == dqflags.group['SATURATED'])
+
+
+def test_nirspec_nrs_readpatt(setup_nrs_nrs_cube):
+    """
+    Tests that the readpatt framework (saturation in grouped data) is working for regular grouped
+    data (code follows all instruments except NIRSpec IRS2)
+    """
+    data, satmap = setup_nrs_nrs_cube()
+
+    pixx, pixy = 1000, 1000
+    data.data[0, 0, pixx, pixy] = 500  # Low signal
+    data.data[0, 1, pixx, pixy] = 50000  # Suspiciously high signal
+    data.data[0, 2, pixx, pixy] = 65000  # Signal exceeds saturation limit of 60000
+    data.data[0, 3, pixx, pixy] = 66000
+    data.data[0, 4, pixx, pixy] = 67000
+
+    # Run saturation detection
+    output = irs2_flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=True)
+
+    # Make sure that group 2 gets flagged
+    assert np.all(output.groupdq[0, 1, pixx, pixy] == dqflags.group['SATURATED'])
+
+
 def test_irs2_zero_frame(setup_nrs_irs2_cube):
     """
     Tests IRS2 ZEROFRAME processing.
@@ -76,7 +117,7 @@ def test_irs2_zero_frame(setup_nrs_irs2_cube):
     ramp.zeroframe[nint, row, col] = 65000.
     ramp.zeroframe[nint, row + 1, col + 1] = -100.
 
-    output = irs2_flag_saturation(ramp, sat, n_pix_grow_sat=0)
+    output = irs2_flag_saturation(ramp, sat, n_pix_grow_sat=0, use_readpatt=False)
 
     check_zframe = np.ones(dims, dtype=ramp.data.dtype) * 1000.
     check_zframe[nint, row, col] = 0.
@@ -111,7 +152,7 @@ def test_ad_floor_flagging(setup_nrc_cube):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Check if the right frames are flagged as saturated
     assert np.all(output.groupdq[0, satindxs, 5, 5]
@@ -146,7 +187,7 @@ def test_ad_floor_and_saturation_flagging(setup_nrc_cube):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Check if the right frames are flagged as ad_floor
     assert np.all(output.groupdq[0, floorindxs, 5, 5]
@@ -180,7 +221,7 @@ def test_signal_fluctuation_flagging(setup_nrc_cube):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Make sure that all groups after first saturated group are flagged
     satindex = np.argmax(output.data[0, :, 5, 5] == satvalue)
@@ -210,7 +251,7 @@ def test_all_groups_saturated(setup_nrc_cube):
     satmap.data[5, 5] = satvalue
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Make sure all groups are flagged
     assert np.all(np.bitwise_and(output.groupdq[0, :, 4:6, 4:6],
@@ -237,7 +278,7 @@ def test_subarray_extraction(setup_miri_cube):
     satmap.dq[550, 100] = dqflags.pixel['NONLINEAR']
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Check for DQ flag in PIXELDQ of subarray image
     assert output.pixeldq[76, 100] == dqflags.pixel['DO_NOT_USE']
@@ -267,7 +308,7 @@ def test_dq_propagation(setup_nrc_cube):
     satmap.dq[5, 5] = dqval2
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Make sure DQ values from data and reference file are added in the output
     assert output.pixeldq[5, 5] == dqval1 + dqval2
@@ -300,7 +341,7 @@ def test_no_sat_check(setup_nrc_cube):
     data.pixeldq[5, 5] = dqflags.pixel['RC']
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Make sure output GROUPDQ does not get flagged as saturated
     # Make sure PIXELDQ is set to NO_SAT_CHECK and original flag
@@ -332,7 +373,7 @@ def test_nans_in_mask(setup_nrc_cube):
     satmap.data[5, 5] = np.nan
 
     # Run the pipeline
-    output = flag_saturation(data, satmap, n_pix_grow_sat=1)
+    output = flag_saturation(data, satmap, n_pix_grow_sat=1, use_readpatt=False)
 
     # Check that output GROUPDQ is not flagged as saturated
     assert np.all(output.groupdq[0, :, 4:6, 4:6] != dqflags.group['SATURATED'])
@@ -501,3 +542,49 @@ def _cube():
 
         return data_model, saturation_model
     return _cube
+
+@pytest.fixture(scope='function')
+def setup_nrs_nrs_cube():
+
+    def _cube():
+
+        # create a JWST datamodel for NIRSPEC NRS-mode data
+        data_model = RampModel((1, 5, 3200, 2048))
+        data_model.data = np.ones(((1, 5, 3200, 2048)))
+        data_model.groupdq = np.zeros(((1, 5, 3200, 2048)))
+        data_model.pixeldq = np.zeros(((3200, 2048)))
+        data_model.meta.instrument.name = 'NIRSPEC'
+        data_model.meta.instrument.detector = 'NRS1'
+        data_model.meta.instrument.filter = 'F100LP'
+        data_model.meta.observation.date = '2019-07-19'
+        data_model.meta.observation.time = '23:23:30.912'
+        data_model.meta.exposure.type = 'NRS_LAMP'
+        data_model.meta.subarray.name = 'FULL'
+        data_model.meta.subarray.xstart = 1
+        data_model.meta.subarray.xsize = 2048
+        data_model.meta.subarray.ystart = 1
+        data_model.meta.subarray.ysize = 2048
+        data_model.meta.exposure.nrs_normal = 16
+        data_model.meta.exposure.nrs_reference = 4
+        data_model.meta.exposure.readpatt = 'NRS'
+        data_model.meta.exposure.nframes = 4
+        data_model.meta.exposure.ngroups = 5
+
+        # create a saturation model for the saturation step
+        saturation_model = SaturationModel((2048, 2048))
+        saturation_model.data = np.ones((2048, 2048)) * 60000  # saturation limit for every pixel is 60000
+        saturation_model.meta.description = 'Fake data.'
+        saturation_model.meta.telescope = 'JWST'
+        saturation_model.meta.reftype = 'SaturationModel'
+        saturation_model.meta.useafter = '2015-10-01T00:00:00'
+        saturation_model.meta.instrument.name = 'NIRSPEC'
+        saturation_model.meta.instrument.detector = 'NRS1'
+        saturation_model.meta.author = 'David'
+        saturation_model.meta.pedigree = 'Dummy'
+        saturation_model.meta.subarray.xstart = 1
+        saturation_model.meta.subarray.xsize = 2048
+        saturation_model.meta.subarray.ystart = 1
+        saturation_model.meta.subarray.ysize = 2048
+
+        return data_model, saturation_model
+    return _cube