Set DQ of ref pixels to DO_NOT_USE after refpix correction

docs/jwst/refpix/description.rst

@@ -77,6 +77,8 @@ NIR Detector Data
        subtracted from the full group on a row-by-row basis.  Note that the ``odd_even_rows``
        parameter is ignored for NIR data when the side reference pixels are processed.
     #. Transform the data back to the JWST focal plane, or DMS, frame.
+    #. Flag reference pixel locations in DQ array as DO_NOT_USE now that the correction is
+       complete, so that they are not used in ramp fitting or any other downstream steps.
 
 MIR Detector Data
 +++++++++++++++++

jwst/refpix/reference_pixels.py

@@ -47,10 +47,13 @@
 from ..datamodels import dqflags
 from ..lib import reffile_utils
 
+
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
 
-#
+DO_NOT_USE = dqflags.pixel["DO_NOT_USE"]
+REFERENCE_PIXEL = dqflags.pixel["REFERENCE_PIXEL"]
+
 # NIR Reference section dictionaries are zero indexed and specify the values
 # to be used in the following slice:
 # (rowstart: rowstop, colstart:colstop)
@@ -216,15 +219,14 @@ def sigma_clip(self, data, dq, low=3.0, high=3.0):
 
         """
 
-        #
         # Only calculate the clipped mean for pixels that don't have the DO_NOT_USE
-        # DQ bit set
-        goodpixels = np.where(np.bitwise_and(dq, dqflags.pixel['DO_NOT_USE']) == 0)
-        #
+        # DQ bit set.  Create a boolean mask.
+        goodpixels = np.bitwise_and(dq, DO_NOT_USE) != DO_NOT_USE
+
         # If there are no good pixels, return None
-        if len(goodpixels[0]) == 0:
+        if not goodpixels.any():
             return None
-        #
+
         # scipy routine fails if the pixels all have exactly the same value
         if np.std(data[goodpixels], dtype=np.float64) != 0.0:
             clipped_ref, lowlim, uplim = stats.sigmaclip(data[goodpixels],
@@ -261,11 +263,10 @@ def get_pixeldq(self):
                 fullcols = 2048
             self.full_shape = (fullrows, fullcols)
             pixeldq = np.zeros(self.full_shape, dtype=self.input_model.pixeldq.dtype)
-            refpixdq_dontuse = dqflags.pixel['DO_NOT_USE'] | dqflags.pixel['REFERENCE_PIXEL']
-            pixeldq[0:4, :] = refpixdq_dontuse
-            pixeldq[fullrows - 4:fullrows, :] = refpixdq_dontuse
-            pixeldq[4:fullrows - 4, 0:4] = refpixdq_dontuse
-            pixeldq[4:fullrows - 4, fullcols - 4:fullcols] = refpixdq_dontuse
+            pixeldq[0:4, :] = DO_NOT_USE | REFERENCE_PIXEL
+            pixeldq[fullrows - 4:fullrows, :] = DO_NOT_USE | REFERENCE_PIXEL
+            pixeldq[4:fullrows - 4, 0:4] = DO_NOT_USE | REFERENCE_PIXEL
+            pixeldq[4:fullrows - 4, fullcols - 4:fullcols] = DO_NOT_USE | REFERENCE_PIXEL
             pixeldq[self.rowstart:self.rowstop, self.colstart:self.colstop] = self.input_model.pixeldq.copy()
         else:
             pixeldq = self.input_model.pixeldq.copy()
@@ -394,21 +395,19 @@ def log_parameters(self):
                 log.info('refpix processing skipped for this mode')
 
     def count_good_side_refpixels(self):
-        donotuse = dqflags.pixel['DO_NOT_USE']
         ngood = 0
         for amplifier in 'AD':
             rowstart, rowstop, colstart, colstop = NIR_reference_sections[amplifier]['side']
-            good = np.where(np.bitwise_and(self.pixeldq[rowstart:rowstop, colstart:colstop], donotuse) != donotuse)
+            good = np.where(np.bitwise_and(self.pixeldq[rowstart:rowstop, colstart:colstop], DO_NOT_USE) != DO_NOT_USE)
             ngood += len(good[0])
         return ngood
 
     def count_good_top_bottom_refpixels(self):
-        donotuse = dqflags.pixel['DO_NOT_USE']
         ngood = 0
         for edge in ['top', 'bottom']:
             for amplifier in 'ABCD':
                 rowstart, rowstop, colstart, colstop = NIR_reference_sections[amplifier][edge]
-                good = np.where(np.bitwise_and(self.pixeldq[rowstart:rowstop, colstart:colstop], donotuse) != donotuse)
+                good = np.where(np.bitwise_and(self.pixeldq[rowstart:rowstop, colstart:colstop], DO_NOT_USE) != DO_NOT_USE)
                 ngood += len(good[0])
         return ngood
 
@@ -825,7 +824,7 @@ def median_filter(self, data, dq, smoothing_length):
             rowstart = i
             rowstop = rowstart + smoothing_length
             goodpixels = np.where(np.bitwise_and(augmented_dq[rowstart:rowstop],
-                                                 dqflags.pixel['DO_NOT_USE']) == 0)
+                                                 DO_NOT_USE) == 0)
             if len(goodpixels[0]) == 0:
                 result[i] = np.nan
             else:
@@ -1012,17 +1011,14 @@ def do_subarray_corrections(self):
         """Do corrections for subarray.  Reference pixel value calculated
         separately for odd and even columns if odd_even_columns is True,
         otherwise a single number calculated from all reference pixels"""
-        #
+
         #  First transform to detector coordinates
         #
-        refdq = dqflags.pixel['REFERENCE_PIXEL']
-        donotuse = dqflags.pixel['DO_NOT_USE']
-        #
         # This transforms the pixeldq array from DMS to detector coordinates,
         # only needs to be done once
         self.DMS_to_detector_dq()
         # Determined refpix indices to use on each group
-        refpixindices = np.where((self.pixeldq & refdq == refdq) & (self.pixeldq & donotuse != donotuse))
+        refpixindices = np.where((self.pixeldq & REFERENCE_PIXEL == REFERENCE_PIXEL) & (self.pixeldq & DO_NOT_USE != DO_NOT_USE))
         nrefpixels = len(refpixindices[0])
         if nrefpixels == 0:
             self.bad_reference_pixels = True
@@ -1931,6 +1927,10 @@ def correct_model(input_model, odd_even_columns,
     input_dataset.log_parameters()
     reference_pixel_correction(input_dataset)
 
+    # Now that they are corrected, flag reference pixels as DO_NOT_USE
+    mask = input_model.pixeldq & REFERENCE_PIXEL == REFERENCE_PIXEL
+    input_model.pixeldq[mask] |= DO_NOT_USE
+
     return REFPIX_OK
 
 
@@ -2044,7 +2044,7 @@ def setup_dataset_for_zeroframe(input_dataset, saved_values):
     gdtype = input_dataset.input_model.groupdq.dtype
     gdq = np.zeros(dims, dtype=gdtype)
     wh_zero = saved_values[-1]
-    gdq[wh_zero] = dqflags.pixel['DO_NOT_USE']
+    gdq[wh_zero] = DO_NOT_USE
     gdq = gdq.reshape(new_dims)
 
     # Setup dataset with ZEROFRAME data

jwst/refpix/tests/test_refpix.py

@@ -6,6 +6,10 @@
 from jwst.refpix.reference_pixels import Dataset, NIRDataset, correct_model, create_dataset
 
 
+REFERENCE_PIXEL = dqflags.pixel["REFERENCE_PIXEL"]
+DO_NOT_USE = dqflags.pixel["DO_NOT_USE"]
+
+
 def test_refpix_subarray():
     '''Check that the correction is skipped for MIR subarray data '''
 
@@ -59,8 +63,8 @@ def test_each_amp():
     im.data[:, 1:, :, 1031] = 4.0
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
 
     # run the step
     out = RefPixStep.call(im)
@@ -103,8 +107,8 @@ def test_firstframe_sub():
     im.data[:, :, :, 1031] = 4.0
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
 
     # run the step
     outim = RefPixStep.call(im)
@@ -148,8 +152,8 @@ def test_odd_even():
     im.data[:, 1:, 0:ysize - 1:2, 1031] = 8.0
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
 
     # run the step
     out = RefPixStep.call(im)
@@ -165,6 +169,11 @@ def test_odd_even():
     assert out.data[0, 5, 101, 6] == 47.0
     assert out.data[0, 5, 101, 7] == 46.0
 
+    # Make sure reference pixel DQs are set to DO_NOT_USE
+    mask1 = out.pixeldq & REFERENCE_PIXEL == REFERENCE_PIXEL
+    mask2 = out.pixeldq & DO_NOT_USE == DO_NOT_USE
+    np.testing.assert_array_equal(mask1, mask1 & mask2)
+
 
 def test_no_odd_even():
     '''Check that odd/even rows are not applied if flag is set to False'''
@@ -200,8 +209,8 @@ def test_no_odd_even():
     im.data[:, 1:, 0:ysize - 1:2, 1031] = 8.0
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
 
     # run the step
     out = RefPixStep.call(im, odd_even_rows=False)
@@ -240,8 +249,8 @@ def test_side_averaging():
     im.data[:, 1:, :, 1028:] = 2.0
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
 
     # run the step
     out = RefPixStep.call(im)
@@ -271,8 +280,8 @@ def test_above_sigma():
     im.data[0, 3, 50, 3] = 35.0
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
 
     # run the step
     out = RefPixStep.call(im)
@@ -305,9 +314,9 @@ def test_nan_refpix():
     im.data[0, 3, 50, 3] = np.nan
 
     # set reference pixels to 'REFERENCE_PIXEL'
-    im.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[:, 1028:] = dqflags.pixel['REFERENCE_PIXEL']
-    im.pixeldq[50, 3] = dqflags.pixel['DO_NOT_USE']
+    im.pixeldq[:, :4] = REFERENCE_PIXEL
+    im.pixeldq[:, 1028:] = REFERENCE_PIXEL
+    im.pixeldq[50, 3] = DO_NOT_USE
 
     # run the step
     out = RefPixStep.call(im)
@@ -341,8 +350,8 @@ def test_do_corrections_subarray_no_oddEven(setup_subarray_cube):
     input_model.data[0, 0, :, :] = dataval
     input_model.data[0, 0, :4, :] = bottom_rpix
     input_model.data[0, 0, :, :4] = left_rpix
-    input_model.pixeldq[:4, :] = dqflags.pixel['REFERENCE_PIXEL']
-    input_model.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
+    input_model.pixeldq[:4, :] = REFERENCE_PIXEL
+    input_model.pixeldq[:, :4] = REFERENCE_PIXEL
 
     init_dataset = create_dataset(input_model,
                                   odd_even_columns,
@@ -383,8 +392,8 @@ def test_do_corrections_subarray(setup_subarray_cube):
     input_model.data[0, 0, :, :] = dataval
     input_model.data[0, 0, :4, :] = bottom_rpix
     input_model.data[0, 0, :, :4] = left_rpix
-    input_model.pixeldq[:4, :] = dqflags.pixel['REFERENCE_PIXEL']
-    input_model.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
+    input_model.pixeldq[:4, :] = REFERENCE_PIXEL
+    input_model.pixeldq[:, :4] = REFERENCE_PIXEL
 
     init_dataset = create_dataset(input_model,
                                   odd_even_columns,
@@ -453,9 +462,9 @@ def test_do_corrections_subarray_4amp(setup_subarray_cube):
     input_model.data[0, 0, :, 1:4:2] = left_rpix + bottom_rpix_a_even
     input_model.data[0, 0, :, -4::2] = right_rpix + bottom_rpix_d_odd
     input_model.data[0, 0, :, -3::2] = right_rpix + bottom_rpix_d_even
-    input_model.pixeldq[:4, :] = dqflags.pixel['REFERENCE_PIXEL']
-    input_model.pixeldq[:, :4] = dqflags.pixel['REFERENCE_PIXEL']
-    input_model.pixeldq[:, -4:] = dqflags.pixel['REFERENCE_PIXEL']
+    input_model.pixeldq[:4, :] = REFERENCE_PIXEL
+    input_model.pixeldq[:, :4] = REFERENCE_PIXEL
+    input_model.pixeldq[:, -4:] = REFERENCE_PIXEL
 
     init_dataset = create_dataset(input_model,
                                   odd_even_columns,
@@ -752,6 +761,12 @@ def test_correct_model(setup_cube, instr, det):
     input_model.data[0, 0, :, :] = rpix
     input_model.data[0, 0, 4:-4, 4:-4] = dataval
 
+    # Populate DQ array with REFERENCE_PIXEL where appropriate
+    input_model.pixeldq[:4,:] = REFERENCE_PIXEL
+    input_model.pixeldq[-4:,:] = REFERENCE_PIXEL
+    input_model.pixeldq[:,:4] = REFERENCE_PIXEL
+    input_model.pixeldq[:,-4:] = REFERENCE_PIXEL
+
     correct_model(input_model,
                   odd_even_columns,
                   use_side_ref_pixels,
@@ -762,6 +777,11 @@ def test_correct_model(setup_cube, instr, det):
     np.testing.assert_almost_equal(np.mean(input_model.data[0, 0, :4, 4:-4]), 0, decimal=0)
     np.testing.assert_almost_equal(np.mean(input_model.data[0, 0, 4:-4, 4:-4]), dataval - rpix, decimal=0)
 
+    # Make sure reference pixel DQs are set to DO_NOT_USE
+    mask1 = input_model.pixeldq & REFERENCE_PIXEL == REFERENCE_PIXEL
+    mask2 = input_model.pixeldq & DO_NOT_USE == DO_NOT_USE
+    np.testing.assert_array_equal(mask1, mask1 & mask2)
+
 
 def test_zero_frame(setup_cube):
     """