JP-3437: Address alternating column noise for NIRSpec IRS2 (#8143)

Co-authored-by: Howard Bushouse <bushouse@stsci.edu>

CHANGES.rst

@@ -141,6 +141,9 @@ refpix
 - Add option for NIRSpec IRS2 to preserve interleaved reference pixels in the
   output file, for calibration and diagnostic purposes. [#8255]
 
+- Add option to correct for mean reference pixel offsets by amplifier and detector
+  column in NIRSpec IRS2 mode. [#8143]
+
 resample
 --------
 

docs/jwst/refpix/description.rst

@@ -176,7 +176,7 @@ reference pixel correction for data read out using the IRS2 readout pattern.
 See the JdoxIRS2_ page for for an overview, and see Rauscher2017_ for
 details.
 
-The raw data include both the science data and interspersed reference
+The raw data include both the science data and interleaved reference
 pixel values.  The time to read out the entire detector includes not only
 the time to read each pixel of science ("normal") data and some of the
 reference pixels, but also time for the transition between reading normal
@@ -185,25 +185,22 @@ each row and between frames.  For example, it takes the same length of time
 to jump from reading normal pixels to reading reference pixels as it does
 to read one pixel value, about ten microseconds.
 
-Before subtracting the reference pixel and reference output values from
-the science data, some processing is done on the reference values, and the
-CRDS reference file factors are applied.  IRS2 readout is only used for
-full-frame data, never for subarrays.  The full detector is read out
-by four separate amplifiers simultaneously, and the reference output is
-read at the same time.  Each of these five readouts is the same size,
-640 by 2048 pixels (for IRS2).  If the CRDS reference file includes a
-DQ (data quality) BINTABLE extension, interleaved reference pixel values
-will be set to zero if they are flagged as bad in the DQ extension.
-
-At this point the algorithm looks for intermittently bad (or suspicious)
+IRS2 readout is only used for full-frame data, never for subarrays. The full
+detector is read out by four separate amplifiers simultaneously, and the
+reference output is read at the same time.  Each of these five readouts is
+the same size, 640 by 2048 pixels, each containing a repeating set of 8
+normal pixel readouts, 4 interleaved reference pixel readouts, and 8 more
+normal pixel readouts.
+
+The first step in processing IRS2 data is to look for intermittently bad
 reference pixels. This is done by calculating the means and standard
 deviations per reference pixel column, as well as the absolute value of the
 difference between readout pairs, across all groups within each integration.
 The robust mean and standard deviation of each of these arrays is then
 computed. Values greater than the robust mean plus the standard
 deviation, times a factor to avoid overcorrection, are flagged as bad
 pixels.  Readout pairs are always flagged together, and are flagged across
-all groups and integrations. Bad values are replaced by values from the
+all groups and integrations. Bad values will be replaced by values from the
 nearest reference group within the same amplifier, respecting parity
 (even/oddness).  The replacement value is the average of upper and lower
 values if both are good, or directly using the upper or lower values if only
@@ -212,6 +209,24 @@ good adjacent neighbor that does not match parity, that value is used.  If
 there are no good replacement values, the bad pixel is set to 0.0 to be
 interpolated over in the IRS2 correction to follow.
 
+After flagging bad reference pixels, the step performs an optional
+correction for overall mean reference pixel offsets by amplifier and
+column parity. The algorithm described above for the traditional NIR readout
+mode is applied to IRS2 data to perform this correction, with two small
+differences:
+
+    #. Side pixel correction is never applied for IRS2 data.
+
+    #. "Even" and "odd" refer to detector column addresses, rather than
+       data array locations, to ensure that interleaved reference pixel
+       columns are accounted for correctly.
+
+After the mean offsets are subtracted and bad pixels are replaced, some processing
+is done on the remaining reference values, and the CRDS reference file
+factors are applied. If the CRDS reference file includes a DQ (data quality)
+BINTABLE extension, interleaved reference pixel values will be set to zero if
+they are flagged as bad in the DQ extension.
+
 The next step in this processing is to
 copy the science data and the reference pixel data separately to temporary
 1-D arrays (both of length 712 * 2048); this is done separately for each

jwst/refpix/irs2_subtract_reference.py

@@ -607,37 +607,37 @@ def flag_bad_refpix(datamodel, n_sigma=3.0, flag_only=False, replace_only=False)
                         rp_stds.append(rp_s)
                         rp_diffs.append(rp_d)
 
-            ref_pix = np.array(ref_pix, dtype=int)
-            rp_diffs = np.array(rp_diffs)
-            rp_means = np.array(rp_means)
-            rp_stds = np.array(rp_stds)
-            pair_pixel = ref_pix + 1
-
-            # clipped stats for all tests
-            mean_of_diffs, _, std_of_diffs = sigma_clipped_stats(rp_diffs, sigma=n_sigma)
-            mean_of_means, _, std_of_means = sigma_clipped_stats(rp_means, sigma=n_sigma)
-            mean_of_stds, _, std_of_stds = sigma_clipped_stats(rp_stds, sigma=n_sigma)
-
-            # find the additional intermittent bad pixels, marking both readouts
-            high_diffs = (rp_diffs - mean_of_diffs) > (n_sigma * std_of_diffs)
-            high_means = (rp_means - mean_of_means) > (n_sigma * std_of_means)
-            high_stds = (rp_stds - mean_of_stds) > (n_sigma * std_of_stds)
-
-            log.debug(f'High diffs={np.sum(high_diffs)}, '
-                      f'high means={np.sum(high_means)}, '
-                      f'high stds={np.sum(high_stds)}')
-
-            int_bad[ref_pix[high_diffs]] = True
-            int_bad[pair_pixel[high_diffs]] = True
-            int_bad[ref_pix[high_means]] = True
-            int_bad[pair_pixel[high_means]] = True
-            int_bad[ref_pix[high_stds]] = True
-            int_bad[pair_pixel[high_stds]] = True
-
-            log.debug(f'{np.sum(int_bad[offset:offset + amplifier])} '
-                      f'suspicious bad reference pixels in '
-                      f'amplifier {k}, integration {j}')
-            mask_bad |= int_bad
+            if not replace_only:
+                ref_pix = np.array(ref_pix, dtype=int)
+                rp_diffs = np.array(rp_diffs)
+                rp_means = np.array(rp_means)
+                rp_stds = np.array(rp_stds)
+                pair_pixel = ref_pix + 1
+
+                # clipped stats for all tests
+                mean_of_diffs, _, std_of_diffs = sigma_clipped_stats(rp_diffs, sigma=n_sigma)
+                mean_of_means, _, std_of_means = sigma_clipped_stats(rp_means, sigma=n_sigma)
+                mean_of_stds, _, std_of_stds = sigma_clipped_stats(rp_stds, sigma=n_sigma)
+
+                # find the additional intermittent bad pixels, marking both readouts
+                high_diffs = (rp_diffs - mean_of_diffs) > (n_sigma * std_of_diffs)
+                high_means = (rp_means - mean_of_means) > (n_sigma * std_of_means)
+                high_stds = (rp_stds - mean_of_stds) > (n_sigma * std_of_stds)
+
+                log.debug(f'High diffs={np.sum(high_diffs)}, '
+                          f'high means={np.sum(high_means)}, '
+                          f'high stds={np.sum(high_stds)}')
+                int_bad[ref_pix[high_diffs]] = True
+                int_bad[pair_pixel[high_diffs]] = True
+                int_bad[ref_pix[high_means]] = True
+                int_bad[pair_pixel[high_means]] = True
+                int_bad[ref_pix[high_stds]] = True
+                int_bad[pair_pixel[high_stds]] = True
+
+                log.debug(f'{np.sum(int_bad[offset:offset + amplifier])} '
+                          f'suspicious bad reference pixels in '
+                          f'amplifier {k}, integration {j}')
+                mask_bad |= int_bad
 
         # replace any flagged pixels if desired
         if not flag_only:
@@ -649,7 +649,10 @@ def flag_bad_refpix(datamodel, n_sigma=3.0, flag_only=False, replace_only=False)
                                offset, offset + amplifier,
                                scipix_n, refpix_r)
 
-    log.info(f'Total bad reference pixels: {np.sum(mask_bad)}')
+    if flag_only:
+        log.info(f'Total bad reference pixels flagged: {np.sum(mask_bad)}')
+    else:
+        log.info(f'Total bad reference pixels replaced: {np.sum(mask_bad)}')
     if pixeldq is not None:
         pixeldq[mask_bad] |= (dqflags.pixel['BAD_REF_PIXEL']
                               | dqflags.pixel['DO_NOT_USE'])