[Bug]: Saturation flags overwritten with do_not_use flags

[ben-cassese]

FAQ check

• Yes, I checked the FAQ and my question has not been addressed.

Instrument

NIRSpec (Stages 1-3)

What happened?

Hello! Thanks to all the developers for their work on this package, this has been a great resource for getting started with JWST data. I've labeled this a "bug", but am not 100% sure that's the case- this might be an intentional design choice, but either way, some clarification about what's happening would be appreciated.

At the end of update_sat in src/eureka/S1_detector_processing/update_saturation.py, we modify the groupdq attribute of the datamodel to reflect the more aggressive saturation flags computed earlier in the function. I think the idea is to add a SATURATED flag for the more aggressive saturation masking, and a DO_NOT_USE flag for any column that contains a pixel that saturates on the first group. However, since in this line we change the shape of new_sat_mask back to (n_integrations, n_groups, n_rows, n_columns), the slice new_sat_mask[0,:,:] on this line no longer checks the first group for saturation, but the first integration. This is in contrast to the slicing a few lines above when checking whether to log the existence of any first-group saturated pixels (if np.count_nonzero(new_sat_mask[0]) > 0:), since this is before the shape change.

The upshot is that any pixel that's marked as saturated ends up flagged as DO_NOT_USE, not as SATURATED. Is this intentional? It looks like there was some discussion in #597 and #651, though I think the goal there was to mark only the first-group saturated pixels/columns as DO_NOT_USE.

Here's a MRE, using the first segment of ERS NIRSpec Prism data of WASP-39b, jw01366004001_04101_00001-seg001_nrs1_uncal.fits:

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits

import stdatamodels.jwst.datamodels as dm
from stdatamodels.jwst.datamodels import dqflags
from jwst.group_scale import GroupScaleStep
from jwst.dq_init import DQInitStep
from jwst.saturation import SaturationStep

f = "jw01366004001_04101_00001-seg001_nrs1_uncal.fits"
model = dm.open(f)
model = GroupScaleStep.call(model)
model = DQInitStep.call(model)
model = SaturationStep.call(model)

# this update_sat is same as eureka.S1_detector_processing.update_saturation,
# just with logger and meta stripped out
updated_sat_model = update_sat(model)

fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(8, 8))
ax[0, 0].imshow(updated_sat_model.groupdq[1000, 0, :, :], vmin=0, vmax=2)
ax[0, 1].imshow(model.groupdq[1000, 0, :, :], vmin=0, vmax=2)
ax[1, 0].imshow(updated_sat_model.groupdq[1000, 1, :, :], vmin=0, vmax=2)
ax[1, 1].imshow(model.groupdq[1000, 1, :, :], vmin=0, vmax=2)
ax[2, 0].imshow(updated_sat_model.groupdq[1000, 2, :, :], vmin=0, vmax=2)
ax[2, 1].imshow(model.groupdq[1000, 2, :, :], vmin=0, vmax=2)

for i in range(3):
    for j in range(2):
        ax[i, j].set(aspect=8)
for i in range(3):
    ax[i, 0].set(title=f"updated saturation, group {i+1}")
    ax[i, 1].set(title=f"jwst default, group {i+1}")

The DO_NOT_USE flag has a value of 1, while SATURATED has a value of 2. Even though no pixels are saturated on the first group here (the SaturationStep of jwst says nothing saturates until group 3, so we would want to start flagging group 2 as suspect), a large number of pixels end up being marked DO_NOT_USE and will end up ignored later in processing.

Let me know if this was the intended behavior, or if not, if I should open a PR to change the definition of full_saturation. Thanks again!

Error traceback output

No errors.

What operating system are you using?

macOS Sequoia 15.1

What version of Python are you running?

3.10.14

What Python packages do you have installed?

No response

Code of Conduct

• I agree to follow this project's Code of Conduct

[kevin218]

@erinmmay, can you verify that update_saturation.py achieves the intended behavior? See details above.