Repository: spacetelescope/jwst

Issue:
Propagate 1D flat errors
_Issue [JP-3429](https://jira.stsci.edu/browse/JP-3429) was created on JIRA by [Melanie Clarke](https://jira.stsci.edu/secure/ViewProfile.jspa?name=mclarke):_

NIRSpec is updating its flat files to include a 1D error array, associated with the fast variation component in the F-flats (JP-3304).  The flat field step needs to be updated to propagate this error from the reference files to the output error arrays.

In particular, the combine_fast_slow function needs to be updated to take an error array as input and interpolate it onto the wavelengths for the observation.  The interpolated error should then be added in quadrature to any error from other components.

Two implementation notes:

1. Both zero and NaN are possible values for all flat error arrays, 2D and 1D.  The S and D flats will not have error arrays for their 1D components, so these values will default to NaN; their 2D error components will be set to zero, since all flat errors are accounted for in the F-flat instead. When combining flat error values, NaN values should be ignored and zero values should be treated as valid numbers.  

2. The 1D flat vector in the combine_fast_slow algorithm is interpolated with a 3-point Gaussian integration.  A rigorous error propagation through that step should take into account cross-correlations between the pixels that are weighted and summed together, but we don't recommend that.  The flat errors are not likely to be known to enough precision to justify the complexity.  Instead, we recommend interpolating the error vector onto the wavelengths with a simple linear interpolation, e.g.:
  error_value = np.interp(wavelengths, tab_wl, tab_flat_error, left=np.nan, right=np.nan)

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/7978?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
Attention: `17 lines` in your changes are missing coverage. Please review.

| [Files](https://app.codecov.io/gh/spacetelescope/jwst/pull/7978?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | |
|---|---|---|
| [jwst/flatfield/flat\_field.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/7978?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9mbGF0ZmllbGQvZmxhdF9maWVsZC5weQ==) | `78.98% <39.28%> (-1.82%)` | :arrow_down: |


:loudspeaker: Thoughts on this report? [Let us know!](https://about.codecov.io/pull-request-comment-report/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

braingram: Thanks! I started regression tests for this PR at:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/975/
hbushouse: Let me know when you're ready for another regression test run and I'll fire one off.
melanieclarke: Thanks @hbushouse - I'm done with my local testing, so I think this is ready for another regression test run.
melanieclarke: > And I assume this is all backwards compatible with old F-flats that don't have the error column in the table?

Yes, this should work for the old F-flats.  The error values are all assigned to NaN when they're read in, which is handled as zero error in the code.  I've tested with the existing CRDS flats and see no change in the output, compared to products reduced with v1.12.1.

hbushouse: I'll start another regtest run as soon as the queue opens up.
hbushouse: Regression tests started at https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/980
hbushouse: Regression test results are at https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/980/testReport/  @melanieclarke can you see those results (not sure if you have permission)?
They show differences in the VAR_FLAT and ERR array values for many NIRSpec fixed-slit and IFU test results, but none for any MOS mode tests (not sure why). The differences are small numerical differences, or differences showing a value where there used to be a zero or NaN, all of which seems expected due to the changes here in the way the various error estimates are propagated (even without using new ref files).
melanieclarke: @hbushouse - I have permission enough to see the test results.  I'm pretty sure those changes are as expected, but I'll take a closer look to confirm.
melanieclarke: For the values where there used to be NaN: this is due to the change in NaN handling for input flat values.  All NaNs are now set to zero so that NaNs in flat variance don't contaminate later error estimates.

For the numerical differences: from the filenames, it looks like the inputs are all pre-flight data, for which there are still some associated flats that do not have all 2D error values set to 0.  The difference comes from a change in assumption about the meaning of the 2D error. Before, the code assumed that the 2D error was the only source of error for each flat type, and it was normalized by the composed slow * fast components before being added in quadrature to the other components.  Now, it is assumed that 2D errors apply only to the 2D component, since 1D errors are provided separately. This changes the normalization for the non-zero errors.

This change is moot for all flight data: 2D errors are currently set to zero in CRDS flats, and we do not expect to provide non-zero 2D errors going forward.  However, this does resolve a long-standing difference in the NIRSpec validation tests, which also use pre-flight data (JP-3103).  Error calculations from this change now match the validation code.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index b416fb8c6..95802375e 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -17,6 +17,12 @@ set_telescope_pointing
 
 - Ensure that CRPIX1/2 are set to default values for Guiding modes that fail pointing determination [#7983]
 
+flat_field
+----------
+
+- Update the ``combine_fast_slow`` function for NIRSpec spectroscopic flats 
+   to use 1D error values provided by F-flat reference files. [#7978]
+
 1.12.1 (2023-09-26)
 ===================
 
diff --git a/jwst/flatfield/flat_field.py b/jwst/flatfield/flat_field.py
index c222b342e..c7b393a97 100644
--- a/jwst/flatfield/flat_field.py
+++ b/jwst/flatfield/flat_field.py
@@ -434,7 +434,7 @@ def nirspec_fs_msa(output_model, f_flat_model, s_flat_model, d_flat_model, dispa
         flat_data_squared = slit_flat.data ** 2
         slit.var_poisson /= flat_data_squared
         slit.var_rnoise /= flat_data_squared
-        slit.var_flat = slit.data ** 2 / flat_data_squared * slit_flat.err ** 2
+        slit.var_flat = (slit.data / slit_flat.data * slit_flat.err) ** 2
         slit.err = np.sqrt(slit.var_poisson + slit.var_rnoise + slit.var_flat)
 
         # Combine the science and flat DQ arrays
@@ -510,7 +510,7 @@ def nirspec_brightobj(output_model, f_flat_model, s_flat_model, d_flat_model, di
     flat_data_squared = interpolated_flat.data ** 2
     output_model.var_poisson /= flat_data_squared
     output_model.var_rnoise /= flat_data_squared
-    output_model.var_flat = output_model.data ** 2 / flat_data_squared * interpolated_flat.err ** 2
+    output_model.var_flat = (output_model.data / interpolated_flat.data * interpolated_flat.err) ** 2
     output_model.err = np.sqrt(
         output_model.var_poisson + output_model.var_rnoise + output_model.var_flat
     )
@@ -576,7 +576,7 @@ def nirspec_ifu(output_model, f_flat_model, s_flat_model, d_flat_model, dispaxis
         flat_data_squared = flat ** 2
         output_model.var_poisson /= flat_data_squared
         output_model.var_rnoise /= flat_data_squared
-        output_model.var_flat = output_model.data ** 2 / flat_data_squared * flat_err ** 2
+        output_model.var_flat = (output_model.data / flat * flat_err) ** 2
         output_model.err = np.sqrt(
             output_model.var_poisson + output_model.var_rnoise + output_model.var_flat
         )
@@ -656,6 +656,7 @@ def create_flat_field(wl, f_flat_model, s_flat_model, d_flat_model,
     s_flat, s_flat_dq, s_flat_err = spectrograph_flat(
         wl, s_flat_model, xstart, xstop, ystart, ystop,
         exposure_type, dispaxis, slit_name)
+
     d_flat, d_flat_dq, d_flat_err = detector_flat(
         wl, d_flat_model, xstart, xstop, ystart, ystop,
         exposure_type, dispaxis, slit_name)
@@ -665,16 +666,17 @@ def create_flat_field(wl, f_flat_model, s_flat_model, d_flat_model,
                          default_shape=flat_2d.shape)
 
     # Combine the uncertainty arrays, excluding the ones that are None
+    # Divide error by flat before squaring to avoid overflow
     sum_var = np.zeros_like(flat_2d)
     if f_flat_err is not None:
         np.place(f_flat, f_flat == 0, 1.)
-        sum_var += f_flat_err ** 2 / f_flat ** 2
+        sum_var += (f_flat_err / f_flat) ** 2
     if s_flat_err is not None:
         np.place(s_flat, s_flat == 0, 1.)
-        sum_var += s_flat_err ** 2 / s_flat ** 2
+        sum_var += (s_flat_err / s_flat) ** 2
     if d_flat_err is not None:
         np.place(d_flat, d_flat == 0, 1.)
-        sum_var += d_flat_err ** 2 / d_flat ** 2
+        sum_var += (d_flat_err / d_flat) ** 2
     flat_err = flat_2d * np.sqrt(sum_var)
 
     mask = np.bitwise_and(flat_dq, dqflags.pixel['DO_NOT_USE'])
@@ -731,8 +733,8 @@ def fore_optics_flat(wl, f_flat_model, exposure_type, dispaxis,
     else:
         quadrant = slit_nt.quadrant - 1  # convert to zero indexed
 
-    (tab_wl, tab_flat) = read_flat_table(f_flat_model, exposure_type,
-                                         slit_name, quadrant)
+    (tab_wl, tab_flat, tab_flat_err) = read_flat_table(f_flat_model, exposure_type,
+                                                       slit_name, quadrant)
     if tab_wl.max() < MICRONS_100:
         log.warning("Wavelengths in f_flat table appear to be in meters")
 
@@ -776,10 +778,10 @@ def fore_optics_flat(wl, f_flat_model, exposure_type, dispaxis,
         tab_flat *= np.interp(tab_wl, image_wl, one_d_flat, 1., 1.)
         f_flat_err = np.interp(wl, image_wl, one_d_err, 0., 0.)
 
+    # Combine 2D and 1D components with error propagation.
     # The shape of the output array is obtained from `wl`.
-
-    f_flat, f_flat_dq = combine_fast_slow(wl, flat_2d, f_flat_dq,
-                                          tab_wl, tab_flat, dispaxis)
+    f_flat, f_flat_dq, f_flat_err = combine_fast_slow(
+        wl, flat_2d, f_flat_dq, f_flat_err, tab_wl, tab_flat, tab_flat_err, dispaxis)
 
     # Find pixels in the flat that have a value of NaN and add to
     # DQ mask, DO_NOT_USE + NO_FLAT_FIELD
@@ -787,7 +789,7 @@ def fore_optics_flat(wl, f_flat_model, exposure_type, dispaxis,
     flat_nan = np.isnan(f_flat)
     f_flat_dq[flat_nan] = np.bitwise_or(f_flat_dq[flat_nan], bad_flag)
 
-    # Find pixels in the flat have have a value of zero, and add to
+    # Find pixels in the flat that have a value of zero, and add to
     # DQ mask,  DO_NOT_USE + NO_FLAT_FIELD
     flat_zero = np.where(f_flat == 0.)
     f_flat_dq[flat_zero] = np.bitwise_or(f_flat_dq[flat_zero], bad_flag)
@@ -856,7 +858,7 @@ def spectrograph_flat(wl, s_flat_model,
     if xstart >= xstop or ystart >= ystop:
         return 1., None
 
-    (tab_wl, tab_flat) = read_flat_table(s_flat_model, exposure_type,
+    (tab_wl, tab_flat, tab_flat_err) = read_flat_table(s_flat_model, exposure_type,
                                          slit_name, quadrant)
     if tab_wl.max() < MICRONS_100:
         log.warning("Wavelengths in s_flat table appear to be in meters")
@@ -902,8 +904,9 @@ def spectrograph_flat(wl, s_flat_model,
     # correction is made
     flat_2d[np.where(flat_bad)] = 1.0
 
-    s_flat, s_flat_dq = combine_fast_slow(wl, flat_2d, s_flat_dq,
-                                          tab_wl, tab_flat, dispaxis)
+    # Combine 2D and 1D components with error propagation
+    s_flat, s_flat_dq, s_flat_err = combine_fast_slow(
+        wl, flat_2d, s_flat_dq, s_flat_err, tab_wl, tab_flat, tab_flat_err, dispaxis)
 
     return s_flat, s_flat_dq, s_flat_err
 
@@ -963,7 +966,7 @@ def detector_flat(wl, d_flat_model,
     if xstart >= xstop or ystart >= ystop:
         return 1., None
 
-    (tab_wl, tab_flat) = read_flat_table(d_flat_model, exposure_type,
+    (tab_wl, tab_flat, tab_flat_err) = read_flat_table(d_flat_model, exposure_type,
                                          slit_name, quadrant)
     if tab_wl.max() < MICRONS_100:
         log.warning("Wavelengths in d_flat table appear to be in meters.")
@@ -1001,8 +1004,9 @@ def detector_flat(wl, d_flat_model,
     # correction is made
     flat_2d[np.where(flat_bad)] = 1.0
 
-    d_flat, d_flat_dq = combine_fast_slow(wl, flat_2d, d_flat_dq,
-                                          tab_wl, tab_flat, dispaxis)
+    # Combine 2D and 1D components with error propagation
+    d_flat, d_flat_dq, d_flat_err = combine_fast_slow(
+        wl, flat_2d, d_flat_dq, d_flat_err, tab_wl, tab_flat, tab_flat_err, dispaxis)
     return d_flat, d_flat_dq, d_flat_err
 
 
@@ -1146,6 +1150,10 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
     tab_flat : ndarray, 1-D, float
         The column of flat_field values read from the fast-variation table.
         `tab_wl` and `tab_flat` should be the same length.
+
+    tab_flat_err : ndarray, 1-D
+        The column of flat_field error values read from the fast-variation
+        table.
     """
 
     if quadrant is not None:  # NRS_MSASPEC
@@ -1163,6 +1171,7 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
         nelem_col = None
     wl_col = data["wavelength"]
     flat_col = data["data"]
+    flat_err_col = data["error"]
 
     nrows = len(wl_col)
     row = None
@@ -1188,6 +1197,7 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
         # Table contains arrays; use the row that was found above.
         tab_wl = wl_col[row].copy()
         tab_flat = flat_col[row].copy()
+        tab_flat_err = flat_err_col[row].copy()
         if nelem_col is not None:
             nelem = nelem_col[row]
     else:
@@ -1196,12 +1206,14 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
             # Table contains arrays, but there should be only one row.
             tab_wl = wl_col[0].copy()
             tab_flat = flat_col[0].copy()
+            tab_flat_err = flat_err_col[0].copy()
             if nelem_col is not None:
                 nelem = nelem_col[0]
         else:
             # Table contains scalar columns.
             tab_wl = wl_col.copy()
             tab_flat = flat_col.copy()
+            tab_flat_err = flat_err_col.copy()
             if nelem_col is not None:
                 nelem = nelem_col[0]  # arbitrary choice of row
     if nelem is not None:
@@ -1213,6 +1225,7 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
         else:
             tab_wl = tab_wl[:nelem]
             tab_flat = tab_flat[:nelem]
+            tab_flat_err = tab_flat_err[:nelem]
     else:
         nelem = len(tab_wl)
 
@@ -1227,6 +1240,7 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
                   "%d NaNs; these have been skipped.", nelem - n1)
         tab_wl = tab_wl[filter]
         tab_flat = tab_flat[filter]
+        tab_flat_err = tab_flat_err[filter]
     del filter1, filter2, filter
     # Skip zero or negative wavelengths, and skip zero flat-field values.
     filter1 = (tab_wl > 0.)
@@ -1239,6 +1253,7 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
                   n1 - n2)
         tab_wl = tab_wl[filter]
         tab_flat = tab_flat[filter]
+        tab_flat_err = tab_flat_err[filter]
     del filter1, filter2, filter
 
     # Check that the wavelengths are increasing.  This is a requirement
@@ -1249,10 +1264,10 @@ def read_flat_table(flat_model, exposure_type, slit_name=None, quadrant=None):
             log.warning("Wavelengths in the fast-variation table "
                         "must be strictly increasing.")
 
-    return tab_wl, tab_flat
+    return tab_wl, tab_flat, tab_flat_err
 
 
-def combine_fast_slow(wl, flat_2d, flat_dq, tab_wl, tab_flat, dispaxis):
+def combine_fast_slow(wl, flat_2d, flat_dq, flat_err, tab_wl, tab_flat, tab_flat_error, dispaxis):
     """Multiply the image by the tabular values.
 
     Parameters
@@ -1271,6 +1286,12 @@ def combine_fast_slow(wl, flat_2d, flat_dq, tab_wl, tab_flat, dispaxis):
         on the fast-variation component, and the updated array will be
         returned.
 
+    flat_err : ndarray or None
+        If not None, the error array corresponding to `flat_2d`.
+        A copy of this will be updated with errors added in quadrature from
+        the fast-variation component, and the updated array will be
+        returned.
+
     tab_wl : ndarray, 1-D
         Wavelengths corresponding to `tab_flat`.
 
@@ -1279,6 +1300,9 @@ def combine_fast_slow(wl, flat_2d, flat_dq, tab_wl, tab_flat, dispaxis):
         is the "fast" variation of the flat, i.e. fast with respect to
         wavelength.
 
+    tab_flat_err : ndarray, 1-D
+        The flat field error from the table part of the reference file.
+
     dispaxis : int
         1 is horizontal, 2 is vertical.
 
@@ -1294,6 +1318,10 @@ def combine_fast_slow(wl, flat_2d, flat_dq, tab_wl, tab_flat, dispaxis):
         within the range of `tab_wl`, NO_FLAT_FIELD will be used to flag
         this condition, and the fast-variation flat field value at that
         pixel will be set to 1.
+
+    flat_error : ndarray, 2-D, float32
+        The `tab_flat_err` values interpolated to the wavelengths of the
+        science image, i.e. `wl`. Missing values are set to 0.
     """
 
     wl_c = clean_wl(wl, dispaxis)
@@ -1304,6 +1332,11 @@ def combine_fast_slow(wl, flat_2d, flat_dq, tab_wl, tab_flat, dispaxis):
     else:
         combined_dq = flat_dq.copy()
 
+    if flat_err is None:
+        combined_err = np.zeros(wl.shape, dtype=np.float64)
+    else:
+        combined_err = flat_err.copy()
+
     if dispaxis == HORIZONTAL:
         dwl[:, 0:-1] = wl_c[:, 1:] - wl_c[:, 0:-1]
         dwl[:, -1] = dwl[:, -2]
@@ -1328,17 +1361,30 @@ def combine_fast_slow(wl, flat_2d, flat_dq, tab_wl, tab_flat, dispaxis):
         values += weight * np.interp(wavelengths, tab_wl, tab_flat,
                                      left=np.nan, right=np.nan)
 
+    # Interpolate error values from reference file using a simple
+    # linear interpolation as these don't have the required precision
+    # to justify a more complex interpolation
+    error_value = np.interp(wl_c, tab_wl, tab_flat_error, left=np.nan, right=np.nan)
+
     # Handle bad wavelength values in un-cleaned wavelength array
     bad = (wl <= 0)
     values[bad] = 1.0
+    error_value[bad] = 0.0
 
     # Handle missing values
     missing = np.isnan(values)
     values[missing] = 1.0
+    error_value[missing] = 0.0
     combined_dq[missing] |= dqflags.pixel['NO_FLAT_FIELD']
     combined_dq[missing] |= dqflags.pixel['DO_NOT_USE']
 
-    return flat_2d * values, combined_dq
+    # Add new 1D errors and input 2D errors in quadrature,
+    # treating NaNs as zeros
+    v1 = np.square(error_value * flat_2d)
+    v2 = np.square(combined_err * values)
+    combined_err = np.sqrt(np.nansum([v1, v2], axis=0))
+
+    return flat_2d * values, combined_dq, combined_err
 
 
 def clean_wl(wl, dispaxis):


Test patch:
diff --git a/jwst/flatfield/tests/test_combine_fast_slow.py b/jwst/flatfield/tests/test_combine_fast_slow.py
index 022b164b8..f23885fa5 100644
--- a/jwst/flatfield/tests/test_combine_fast_slow.py
+++ b/jwst/flatfield/tests/test_combine_fast_slow.py
@@ -2,6 +2,7 @@
 Test for flat_field.combine_fast_slow
 """
 import numpy as np
+import pytest
 from scipy.integrate import quad
 from astropy.modeling import polynomial
 from stdatamodels.jwst.datamodels import dqflags
@@ -9,7 +10,12 @@ from stdatamodels.jwst.datamodels import dqflags
 from jwst.flatfield.flat_field import combine_fast_slow
 
 
-def test_combine_fast_slow():
+@pytest.mark.parametrize('flat_err_1,flat_err_2',
+                         [(0.0, None), (np.nan, None), (np.nan, np.nan),
+                          (0.0, 0.1), (np.nan, 0.1), (0, np.nan),
+                          (0.1, 0.0), (0.1, np.nan),
+                          (0.1, None), (0.1, 0.1)])
+def test_combine_fast_slow(flat_err_1, flat_err_2):
     """Assign an array of values to tab_wl; these are strictly increasing,
     but they're not uniformly spaced.
     The abscissas and weights in combine_fast_slow are for three-point
@@ -34,6 +40,7 @@ def test_combine_fast_slow():
              'c5': 2.0e-3}
     poly = polynomial.Polynomial1D(degree=5, **coeff)
     tab_flat = poly(tab_wl)
+    tab_flat_err = np.full(tab_flat.shape, flat_err_1, dtype=tab_flat.dtype)
 
     wl0 = 7.37
     dwl0 = 2.99
@@ -50,21 +57,47 @@ def test_combine_fast_slow():
     wl = np.tile(wl_1d, (10, 1))
     flat_2d = np.ones((10, 10), dtype=float)
     flat_dq = np.zeros((10, 10), dtype=np.uint32)
+
+    # input 2d flat error
+    if flat_err_2 is None:
+        flat_err = None
+        v2 = 0
+    else:
+        flat_err = np.full((10, 10), flat_err_2)
+        v2 = flat_err_2**2
+
+    # Expected output variance is:
+    #   v = v1 * f2^2 + v2 * f1^2, where f2 is always 1
+    err_result = np.sqrt(np.nansum([flat_err_1**2, v2 * correct_value**2]))
+
     dispaxis = 1
+    value, new_dq, new_err = combine_fast_slow(wl, flat_2d, flat_dq, flat_err, tab_wl, tab_flat, tab_flat_err, dispaxis)
 
-    value, new_dq = combine_fast_slow(wl, flat_2d, flat_dq, tab_wl,
-                                      tab_flat, dispaxis)
+    # Check output shapes
+    assert value.shape == wl.shape
+    assert new_dq.shape == wl.shape
+    assert new_err.shape == wl.shape
 
     # Column 5 is the expected value
     assert np.allclose(value[:, 5], correct_value, rtol=1.e-8, atol=1.e-8)
     assert np.all(new_dq[:, 5] == 0)
+    assert np.allclose(new_err[:, 5], err_result, rtol=1.e-8, atol=1.e-8)
 
     # Columns 0-2 are bad (negative wavelengths, not marked in DQ)
     assert np.all(value[:, :3] == 1)
     assert np.all(new_dq[:, :3] == 0)
+    if flat_err is None or np.isnan(flat_err_2):
+        assert np.all(new_err[:, :3] == 0)
+    else:
+        assert np.allclose(new_err[:, :3], flat_err[:, :3], equal_nan=True)
 
     # Columns 3, 4, 6-9 are not covered by the tabular data
     # (missing values, marked in DQ)
     bad_value = dqflags.pixel['NO_FLAT_FIELD'] | dqflags.pixel['DO_NOT_USE']
     assert np.all(value[:, (3, 4, 6, 7, 8, 9)] == 1)
     assert np.all(new_dq[:, (3, 4, 6, 7, 8, 9)] == bad_value)
+    if flat_err is None or np.isnan(flat_err_2):
+        assert np.all(new_err[:, (3, 4, 6, 7, 8, 9)] == 0)
+    else:
+        assert np.allclose(new_err[:, (3, 4, 6, 7, 8, 9)],
+                           flat_err[:, (3, 4, 6, 7, 8, 9)])