Repository: spacetelescope/jwst

Issue:
Add SCI/ERR/DQ check to resample/cube_build 
_Issue [JP-3570](https://jira.stsci.edu/browse/JP-3570) was created on JIRA by [David Law](https://jira.stsci.edu/secure/ViewProfile.jspa?name=dlaw):_

As seen in tickets such as <https://jira.stsci.edu/browse/JP-3250],> unexpected NaNs in the ERR or VAR arrays can cause problems downstream since those NaNs can be propagated to the resampled products.

This ticket is to therefore catch such issues earlier, by adding a check to the start of resample and cube_build to ensure that all input data with SCI=NaN also have ERR=NaN and DQ DO_NOT_USE bit set.  Likewise, ERR=NaN should have SCI=NaN and DQ DO_NOT_USE bit set.  I believe cube_build may already have something like this implemented, as it is fairly robust against NaN-valued inputs.  An equivalent may also be needed elsewhere in TSO3 processing since that does not use the resample step.

Once implemented, this will need some fairly thorough instrument-mode testing prior to being merged as it may reveal issues with reference files that need to be addressed first.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
Attention: Patch coverage is `87.34177%` with `10 lines` in your changes missing coverage. Please review.
> Project coverage is 58.61%. Comparing base [(`b7e0b10`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/b7e0b107d54ee411b252518e149bbf5501e10403?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to head [(`f78fccf`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/f78fccfeec99d7bcd09cc5a7f19d6c584a776949?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> Report is 10 commits behind head on master.

| [Files](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Patch % | Lines |
|---|---|---|
| [jwst/pathloss/pathloss.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&filepath=jwst%2Fpathloss%2Fpathloss.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9wYXRobG9zcy9wYXRobG9zcy5weQ==) | 16.66% | [5 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/flatfield/flat\_field.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&filepath=jwst%2Fflatfield%2Fflat_field.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9mbGF0ZmllbGQvZmxhdF9maWVsZC5weQ==) | 40.00% | [3 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/outlier\_detection/outlier\_detection.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&filepath=jwst%2Foutlier_detection%2Foutlier_detection.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbi5weQ==) | 75.00% | [1 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |
| [jwst/outlier\_detection/outlier\_detection\_ifu.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&filepath=jwst%2Foutlier_detection%2Foutlier_detection_ifu.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9vdXRsaWVyX2RldGVjdGlvbi9vdXRsaWVyX2RldGVjdGlvbl9pZnUucHk=) | 50.00% | [1 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #8557      +/-   ##
==========================================
+ Coverage   58.02%   58.61%   +0.58%     
==========================================
  Files         388      388              
  Lines       38977    38999      +22     
==========================================
+ Hits        22617    22858     +241     
+ Misses      16360    16141     -219     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/8557?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

melanieclarke: Started regression tests here:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1543/

I predict many new NaNs!
melanieclarke: Now that we are past build 11.0, I think this is ready for discussion and testing.

Regression test results:
- Many new NaNs in intermediate products in SCI/ERR/VAR arrays, most notably outputs from outlier_detection, flat_field, and photom.
- Imaging source catalogs show differences because NaNs are not ignored in photometry (e.g. [test_miri_image3_catalog](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1552/testReport/junit/jwst.regtest/test_miri_image/_stable_deps__test_miri_image3_catalog/))
- NIRSpec IFU shows changes in the DQ array after flat_field: NON_SCIENCE areas are now also flagged DO_NOT_USE.
- There are some small changes to values in the SURF_BRIGHT column only in x1d products (e.g. [test_nirspec_mos_spec3[v000000053-x1d]](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1552/testReport/junit/jwst.regtest/test_nirspec_mos_spec3/_stable_deps__test_nirspec_mos_spec3_v000000053_x1d_/))
- There are some small changes to weight values in s2d files (e.g. [test_nirspec_mos_spec3[v000000053-s2d]](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1552/testReport/junit/jwst.regtest/test_nirspec_mos_spec3/_stable_deps__test_nirspec_mos_spec3_v000000053_s2d_/))
- There are some changes to values in NIRCam coronography final products (e.g. [test_nircam_coron3_product[i2d]](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1552/testReport/junit/jwst.regtest/test_nircam_coron3/_stable_deps__test_nircam_coron3_product_i2d_/)).
melanieclarke: Thanks for reviewing @emolter.

> is the idea that the NaNs in the data and variance arrays are _always_ synced with the DO_NOT_USE flags in the DQ array?  If so, should we consider making some kind of check for this generic to all pipeline steps?  I guess it isn't obvious to me that the steps where you added this are the major/only places this would be helpful.

I think that's the idea, but I was reluctant to add it to every single pipeline step individually in the first pass.  I added it to steps that I know either modify DO_NOT_USE or set NaNs in data arrays.  I can add it to more if needed.

I didn't consider whether it could or should be done more generically, in a core function or class.  It does need some knowledge of the output model composition, so that it does the updates at the right level.  Do you have a better sense of where that might go?

> * I'm glad you checked edge cases where certain extensions are not found - not sure I would have thought of that - but why is the preferred behavior in those cases to return a copy of the datamodel instead of failing?

I was thinking that it should be as minimally invasive as possible, so it can be called on all kinds of datamodels, without assuming that they have all the extensions at all times.  As is, it could be used on partial copies of datamodels as well as on the final products that should have all the expected extensions.
emolter: > I didn't consider whether it could or should be done more generically, in a core function or class. It does need some knowledge of the output model composition, so that it does the updates at the right level. Do you have a better sense of where that might go?

I don't know if I have the expertise or experience to answer the _should_ question, but I am sure it _could_ be implemented. 
 If that were the case I'd suggest putting it somewhere into `jwst.stpipe`, maybe `jwst.stpipe.utils`, and then calling it from the core JWST Step class.  It reminds me a bit of the `query_step_status` function we just implemented, and with regards to knowledge of the model composition, is it anything more complicated than what is implemented in that `query_step_status` function?  That is, do something to the model if it's not a Sequence, loop over all models in the sequence if it is?
melanieclarke: Per @drlaw1558 - we should keep the implementation limited to these specific steps for now, to avoid scope creep and possible unintended consequences.
melanieclarke: Thanks, @drlaw1558 for coordinating the testing!  I'll rebase and fix up the conflicts with later development and see if we can get this in.
melanieclarke: Outlier detection changed significantly after I initially wrote this PR, so I'm re-requesting review from @emolter and @braingram, if you have a chance.

I'll re-run regression tests shortly.
melanieclarke: Regression test results are inconclusive.  There are many new failures introduced by this branch, as expected, but it looks like there are also some unrelated failures that need okifying.  I will run again when things are more stable.
melanieclarke: I think the regression test output is as expected, compared to the previous run.

One thing I noticed in addition to the differences noted above: i2d and crf output for image3 is sometimes different (e.g. [test_miri_image3_i2d](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1552/testReport/junit/jwst.regtest/test_miri_image/_stable_deps__test_miri_image3_i2d/)) because tweakreg is seeing a different set of sources to cross match on.  It might be worth looking into how the starfinder is handling NaNs: I wouldn’t have expected to see much difference here, since the DO_NOT_USE mask is already passed.
melanieclarke: > It looks like `match_nans_and_flags` will not get called during the outlier detection step if `resample_data` is `False` for imaging or spectroscopic modes. Is that true, and if so, is that the intended behavior?

It should still be called in that case, via `flag_model_crs`.  I'll test.
melanieclarke: @emolter - confirmed, it's definitely called for image and spec outlier detection when resample_data is False.

@drlaw1558 - the new ref files are definitely confusing things in comparing to the current truth files, but I think there is additionally a change in tweakreg due to this branch. I ran the regtest locally with the same CRDS context on the current main and PR branches and compared the log output.

On main, tweakreg reports:
```
Detected 171 sources in jw01024001001_04101_00002_mirimage_cal.fits.
Detected 175 sources in jw01024001001_04101_00003_mirimage_cal.fits.
Detected 167 sources in jw01024001001_04101_00004_mirimage_cal.fits.
Detected 150 sources in jw01024001001_04101_00001_mirimage_cal.fits.
...
Final solution based on 111 objects.
```

On this branch, tweakreg reports:
```
Detected 155 sources in jw01024001001_04101_00002_mirimage_cal.fits.
Detected 161 sources in jw01024001001_04101_00003_mirimage_cal.fits.
Detected 153 sources in jw01024001001_04101_00004_mirimage_cal.fits.
Detected 149 sources in jw01024001001_04101_00001_mirimage_cal.fits.
...
Final solution based on 109 objects.
```

It's not a big change, and I don't think it should hold up this PR, but we might want to follow up on it later to make sure the star finder isn't handling NaNs badly.
melanieclarke: Regression tests one more time, to test the changes with the ModelLibrary work:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1689/
melanieclarke: Cross matching regression tests to the previous run, I don't see anything new or unexpected, so I think this is ready to go.

@tapastro, over to you!

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 1c8995df2..661f78510 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -40,6 +40,9 @@ cube_build
   
 - Fixed a bug when ``cube_build`` was called from the ``mrs_imatch`` step. [#8728]
 
+- Ensure that NaNs and DO_NOT_USE flags match up in all input data before
+  building a cube. [#8557]
+
 datamodels
 ----------
 
@@ -57,6 +60,12 @@ emicorr
 - Fixed a bug where MIRI EMI correction step would return NaNs when it was unable
   to compute a correction. [#8675]
 
+flat_field
+----------
+
+- Ensure that NaNs and DO_NOT_USE flags match up in all science, error,
+  variance, and DQ extensions for all modes. [#8557]
+
 general
 -------
 
@@ -99,10 +108,20 @@ outlier_detection
 - For imaging modes, step now uses `ModelLibrary` to handle accessing models consistently
   whether they are in memory or on disk. [#8683]
 
-set_telescope_pointing
-----------------------
+- Ensure that NaNs and DO_NOT_USE flags match up in all output science, error,
+  variance, and DQ extensions. [#8557]
+
+pathloss
+--------
+
+- Ensure that NaNs and DO_NOT_USE flags match up in all output science, error,
+  variance, and DQ extensions. [#8557]
+
+photom
+------
 
-- replace usage of ``copy_arrays=True`` with ``memmap=False`` [#8660]
+- Ensure that NaNs and DO_NOT_USE flags match up in all output science, error,
+  variance, and DQ extensions. [#8557]
 
 pipeline
 --------
@@ -114,9 +133,16 @@ pipeline
 resample
 --------
 
+- Fixed a typo in ``load_custom_wcs`` from ``array_shape`` to ``pixel_shape`` and
+  changed to use values in the top-level ASDF structure if the values in the WCS
+  are ``None``. [#8698]
+
 - Step now uses `ModelLibrary` to handle accessing models consistently
   whether they are in memory or on disk. [#8683]
 
+- Ensure that NaNs and DO_NOT_USE flags match up in all input data before
+  resampling. [#8557]
+
 resample_spec
 -------------
 
@@ -138,12 +164,8 @@ resample_spec
   so that the spectral resampling step only exposes parameters that are appropriate
   for spectral data. [#8622]
 
-resample_step
--------------
-
-- Fixed a typo in ``load_custom_wcs`` from ``array_shape`` to ``pixel_shape`` and
-  changed to use values in the top-level ASDF structure if the values in the WCS
-  are ``None``. [#8698]
+- Ensure that NaNs and DO_NOT_USE flags match up in all input data before
+  resampling. [#8557]
 
 scripts
 -------
@@ -151,6 +173,11 @@ scripts
 - Removed many non-working and out-dated scripts. Including
   many scripts that were replaced by ``strun``. [#8619]
 
+set_telescope_pointing
+----------------------
+
+- Replace usage of ``copy_arrays=True`` with ``memmap=False`` [#8660]
+
 skymatch
 --------
 
diff --git a/jwst/cube_build/cube_build_step.py b/jwst/cube_build/cube_build_step.py
index 7801169ba..c675ae338 100755
--- a/jwst/cube_build/cube_build_step.py
+++ b/jwst/cube_build/cube_build_step.py
@@ -1,14 +1,16 @@
 """ This is the main ifu spectral cube building routine.
 """
 
+import time
+
 from jwst.datamodels import ModelContainer
+from jwst.lib.pipe_utils import match_nans_and_flags
 
-from ..stpipe import Step, record_step_status
 from . import cube_build
 from . import ifu_cube
 from . import data_types
 from ..assign_wcs.util import update_s_region_keyword
-import time
+from ..stpipe import Step, record_step_status
 
 __all__ = ["CubeBuildStep"]
 
@@ -281,6 +283,10 @@ class CubeBuildStep (Step):
 # ________________________________________________________________________________
 # create an instance of class CubeData
 
+        # Make sure all input models have consistent NaN and DO_NOT_USE values
+        for model in self.input_models:
+            match_nans_and_flags(model)
+
         cubeinfo = cube_build.CubeData(
             self.input_models,
             par_filename,
diff --git a/jwst/flatfield/flat_field.py b/jwst/flatfield/flat_field.py
index 0b653d8d2..0ce1fcf80 100644
--- a/jwst/flatfield/flat_field.py
+++ b/jwst/flatfield/flat_field.py
@@ -10,7 +10,7 @@ import numpy as np
 from stdatamodels.jwst import datamodels
 from stdatamodels.jwst.datamodels import dqflags
 
-from ..lib import reffile_utils, wcs_utils
+from ..lib import pipe_utils, reffile_utils, wcs_utils
 from ..assign_wcs import nirspec
 
 log = logging.getLogger(__name__)
@@ -234,6 +234,9 @@ def apply_flat_field(science, flat, inverse=False):
     flag_nonsci = np.bitwise_and(science.dq, dqflags.pixel['NON_SCIENCE']).astype(bool)
     science.dq[flag_nonsci] = np.bitwise_or(science.dq[flag_nonsci], dqflags.pixel['DO_NOT_USE'])
 
+    # Make sure all NaNs and flags match up in the output model
+    pipe_utils.match_nans_and_flags(science)
+
 
 #
 # The following functions are for NIRSpec spectrographic data.
@@ -438,14 +441,8 @@ def nirspec_fs_msa(output_model, f_flat_model, s_flat_model, d_flat_model, dispa
         # Combine the science and flat DQ arrays
         slit.dq |= slit_flat.dq
 
-        # Make sure all DO_NOT_USE pixels are set to NaN,
-        # including those flagged by this step
-        dnu = np.where(slit.dq & dqflags.pixel['DO_NOT_USE'])
-        slit.data[dnu] = np.nan
-        slit.err[dnu] = np.nan
-        slit.var_poisson[dnu] = np.nan
-        slit.var_rnoise[dnu] = np.nan
-        slit.var_flat[dnu] = np.nan
+        # Make sure all NaNs and flags match up in the output model
+        pipe_utils.match_nans_and_flags(slit)
 
         any_updated = True
 
@@ -522,14 +519,8 @@ def nirspec_brightobj(output_model, f_flat_model, s_flat_model, d_flat_model, di
         output_model.var_poisson + output_model.var_rnoise + output_model.var_flat
     )
 
-    # Make sure all DO_NOT_USE pixels are set to NaN,
-    # including those flagged by this step
-    dnu = np.where(output_model.dq & dqflags.pixel['DO_NOT_USE'])
-    output_model.data[dnu] = np.nan
-    output_model.err[dnu] = np.nan
-    output_model.var_poisson[dnu] = np.nan
-    output_model.var_rnoise[dnu] = np.nan
-    output_model.var_flat[dnu] = np.nan
+    # Make sure all NaNs and flags match up in the output model
+    pipe_utils.match_nans_and_flags(output_model)
 
     output_model.meta.cal_step.flat_field = 'COMPLETE'
 
@@ -597,14 +588,8 @@ def nirspec_ifu(output_model, f_flat_model, s_flat_model, d_flat_model, dispaxis
             output_model.var_poisson + output_model.var_rnoise + output_model.var_flat
         )
 
-        # Make sure all DO_NOT_USE pixels are set to NaN,
-        # including those flagged by this step
-        dnu = np.where(output_model.dq & dqflags.pixel['DO_NOT_USE'])
-        output_model.data[dnu] = np.nan
-        output_model.err[dnu] = np.nan
-        output_model.var_poisson[dnu] = np.nan
-        output_model.var_rnoise[dnu] = np.nan
-        output_model.var_flat[dnu] = np.nan
+        # Make sure all NaNs and flags match up in the output model
+        pipe_utils.match_nans_and_flags(output_model)
 
         output_model.meta.cal_step.flat_field = 'COMPLETE'
 
diff --git a/jwst/lib/basic_utils.py b/jwst/lib/basic_utils.py
index 92e77dcf4..28971ddb8 100644
--- a/jwst/lib/basic_utils.py
+++ b/jwst/lib/basic_utils.py
@@ -1,31 +1,5 @@
 """General utility objects"""
 
-from stdatamodels.jwst.datamodels import dqflags
-import numpy as np
-
-
-def set_nans_to_donotuse(data, dq):
-    """Set all NaN values in the data that have an even value to
-    DO_NOT_USE.
-
-    Parameters
-    ----------
-    data : numpy array
-        The science data array to find NaN values and
-        check of these have a DQ flag=DO_NOT_USE, or
-        set it if not.
-
-    dq : numpy array
-        The DQ array to be checked.
-
-    Returns
-    -------
-    dq : numpy array
-        The updated DQ array.
-    """
-    dq[np.isnan(data)] |= dqflags.pixel['DO_NOT_USE']
-    return dq
-
 
 class LoggingContext:
     """Logging context manager
diff --git a/jwst/lib/pipe_utils.py b/jwst/lib/pipe_utils.py
index 168fa29d5..940b3bb1b 100644
--- a/jwst/lib/pipe_utils.py
+++ b/jwst/lib/pipe_utils.py
@@ -1,8 +1,16 @@
 """Pipeline utilities objects"""
 
+import logging
+
 import numpy as np
+from stdatamodels.properties import ObjectNode
+from stdatamodels.jwst.datamodels import dqflags, JwstDataModel
+
+from jwst.associations.lib.dms_base import TSO_EXP_TYPES
+
 
-from ..associations.lib.dms_base import TSO_EXP_TYPES
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
 
 
 def is_tso(model):
@@ -76,3 +84,70 @@ def is_irs2(model):
         return True
     else:
         return False
+
+
+def match_nans_and_flags(input_model):
+    """Ensure data, error, variance, and DQ are marked consistently for invalid data.
+
+    Invalid data is assumed to be any pixel set to NaN in any one of the
+    data, error, or variance arrays, or else set to the DO_NOT_USE flag
+    in the DQ array.
+
+    The input model is updated in place with NaNs or DO_NOT_USE flags, as
+    appropriate, at all invalid data locations.
+
+    Parameters
+    ----------
+    input_model : DataModel
+        Input model containing some combination of data, dq, err, var_rnoise,
+        var_poisson, and var_flat extensions. These extensions must all have
+        matching dimensions if present.
+    """
+    # Check for datamodel input or slit instance
+    if (not isinstance(input_model, JwstDataModel)
+            and not isinstance(input_model, ObjectNode)):
+        raise ValueError(f"Input {type(input_model)} is not a datamodel.")
+
+    # Build up the invalid data flags from each available data extension.
+    is_invalid = None
+    data_shape = None
+    nan_extensions = ['data', 'err', 'var_rnoise', 'var_poisson', 'var_flat']
+    for extension in nan_extensions:
+        if not hasattr(input_model, extension):
+            continue
+        data = getattr(input_model, extension)
+        if is_invalid is None:
+            is_invalid = np.isnan(data)
+            data_shape = data.shape
+        else:
+            if data.shape != data_shape:
+                log.warning(f"Mismatched data shapes; skipping invalid data "
+                            f"updates for extension '{extension}'")
+                continue
+            is_invalid |= np.isnan(data)
+
+    # Nothing to do if no extensions were found to update
+    if is_invalid is None:
+        return
+
+    # Add in invalid flags from the DQ extension if present
+    if hasattr(input_model, 'dq'):
+        do_not_use = (input_model.dq & dqflags.pixel['DO_NOT_USE']).astype(bool)
+        if input_model.dq.shape != data_shape:
+            log.warning("Mismatched data shapes; skipping invalid data "
+                        "updates for extension 'dq'")
+        else:
+            is_invalid |= do_not_use
+
+    # Update all the data extensions
+    for extension in nan_extensions:
+        if not hasattr(input_model, extension):
+            continue
+        data = getattr(input_model, extension)
+        if data.shape != data_shape:
+            continue
+        data[is_invalid] = np.nan
+
+    # Update the DQ extension
+    if input_model.dq.shape == data_shape:
+        input_model.dq[is_invalid] |= dqflags.pixel['DO_NOT_USE']
diff --git a/jwst/outlier_detection/ifu.py b/jwst/outlier_detection/ifu.py
index dfd7c7f6e..6ab12a015 100644
--- a/jwst/outlier_detection/ifu.py
+++ b/jwst/outlier_detection/ifu.py
@@ -23,7 +23,7 @@ This routine performs the following operations::
   3. For each input file store the  minimum of the pixel neighbor differences
   4. Comparing all the differences from all the input data find the minimum neighbor difference
   5. Normalize minimum difference to local median of difference array
-  6. select outliers by flagging those normailzed minimum values > threshold_percent
+  6. Select outliers by flagging those normalized minimum values > threshold_percent
   7. Updates input ImageModel DQ arrays with mask of detected outliers.
 """
 
@@ -32,6 +32,7 @@ import logging
 import numpy as np
 
 from jwst.datamodels import ModelContainer
+from jwst.lib.pipe_utils import match_nans_and_flags
 from jwst.stpipe.utilities import record_step_status
 from stdatamodels.jwst import datamodels
 from stdatamodels.jwst.datamodels import dqflags
@@ -105,6 +106,7 @@ def detect_outliers(
                       make_output_path)
     return input_models
 
+
 def flag_outliers(input_models, idet, uq_det, ndet_files,
                   diffaxis, nx, ny,
                   kern_size, threshold_percent,
@@ -258,6 +260,11 @@ def flag_outliers(input_models, idet, uq_det, ndet_files,
             total_bad = num_above + nadditional
             percent_cr = total_bad / (model.data.shape[0] * model.data.shape[1]) * 100
             log.info(f"Total #  pixels flagged as outliers: {total_bad} ({percent_cr:.2f}%)")
+
+            # Make sure all error and variance arrays also have matching
+            # NaNs and DQ flags
+            match_nans_and_flags(model)
+
             # update model
             input_models[i] = model
 
@@ -300,4 +307,3 @@ def create_optional_results_model(opt_info):
     opt_model.meta.kernel_ysize = kernsize_y
     opt_model.meta.threshold_percent = threshold_percent
     return opt_model
-
diff --git a/jwst/outlier_detection/utils.py b/jwst/outlier_detection/utils.py
index b14ab125c..649141571 100644
--- a/jwst/outlier_detection/utils.py
+++ b/jwst/outlier_detection/utils.py
@@ -5,6 +5,7 @@ import warnings
 
 import numpy as np
 
+from jwst.lib.pipe_utils import match_nans_and_flags
 from jwst.resample.resample import compute_image_pixel_area
 from stcal.outlier_detection.utils import compute_weight_threshold, gwcs_blot, flag_crs, flag_resampled_crs
 from stdatamodels.jwst import datamodels
@@ -245,6 +246,10 @@ def _flag_resampled_model_crs(
     input_model.dq |= cr_mask * np.uint32(DO_NOT_USE | OUTLIER)
     log.info(f"{np.count_nonzero(cr_mask)} pixels marked as outliers")
 
+    # Make sure all data, error, and variance arrays have
+    # matching NaNs and DQ flags
+    match_nans_and_flags(input_model)
+
 
 def flag_crs_in_models_with_resampling(
     input_models,
@@ -262,6 +267,12 @@ def flag_crs_in_models_with_resampling(
 
 def flag_model_crs(image, blot, snr):
     cr_mask = flag_crs(image.data, image.err, blot, snr)
-    # update dq array in-place
+
+    # Update dq array in-place
     image.dq |= cr_mask * np.uint32(DO_NOT_USE | OUTLIER)
+
+    # Make sure all data, error, and variance arrays have
+    # matching NaNs and DQ flags
+    match_nans_and_flags(image)
+
     log.info(f"{np.count_nonzero(cr_mask)} pixels marked as outliers")
diff --git a/jwst/pathloss/pathloss.py b/jwst/pathloss/pathloss.py
index 9586bc320..c3aba3b32 100644
--- a/jwst/pathloss/pathloss.py
+++ b/jwst/pathloss/pathloss.py
@@ -10,7 +10,7 @@ import stdatamodels.jwst.datamodels as datamodels
 
 from jwst.assign_wcs import nirspec, util
 from jwst.lib.wcs_utils import get_wavelengths
-from jwst.lib.basic_utils import set_nans_to_donotuse
+from jwst.lib.pipe_utils import match_nans_and_flags
 
 
 log = logging.getLogger(__name__)
@@ -547,8 +547,8 @@ def do_correction_mos(data, pathloss, inverse=False, source_type=None, correctio
         slit.pathloss_point = correction.pathloss_point
         slit.pathloss_uniform = correction.pathloss_uniform
 
-        # check the dq flags have the correct value
-        slit.dq = set_nans_to_donotuse(slit.data, slit.dq)
+        # Make sure all NaNs and flags match up in the output slit model
+        match_nans_and_flags(slit)
 
     # Set step status to complete
     data.meta.cal_step.pathloss = 'COMPLETE'
@@ -613,8 +613,8 @@ def do_correction_fixedslit(data, pathloss, inverse=False, source_type=None, cor
         slit.pathloss_point = correction.pathloss_point
         slit.pathloss_uniform = correction.pathloss_uniform
 
-        # check the dq flags have the correct value
-        slit.dq = set_nans_to_donotuse(slit.data, slit.dq)
+        # Make sure all NaNs and flags match up in the output slit model
+        match_nans_and_flags(slit)
 
     # Set step status to complete
     data.meta.cal_step.pathloss = 'COMPLETE'
@@ -669,8 +669,8 @@ def do_correction_ifu(data, pathloss, inverse=False, source_type=None, correctio
     # This might be useful to other steps
     data.wavelength = correction.wavelength
 
-    # check the dq flags have the correct value
-    data.dq = set_nans_to_donotuse(data.data, data.dq)
+    # Make sure all NaNs and flags match up in the output data model
+    match_nans_and_flags(data)
 
     # Set the step status to complete
     data.meta.cal_step.pathloss = 'COMPLETE'
@@ -713,8 +713,8 @@ def do_correction_lrs(data, pathloss, user_slit_loc):
     # This might be useful to other steps
     data.wavelength = correction.wavelength
 
-    # check the dq flags have the correct value
-    data.dq = set_nans_to_donotuse(data.data, data.dq)
+    # Make sure all NaNs and flags match up in the output data model
+    match_nans_and_flags(data)
 
     # Set the step status to complete
     data.meta.cal_step.pathloss = 'COMPLETE'
@@ -808,8 +808,8 @@ def do_correction_soss(data, pathloss):
         data.var_flat /= pathloss_2d**2
     data.pathloss_point = pathloss_2d
 
-    # check the dq flags have the correct value
-    data.dq = set_nans_to_donotuse(data.data, data.dq)
+    # Make sure all NaNs and flags match up in the output data model
+    match_nans_and_flags(data)
 
     # Set step status to complete
     data.meta.cal_step.pathloss = 'COMPLETE'
diff --git a/jwst/photom/photom.py b/jwst/photom/photom.py
index 7628582e3..f8ec16567 100644
--- a/jwst/photom/photom.py
+++ b/jwst/photom/photom.py
@@ -7,6 +7,7 @@ from astropy import units as u
 from stdatamodels.jwst import datamodels
 from stdatamodels.jwst.datamodels import dqflags
 
+from .. lib.pipe_utils import match_nans_and_flags
 from .. lib.wcs_utils import get_wavelengths
 from .. lib.dispaxis import get_dispersion_direction
 from . import miri_mrs
@@ -940,6 +941,9 @@ class DataSet():
                 self.input.meta.bunit_data = 'DN/s'
                 self.input.meta.bunit_err = 'DN/s'
 
+            # Make sure output model has consistent NaN and DO_NOT_USE values
+            match_nans_and_flags(slit)
+
         elif isinstance(self.input, datamodels.MultiSpecModel):
             # Does this block need to address SB columns as well, or will
             # they (presumably) never be populated for SOSS?
@@ -993,6 +997,9 @@ class DataSet():
                 self.input.meta.bunit_data = 'DN/s'
                 self.input.meta.bunit_err = 'DN/s'
 
+            # Make sure output model has consistent NaN and DO_NOT_USE values
+            match_nans_and_flags(self.input)
+
         return
 
     def create_2d_conversion(self, model, exptype, conversion, waves, relresps,
diff --git a/jwst/resample/resample_spec_step.py b/jwst/resample/resample_spec_step.py
index bcc3d1a95..7aa406091 100755
--- a/jwst/resample/resample_spec_step.py
+++ b/jwst/resample/resample_spec_step.py
@@ -4,11 +4,13 @@ from stdatamodels.jwst import datamodels
 from stdatamodels.jwst.datamodels import MultiSlitModel, ImageModel
 
 from jwst.datamodels import ModelContainer, ModelLibrary
+from jwst.lib.pipe_utils import match_nans_and_flags
+from jwst.lib.wcs_utils import get_wavelengths
+
 from . import resample_spec, ResampleStep
 from ..exp_to_source import multislit_to_container
 from ..assign_wcs.util import update_s_region_spectral
 from ..stpipe import Step
-from jwst.lib.wcs_utils import get_wavelengths
 
 
 # Force use of all DQ flagged data except for DO_NOT_USE and NON_SCIENCE
@@ -132,6 +134,10 @@ class ResampleSpecStep(Step):
 
         pscale_ratio = None
         for container in containers.values():
+            # Make sure all input models have consistent NaN and DO_NOT_USE values
+            for model in container:
+                match_nans_and_flags(model)
+
             resamp = resample_spec.ResampleSpecData(container, **self.drizpars)
 
             library = ModelLibrary(container, on_disk=False)
@@ -204,6 +210,9 @@ class ResampleSpecStep(Step):
         result : `~jwst.datamodels.SlitModel`
             The resampled output
         """
+        # Make sure all input models have consistent NaN and DO_NOT_USE values
+        for model in input_models:
+            match_nans_and_flags(model)
 
         resamp = resample_spec.ResampleSpecData(input_models, **self.drizpars)
 
@@ -243,4 +252,4 @@ class ResampleSpecStep(Step):
                 pass
             else:
                 if val is not None:
-                    setattr(model, attr, val)
\ No newline at end of file
+                    setattr(model, attr, val)
diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py
index a2fd35db9..5790fb4f5 100755
--- a/jwst/resample/resample_step.py
+++ b/jwst/resample/resample_step.py
@@ -5,6 +5,7 @@ from copy import deepcopy
 import asdf
 
 from jwst.datamodels import ModelLibrary, ImageModel
+from jwst.lib.pipe_utils import match_nans_and_flags
 
 from . import resample
 from ..stpipe import Step
@@ -88,6 +89,12 @@ class ResampleStep(Step):
                 raise RuntimeError(f"Input {example_model} is not a 2D image.")
             del example_model
 
+            # Make sure all input models have consistent NaN and DO_NOT_USE values
+            for model in input_models:
+                match_nans_and_flags(model)
+                input_models.shelve(model)
+            del model
+
         # Setup drizzle-related parameters
         kwargs = self.get_drizpars()
 


Test patch:
diff --git a/jwst/flatfield/tests/test_flatfield.py b/jwst/flatfield/tests/test_flatfield.py
index a332632bf..71a5f8b43 100644
--- a/jwst/flatfield/tests/test_flatfield.py
+++ b/jwst/flatfield/tests/test_flatfield.py
@@ -60,7 +60,18 @@ def test_flatfield_step_interface(instrument, exptype):
     finally:
         FlatFieldStep.reference_file_types = previous_reference_file_types
 
-    assert (result.data == data.data).all()
+    # flat is unity, so data should match up except where invalid
+    is_nan = np.isnan(result.data)
+    assert np.all(result.data[~is_nan] == data.data[~is_nan])
+
+    # where data is invalid, error and variance are also invalid
+    # and DQ is set to DO_NOT_USE
+    assert np.all(np.isnan(result.err[is_nan]))
+    assert np.all(np.isnan(result.var_rnoise[is_nan]))
+    assert np.all(np.isnan(result.var_poisson[is_nan]))
+    assert np.all(np.isnan(result.var_flat[is_nan]))
+    assert np.all(result.dq[is_nan] | datamodels.dqflags.pixel['DO_NOT_USE'])
+
     assert result.var_flat.shape == shape
     assert result.meta.cal_step.flat_field == 'COMPLETE'
 
diff --git a/jwst/lib/tests/test_pipe_utils.py b/jwst/lib/tests/test_pipe_utils.py
index 00a87688b..cdd840ca3 100644
--- a/jwst/lib/tests/test_pipe_utils.py
+++ b/jwst/lib/tests/test_pipe_utils.py
@@ -4,6 +4,7 @@ import pytest
 import numpy as np
 
 from stdatamodels.jwst import datamodels
+from stdatamodels.properties import ObjectNode
 
 from jwst.lib import pipe_utils
 from jwst.associations.lib import dms_base
@@ -84,3 +85,134 @@ def test_is_irs2_2():
     x = np.ones((3200, 2), dtype=np.float32)
     model = datamodels.ImageModel(data=x)
     assert pipe_utils.is_irs2(model)
+
+
+@pytest.mark.parametrize('shape', [(10,), (10, 10), (10, 10, 10)])
+@pytest.mark.parametrize('extensions',
+                         [('data',), ('data', 'dq'),
+                          ('data', 'err', 'var_rnoise'),
+                          ('data', 'err', 'var_rnoise', 'var_poisson',
+                           'var_flat', 'dq')])
+def test_match_nans_and_flags(shape, extensions):
+    # Set up a model with matching data shapes, variable extensions
+    dnu = datamodels.dqflags.pixel['DO_NOT_USE']
+    model = datamodels.SlitModel()
+    invalid = np.full(shape, False)
+    for i, extname in enumerate(extensions):
+        if extname == 'dq':
+            data = np.zeros(shape, dtype=np.uint32)
+            data.flat[i] = dnu
+        else:
+            data = np.ones(shape)
+            data.flat[i] = np.nan
+        setattr(model, extname, data)
+        invalid.flat[i] = True
+
+    # Match flags and NaNs across all extensions
+    pipe_utils.match_nans_and_flags(model)
+
+    # Check that all extensions have all invalid flags
+    for extname in extensions:
+        data = getattr(model, extname)
+        if extname == 'dq':
+            assert np.all(data[invalid] == dnu)
+            assert np.all(data[~invalid] == 0)
+        else:
+            assert np.all(np.isnan(data[invalid]))
+            assert np.all(data[~invalid] == 1)
+
+    model.close()
+
+
+def test_match_nans_and_flags_dq_only():
+    # Set up a model with only a DQ array, to test edge conditions
+    dnu = datamodels.dqflags.pixel['DO_NOT_USE']
+    model = datamodels.MaskModel()
+    shape = (10, 10)
+    invalid = np.full(shape, False)
+
+    model.dq = np.zeros(shape, dtype=np.uint32)
+    model.dq.flat[5] = dnu
+    invalid.flat[5] = True
+
+    # Match flags and NaNs across all extensions
+    dq_copy = model.dq.copy()
+    pipe_utils.match_nans_and_flags(model)
+
+    # Check that DQ extension is unchanged
+    assert np.all(model.dq == dq_copy)
+
+    model.close()
+
+
+def test_match_nans_and_flags_no_data(caplog):
+    # Set up a model with no data, to test edge conditions
+    model = datamodels.JwstDataModel()
+    model_copy = model.copy()
+    pipe_utils.match_nans_and_flags(model)
+
+    # nothing happens
+    assert model == model_copy
+    assert 'WARNING' not in caplog.text
+    assert 'ERROR' not in caplog.text
+
+    model.close()
+    model_copy.close()
+
+
+def test_match_nans_and_flags_not_a_model():
+    # Make sure the function throws an error if input is not a model
+    model = "not_a_model"
+    with pytest.raises(ValueError, match="not a datamodel"):
+        pipe_utils.match_nans_and_flags(model)
+
+    # empty JwstDataModel is okay
+    model = datamodels.JwstDataModel()
+    pipe_utils.match_nans_and_flags(model)
+
+    # a datamodel object node is also okay
+    multislit = datamodels.MultiSlitModel()
+    multislit.slits.append(datamodels.SlitModel())
+    assert not isinstance(multislit.slits[0], datamodels.JwstDataModel)
+    assert isinstance(multislit.slits[0], ObjectNode)
+    pipe_utils.match_nans_and_flags(multislit.slits[0])
+
+
+def test_match_nans_and_flags_shape_mismatch():
+    # Set up a model with mismatched data shapes
+    dnu = datamodels.dqflags.pixel['DO_NOT_USE']
+    model = datamodels.SlitModel()
+
+    shape = (10, 10)
+    bad_shape = (10, 5)
+    invalid = np.full(shape, False)
+    extensions = ['data', 'err', 'dq', 'var_poisson']
+    for i, extname in enumerate(extensions):
+        if extname == 'dq':
+            data = np.zeros(bad_shape, dtype=np.uint32)
+            data.flat[i] = dnu
+        elif extname == 'var_poisson':
+            data = np.ones(bad_shape, dtype=np.uint32)
+            data.flat[i] = dnu
+        else:
+            data = np.ones(shape)
+            data.flat[i] = np.nan
+            invalid.flat[i] = True
+        setattr(model, extname, data)
+
+    # Match flags and NaNs across all extensions:
+    # will warn for data mismatch
+    model_copy = model.copy()
+    pipe_utils.match_nans_and_flags(model)
+
+    # Only data and err are updated
+    for extname in extensions:
+        data = getattr(model, extname)
+        if extname == 'data' or extname == 'err':
+            assert np.all(np.isnan(data[invalid]))
+            assert np.all(data[~invalid] == 1)
+        else:
+            assert np.allclose(data, getattr(model_copy, extname))
+
+    model.close()
+    model_copy.close()
diff --git a/jwst/outlier_detection/tests/test_outlier_detection.py b/jwst/outlier_detection/tests/test_outlier_detection.py
index 7cb74f5a8..b3b8c9b80 100644
--- a/jwst/outlier_detection/tests/test_outlier_detection.py
+++ b/jwst/outlier_detection/tests/test_outlier_detection.py
@@ -88,7 +88,10 @@ def test_flag_cr(sci_blot_image_pair):
     )
 
     # Make sure science data array is unchanged after flag_cr()
-    np.testing.assert_allclose(sci.data, data_copy)
+    # except outliers are NaN
+    dnu = (sci.dq & OUTLIER_DO_NOT_USE).astype(bool)
+    assert np.all(np.isnan(sci.data[dnu]))
+    assert np.allclose(sci.data[~dnu], data_copy[~dnu])
 
     # Verify that both DQ flags are set in the DQ array for all outliers
     assert sci.dq[3, 3] == OUTLIER_DO_NOT_USE
@@ -224,21 +227,21 @@ def test_outlier_step_base(we_three_sci, tmp_cwd):
     assert len(median_files) == 0
 
     # Save all the data into a separate array before passing into step
-    data_as_cube = container.map_function(lambda model, index: model.data.copy(), modify=False)
+    data_as_cube = list(container.map_function(
+        lambda model, index: model.data.copy(), modify=False))
 
     result = OutlierDetectionStep.call(
         container, save_results=True, save_intermediate_results=True
     )
 
-    # Make sure nothing changed in SCI array
-    result.map_function(
-        lambda model, index: np.testing.assert_allclose(data_as_cube[index], model.data),
-        modify=False,
-    )
-
-    # Verify source is not flagged
     with result:
-        for r in result:
+        for i, r in enumerate(result):
+            # Make sure nothing changed in SCI array except outliers are NaN
+            dnu = (r.dq & OUTLIER_DO_NOT_USE).astype(bool)
+            assert np.all(np.isnan(r.data[dnu]))
+            assert np.allclose(data_as_cube[i][~dnu], r.data[~dnu])
+
+            # Verify source is not flagged
             assert r.dq[7, 7] == datamodels.dqflags.pixel["GOOD"]
             result.shelve(r, modify=False)
 
@@ -378,21 +381,21 @@ def test_outlier_step_on_disk(three_sci_as_asn, tmp_cwd):
     container = ModelLibrary(three_sci_as_asn, on_disk=True)
 
     # Save all the data into a separate array before passing into step
-    data_as_cube = container.map_function(lambda model, index: model.data.copy(), modify=False)
+    data_as_cube = list(container.map_function(
+        lambda model, index: model.data.copy(), modify=False))
 
     result = OutlierDetectionStep.call(
         container, save_results=True, save_intermediate_results=True, in_memory=False
     )
 
-    # Make sure nothing changed in SCI array
-    result.map_function(
-        lambda model, index: np.testing.assert_allclose(data_as_cube[index], model.data),
-        modify=False,
-    )
-
-    # Verify source is not flagged
     with result:
-        for r in result:
+        for i, r in enumerate(result):
+            # Make sure nothing changed in SCI array except outliers are NaN
+            dnu = (r.dq & OUTLIER_DO_NOT_USE).astype(bool)
+            assert np.all(np.isnan(r.data[dnu]))
+            assert np.allclose(data_as_cube[i][~dnu], r.data[~dnu])
+
+            # Verify source is not flagged
             assert r.dq[7, 7] == datamodels.dqflags.pixel["GOOD"]
             result.shelve(r, modify=False)
 
@@ -440,7 +443,7 @@ def test_outlier_step_square_source_no_outliers(we_three_sci, tmp_cwd):
 
 
 @pytest.mark.parametrize("exptype", IMAGE_MODES)
-def test_outlier_step_image_weak_CR_dither(exptype, tmp_cwd):
+def test_outlier_step_image_weak_cr_dither(exptype, tmp_cwd):
     """Test whole step with an outlier for imaging modes"""
     bkg = 1.5
     sig = 0.02
@@ -464,15 +467,14 @@ def test_outlier_step_image_weak_CR_dither(exptype, tmp_cwd):
 
     result = OutlierDetectionStep.call(container, in_memory=True)
 
-    # Make sure nothing changed in SCI array
-    result.map_function(
-        lambda model, index: np.testing.assert_allclose(data_as_cube[index], model.data),
-        modify=False,
-    )
-
-    # Verify source is not flagged
     with result:
-        for r in result:
+        for i, r in enumerate(result):
+            # Make sure nothing changed in SCI array except outliers are NaN
+            dnu = (r.dq & OUTLIER_DO_NOT_USE).astype(bool)
+            assert np.all(np.isnan(r.data[dnu]))
+            assert np.allclose(data_as_cube[i][~dnu], r.data[~dnu])
+
+            # Verify source is not flagged
             assert r.dq[7, 7] == datamodels.dqflags.pixel["GOOD"]
             result.shelve(r, modify=False)
 
@@ -484,7 +486,7 @@ def test_outlier_step_image_weak_CR_dither(exptype, tmp_cwd):
 
 
 @pytest.mark.parametrize("exptype, tsovisit", exptypes_coron)
-def test_outlier_step_image_weak_CR_coron(exptype, tsovisit, tmp_cwd):
+def test_outlier_step_image_weak_cr_coron(exptype, tsovisit, tmp_cwd):
     """Test whole step with an outlier for coronagraphic modes"""
     bkg = 1.5
     sig = 0.02
@@ -503,15 +505,19 @@ def test_outlier_step_image_weak_CR_coron(exptype, tsovisit, tmp_cwd):
 
     result = OutlierDetectionStep.call(cube)
 
-    # Make sure nothing changed in SCI array
+    # Make sure nothing changed in SCI array except that
+    # outliers are NaN
     for i, image in enumerate(container):
-        np.testing.assert_allclose(image.data, result.data[i])
+        dnu = (result.dq[i] & OUTLIER_DO_NOT_USE).astype(bool)
+        assert np.all(np.isnan(result.data[i][dnu]))
+        assert np.allclose(image.data[~dnu], result.data[i][~dnu])
 
     # Verify source is not flagged
     assert np.all(result.dq[:, 7, 7] == datamodels.dqflags.pixel["GOOD"])
 
     # Verify CR is flagged
-    assert np.all(result.dq[0, 12, 12] == OUTLIER_DO_NOT_USE)
+    assert result.dq[0, 12, 12] == OUTLIER_DO_NOT_USE
+    assert np.isnan(result.data[0, 12, 12])
 
 
 @pytest.mark.parametrize("exptype, tsovisit", exptypes_tso)
@@ -543,9 +549,12 @@ def test_outlier_step_weak_cr_tso(exptype, tsovisit):
 
     result = OutlierDetectionStep.call(cube, rolling_window_width=rolling_window_width)
 
-    # Make sure nothing changed in SCI array
+    # Make sure nothing changed in SCI array except
+    # that outliers are NaN
     for i, model in enumerate(im):
-        np.testing.assert_allclose(model.data, result.data[i])
+        dnu = (result.dq[i] & OUTLIER_DO_NOT_USE).astype(bool)
+        assert np.all(np.isnan(result.data[i][dnu]))
+        assert np.allclose(model.data[~dnu], result.data[i][~dnu])
 
     # Verify source is not flagged
     assert np.all(result.dq[:, 7, 7] == datamodels.dqflags.pixel["GOOD"])
diff --git a/jwst/regtest/test_nirspec_image2.py b/jwst/regtest/test_nirspec_image2.py
index 5bb46cd75..fba5f82ea 100644
--- a/jwst/regtest/test_nirspec_image2.py
+++ b/jwst/regtest/test_nirspec_image2.py
@@ -45,7 +45,12 @@ def test_ff_inv(rtdata, fitsdiff_default_kwargs):
         flatted = FlatFieldStep.call(data)
         unflatted = FlatFieldStep.call(flatted, inverse=True)
 
-        assert np.allclose(data.data, unflatted.data), 'Inversion failed'
+        # flat fielding may set some new NaN values - ignore these in test
+        is_nan = np.isnan(unflatted.data)
+        assert np.allclose(data.data[~is_nan], unflatted.data[~is_nan]), 'Inversion failed'
+
+        # make sure NaNs are only at do_not_use pixels
+        assert np.all(unflatted.dq[is_nan] & dm.dqflags.pixel['DO_NOT_USE'])
 
 
 @pytest.mark.bigdata
@@ -62,4 +67,9 @@ def test_correction_pars(rtdata, fitsdiff_default_kwargs):
         step.use_correction_pars = True
         reflatted = step.run(data)
 
-    assert np.allclose(flatted.data, reflatted.data), 'Re-run with correction parameters failed'
+        # flat fielding may set some new NaN values - ignore these in test
+        is_nan = np.isnan(reflatted.data)
+        assert np.allclose(flatted.data[~is_nan], reflatted.data[~is_nan]), 'Re-run with correction parameters failed'
+
+        # make sure NaNs are only at do_not_use pixels
+        assert np.all(reflatted.dq[is_nan] & dm.dqflags.pixel['DO_NOT_USE'])