Adjust saturation threshold depending on the read pattern

CHANGES.rst

@@ -43,6 +43,11 @@ refpix
 
 - Add initial reference pixel correction step implementation. [#704]
 
+saturation
+----------
+
+- Add read_pattern argument to flag_saturated_pixels. [#836]
+
 general
 -------
 

docs/roman/saturation/description.rst

@@ -8,7 +8,8 @@ below the A/D floor if they have a value of zero DN.
 
 This step examines the data group-by-group, comparing the pixel values in the data array with defined
 saturation thresholds for each pixel. When it finds a pixel value in a given
-group that is above the saturation threshold (high saturation), it sets the
+group that is above the saturation threshold (high saturation) times a
+dilution factor, it sets the
 "SATURATED" flag in the corresponding location of the "groupdq" array in the
 science exposure.  When it finds a pixel in a given group that has a zero or
 negative value (below  the A/D floor), it sets the "AD_FLOOR" and "DO_NOT_USE"
@@ -24,3 +25,13 @@ of 65535 and hence will only be flagged as saturated if the pixel reaches that
 hard limit in at least one group. The "NO_SAT_CHECK" flag is propagated to the
 PIXELDQ array in the output science data to indicate which pixels fall into
 this category.
+
+The "dilution factor" is intended to account for the fact that Roman
+downlinks resultants to the ground, which are usually averages over
+several reads.  The saturation threshold corresponds to the number of
+counts at which a detector pixel saturates.  Because a resultant is
+the average of a number of reads, later reads in a resultant can
+saturate, but if earlier reads are unsaturated, the value of the
+resultant can fall under the saturation threshold.  The dilution
+factor varies resultant-by-resultant and is given by
+:math:`<t>/max(t)` for all times :math:`t` entering a resultant.

pyproject.toml

@@ -28,7 +28,8 @@ dependencies = [
     'roman_datamodels ==0.17.0',
 #    'roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git@main',
     'scipy >=1.11',
-    'stcal >=1.4.0',
+#    'stcal >=1.4.0',
+    'stcal @ git+https://github.com/schlafly/stcal.git@saturation-read-pattern',
     'stpipe >=0.5.0',
     'tweakwcs >=0.8.0',
     'spherical-geometry >= 1.2.22',

requirements-dev.txt

@@ -10,5 +10,5 @@ git+https://github.com/spacetelescope/[email protected]
 git+https://github.com/spacetelescope/[email protected]
 scipy>=0.0.dev0
 git+https://github.com/spacetelescope/stpipe
-git+https://github.com/spacetelescope/stcal
+git+https://github.com/schlafly/stcal@saturation-read-pattern
 git+https://github.com/spacetelescope/tweakwcs

romancal/photom/tests/test_photom.py

@@ -60,7 +60,7 @@ def create_photom_wfi_image(min_r=3.1, delta=0.1):
 
     # Create sample area keyword values
     area_ster = 2.31307642258977e-14 * u.steradian
-    pixelareasr = np.ones(nrows, dtype=np.float32) * area_ster
+    pixelareasr = np.ones(nrows, dtype=np.float64) * area_ster
 
     # Bundle values into a list
     values = list(zip(photmjsr, uncertainty, pixelareasr))

romancal/saturation/saturation.py

@@ -49,7 +49,15 @@ def flag_saturation(input_model, ref_model):
     # The third variable is the processed ZEROFRAME, which is not
     # used in romancal, so is always None.
     gdq_new, pdq_new, _ = flag_saturated_pixels(
-        data, gdq, pdq, sat_thresh, sat_dq, ATOD_LIMIT, dqflags.pixel, n_pix_grow_sat=0
+        data,
+        gdq,
+        pdq,
+        sat_thresh,
+        sat_dq,
+        ATOD_LIMIT,
+        dqflags.pixel,
+        n_pix_grow_sat=0,
+        read_pattern=input_model.meta.exposure.read_pattern,
     )
 
     # Save the flags in the output GROUPDQ array

romancal/saturation/tests/test_saturation.py

@@ -41,6 +41,49 @@ def test_basic_saturation_flagging(setup_wfi_datamodels):
     assert np.all(output.groupdq[satindex:, 5, 5] == dqflags.group["SATURATED"])
 
 
+def test_read_pattern_saturation_flagging(setup_wfi_datamodels):
+    """Check that the saturation threshold varies depending on how the reads
+    are allocated into resultants."""
+
+    # Create inputs, and data and saturation models
+    ngroups = 5
+    nrows = 20
+    ncols = 20
+    satvalue = 60000
+    ramp, satmap = setup_wfi_datamodels(ngroups, nrows, ncols)
+
+    # Add ramp values up to the saturation limit
+    ramp.data[0, 5, 5] = 0 * ramp.data.unit
+    ramp.data[1, 5, 5] = 20000 * ramp.data.unit
+    ramp.data[2, 5, 5] = 40000 * ramp.data.unit
+    ramp.data[3, 5, 5] = 60000 * ramp.data.unit  # Signal reaches saturation limit
+    ramp.data[4, 5, 5] = 62000 * ramp.data.unit
+
+    # set read_pattern to have many reads in the third resultant, so that
+    # its mean exposure time is much smaller than its last read time
+    # (in this case, the ratio is 13 / 20).
+    # This means that the effective saturation for the third resultant
+    # is 60000 * 13 / 20 = 39000 and the third resultant should be marked
+    # saturated.
+    ramp.meta.exposure.read_pattern = [
+        [1],
+        [2],
+        [3, 4, 5, 6, 7, 8, 9, 10],
+        [11],
+        [12],
+        [13],
+    ]
+
+    # Set saturation value in the saturation model
+    satmap.data[5, 5] = satvalue * satmap.data.unit
+
+    # Run the pipeline
+    output = flag_saturation(ramp, satmap)
+
+    # Make sure that groups after the third get flagged
+    assert np.all(output.groupdq[2:, 5, 5] == dqflags.group["SATURATED"])
+
+
 def test_ad_floor_flagging(setup_wfi_datamodels):
     """Check that the ad_floor flag is set when a pixel value is zero or
     negative."""