JP-3157: Removal of NRS IRS2 intermittently bad pixels (#7685)

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

CHANGES.rst

@@ -24,6 +24,12 @@ jump
 - Added a test to prevent a divide by zero when the numger of usable
   groups is less than one. [#7723]
 
+refpix
+------
+
+- Replace intermittently bad pixels with nearest good reference pixel
+  for NIRSpec IRS2 mode. [#7685]
+
 tweakreg
 --------
 

docs/jwst/refpix/description.rst

@@ -194,6 +194,15 @@ 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)
+reference pixels. This is done by calculating the means and standard
+deviations per reference pixel column, as well as the difference between
+even and odd pairs; then calculates the mean of each of these arrays (the
+mean of the absolute values for the differences array), and flag all
+values greater than the corresponding mean. All suspicious pixels will
+be replaced by their nearest good reference pixel.
+
 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

@@ -135,6 +135,8 @@ def correct_model(input_model, irs2_model,
         # Set interleaved reference pixel values to zero if they are flagged
         # as bad in the DQ extension of the CRDS reference file.
         clobber_ref(data, output, odd_even, mask)
+        # Replace intermittently bad pixels
+        rm_intermittent_badpix(data, scipix_n, refpix_r)
     else:
         log.warning("DQ extension not found in reference file")
 
@@ -370,11 +372,8 @@ def decode_mask(output, mask):
     """
 
     # The bit number corresponds to a count of groups of reads of the
-    # interleaved reference pixels.  It wasn't clear to us whether bit
-    # number increases from left to right or right to left within a
-    # 32-bit unsigned integer.  It also wasn't clear whether the direction
-    # should follow the direction of reading within an amplifier output.
-    # The following is our interpretation of the description.
+    # interleaved reference pixels. The 32-bit unsigned integer encoding
+    # has increasing index, following the amplifier readout direction.
 
     flags = np.array([2**n for n in range(32)], dtype=np.uint32)
     temp = np.bitwise_and(flags, mask)
@@ -387,6 +386,118 @@ def decode_mask(output, mask):
     return bits
 
 
+def rm_intermittent_badpix(data, scipix_n, refpix_r):
+    """Find pixels that are intermittently bad and replace with nearest
+    good value.
+
+    Parameters
+    ----------
+    data : 4-D ndarray
+        The data array in detector orientation.  This includes both the
+        science and interleaved reference pixel values.  `data` will be
+        modified in-place. The science data values will not be modified.
+
+    scipix_n : int
+        Number of regular (science) samples before stepping out to collect
+        reference samples.
+
+    refpix_r : int
+        Number of reference samples before stepping back in to collect
+        regular samples.
+
+    Returns
+    -------
+    data : 4-D ndarray
+        The data array in detector orientation.  This includes both the
+        science and interleaved reference pixel values.  The intermittently
+        bad pixels are now set to the nearest good reference pixel value.
+    """
+    # The intermittently bad pixels will be replaced for all integrations
+    # and all groups. The last group will be used to identify them
+    nints, ngroups, ny, nx = np.shape(data)
+    total_rp2replace = []
+    # calculate differences of even and odd pairs per amplifier
+    amplifier = nx // 5                 # 640
+    for k in range(5):
+        diffs, rp2check = [], []
+        offset = int(k * amplifier)
+        # jump from the start of the reference pixel sequence to the next
+        ref_pix, rp2check, pair = [], [], 0
+        rp_means, rp_stds = [], []
+        for rpstart in range(scipix_n//2, amplifier, scipix_n+refpix_r):
+            rpstart += offset
+            # go through the 4 reference pixels
+            for ri in range(4):
+                ri = rpstart + ri
+                rp_m = np.mean(data[nints-1, ngroups-1, :, ri])
+                rp_s = np.std(data[nints-1, ngroups-1, :, ri])
+                # exclude ref pix flagged in the reference file
+                if rp_m != 0.:
+                    ref_pix.append(ri)
+                    rp_means.append(rp_m)
+                    rp_stds.append(rp_s)
+                if ri % 2 != 0:
+                    odd_pix = ri
+                    rp_odd = rp_m
+                else:
+                    even_pix = ri
+                    rp_even = rp_m
+                pair += 1
+                if pair == 2:
+                    # only do difference if the pixel has not already been flagged as bad
+                    if rp_odd != 0. and rp_even != 0.:
+                        diffs.append(rp_odd - rp_even)
+                        rp2check.append(even_pix)
+                        rp2check.append(odd_pix)
+                    pair = 0
+        diff_m = np.mean(np.abs(diffs))
+        mean_mean = np.mean(rp_means)
+        mean_std = np.mean(rp_stds)
+
+        # order indexes increasing from left to right
+        rp2check.sort()
+
+        # find the additional intermittent bad pixels
+        high_diffs = np.where(np.abs(diffs) > diff_m)[0]
+        hd_rp2replace = []
+        for j in high_diffs:
+            rp2r = rp2check[int(diffs.index(diffs[j]) * 2)]
+            # include both even and odd
+            hd_rp2replace.append(rp2r)
+            hd_rp2replace.append(rp2r+1)
+        high_means_idx = np.where(np.array(rp_means) > mean_mean)[0]
+        high_std_idx = np.where(np.array(rp_stds) > mean_std)[0]
+        ref_pix = np.array(ref_pix)
+        rp2replace = []
+        for rp in ref_pix:
+            if rp in ref_pix[high_means_idx] or rp in ref_pix[high_std_idx] or rp in hd_rp2replace:
+                rp2replace.append(rp)
+        rp2replace.sort()
+        rp2replace = np.unique(rp2replace)
+        total_rp2replace.extend(rp2replace)
+        log.info('{} supicious bad reference pixels in amplifier {}'.format(len(rp2replace), k))
+
+        remaining_rp_even, remaining_rp_odd = [], []
+        for rp in ref_pix:
+            if rp not in rp2replace:
+                if rp % 2 == 0:
+                    remaining_rp_even.append(rp)
+                else:
+                    remaining_rp_odd.append(rp)
+        remaining_rp_even, remaining_rp_odd = np.array(remaining_rp_even), np.array(remaining_rp_odd)
+        for bad_pix in rp2replace:
+            # find nearest even/odd good reference pixel
+            if bad_pix % 2 == 0:
+                remaining_rp = remaining_rp_even
+            else:
+                remaining_rp = remaining_rp_odd
+            good_idx = (np.abs(remaining_rp - bad_pix)).argmin()
+            good_pix = remaining_rp[good_idx]
+            data[..., bad_pix] = data[..., good_pix]
+            log.debug('   Pixel {}'.format(bad_pix))
+    log.info('Total intermittent bad reference pixels: {}'.format(len(total_rp2replace)))
+
+
 def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
     """Subtract reference output and pixels for the current integration.
 
@@ -761,7 +872,6 @@ def replace_bad_pixels(data0, ngroups, ny, row):
         dat = numerator / denominator
         dat = dat.reshape(ny, row)
         mask = mask.reshape(ny, row)
-        # xxx why '+=' instead of just '=' ?
         data0[kk, jj, :, :] += dat * (1. - mask)
 
 
@@ -801,7 +911,7 @@ def fill_bad_regions(data0, ngroups, ny, nx, row, scipix_n, refpix_r, pad, hnorm
 
     # IDL:  aa = a # replicate(1, s[3]) ; for application to the data
     # In IDL, aa is a 2-D array with one column of `a` for each group.  In
-    # Python, numpy broadcasting should take care of this.
+    # Python, numpy broadcasting takes care of this.
 
     n_iter_norm = 3
     dd0 = data0[0, :, :, :]

jwst/refpix/tests/test_rm_intermittent_badpix.py

@@ -0,0 +1,25 @@
+import numpy as np
+
+from jwst.refpix.irs2_subtract_reference import rm_intermittent_badpix
+
+
+def test_rm_intermittent_badpix():
+    data = np.ones((2, 3, 5, 3200), dtype=np.float32)
+    # set the bad reference pixels flagged in the ref file
+    data[..., 688] = 0.
+    data[..., 2110] = 0.
+    # set the intermittent bad pixels
+    data[..., 648] = 10.
+    data[..., 988] = 11.
+    data[..., 1369] = 7.
+    data[..., 2150] = 13.
+    data[..., 3128] = 15.
+
+    scipix_n, refpix_r = 16, 4
+    rm_intermittent_badpix(data, scipix_n, refpix_r)
+
+    compare = np.ones((2, 3, 5, 3200), dtype=np.float32)
+    compare[..., 688] = 0.
+    compare[..., 2110] = 0.
+
+    assert np.allclose(data, compare)