Repository: spacetelescope/jwst

Issue:
Revise how resample handles VAR arrays
_Issue [JP-3569](https://jira.stsci.edu/browse/JP-3569) was created on JIRA by [David Law](https://jira.stsci.edu/secure/ViewProfile.jspa?name=dlaw):_

At the moment, the resample step (for both imaging and spectroscopy) is using inverse variance weighting to combine the variance arrays (RNOISE, POISSON, FLAT).  It's doing so by first drizzling each variance from a single exposure to the output frame using uniform weighting, and then doing an inverse-variance weighted combination across all of the exposures to be combined together.  Crucially though, these arrays use themselves as the inverse variance weights.  I.e., for VAR_FLAT, 1/VAR_FLAT is being used as the weighting to combine together the individual exposures to create the final VAR_FLAT.  Likewise, VAR_POISSON is being used as weights for creating the final VAR_POISSON.

In general, for the weighted combination of two images

f3 = (w1 * f1 + w2 * x2)/(w1 + w2)

and the variance V3 of the output image is

V3 = (w1*w1*V1 + w2*w2*V2)/((w1+w2)*(w1+w2))

It's only in the special case when weighting by the variance of these inputs that this reduces to

V3 = 1/( (1/V1) + (1/V2))

Generally we don't want to do that as the weights should always be derived from the read noise variance, not whichever variance is being propagated.

We should therefore modify resample.py in two ways:

1) When using ivm weighting, resample_variance_array should be sure to use the read noise variance to weight the combination of the drizzled variances (Poisson, RN, and Flat) from different exposures.

2) When using exposure time weighting, resample_variance_array should actually use the exposure times (e.g., output_model.meta.exposure.duration) for weighting rather doing an inverse variance weight.

This should address the outstanding issue of zero-valued flatfield variances, and allow such values in reference files to propagate as expected to output products.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437?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 `81.92771%` with `15 lines` in your changes are missing coverage. Please review.
> Project coverage is 56.28%. Comparing base [(`2fb073e`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/2fb073e828d06f18ef7bd163445007c463b45136?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to head [(`1cee9f8`)](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> Report is 38 commits behind head on master.

| [Files](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437?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/resample/resample.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437?src=pr&el=tree&filepath=jwst%2Fresample%2Fresample.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9yZXNhbXBsZS9yZXNhbXBsZS5weQ==) | 81.48% | [15 Missing :warning: ](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437?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    #8437       +/-   ##
===========================================
- Coverage   75.31%   56.28%   -19.03%     
===========================================
  Files         474      390       -84     
  Lines       38965    38920       -45     
===========================================
- Hits        29345    21907     -7438     
- Misses       9620    17013     +7393     
```

| [Flag](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | |
|---|---|---|
| [nightly](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | `?` | |


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/8437?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).

hbushouse: @melanieclarke Pro tip: If you include the name of the JP ticket in the PR title, it'll automatically create links between the two. e.g. "JP-3569: Update variance ..."
melanieclarke: I think this is ready for review now, if someone can start a regression test run for me.
drlaw1558: LG2M from the science side; thanks for adding the doc update.
melanieclarke: Regression tests started here:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1413/
melanieclarke: Reviewing the regression tests, it looks like all changes are to post-resampling products (i2d, s2d, x1d, or catalog), in the error or variance arrays only as expected:
- In the ERR array in s2d and i2d files, there are now NaNs where there used to be 0s around the edges of the images. This is because I added some handling to catch when all input components are NaN and return NaN instead of 0 when they are added with np.nansum.
- In the variance arrays in s2d and i2d images, away from the edges, there are now 0s where there used to be NaNs.  This is due to the change in weighting in the variance combination: since inverse variances are no longer used in the propagation, 0 is an allowed value.
- In the error arrays in the x1d images, there are now values where there used to be NaNs.  This is due to the above change allowing 0 values.
- Most variance and error values are a little smaller, due to the change in the drizzling method for the variance components.

One test result I don't understand and will look into a little closer:
- The catalog output error columns contain some new NaNs where they used to have values, for both MIRI and NIRCam.

The NIRISS filter rotation failure in assign_wcs looks unrelated.
melanieclarke: It looks like the catalog NaNs correspond to sources with saturated cores, for which the variance components are all NaN in the i2d file.  The ERR is now set to NaN at these locations, instead of 0.  The aperture photometry in the source_catalog does not ignore NaNs, which is why NaN values now appear in the error columns in the output catalog.   

They do not also appear in the  flux columns because fillval=INDEF in resample, by default, so missing values are set to zero in the SCI data.  If the fillval default is set to NaN, saturated stars will also report flux = NaN.

@hbushouse - Should I add a mask for NaN values to the source_catalog aperture photometry, leave it as is, or file another ticket?
melanieclarke: Per @drlaw1558 and the JP coordination meeting today - the issue with NaN values in aperture photometry needs a separate ticket.  I will leave the code as is in the catalog step.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 7c45d080d..3f7018eb6 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -44,6 +44,10 @@ extract_1d
 
 - Replaced deprecated ``np.trapz`` with ``np.trapezoid()``. [#8415]
 
+- Fix a crash in ``extract_1d`` encountered when multiple background or source
+  regions are specified and the lower and upper limits for one of them are
+  outside the valid area for some data range. [#8433]
+
 flat_field
 ----------
 
@@ -92,6 +96,9 @@ resample
 - Remove sleep in median combination added in 8305 as it did not address
   the issue in operation [#8419]
 
+- Update variance handling to propagate resampled variance components with
+  weights that match the science `weight_type`. [#8437]
+
 resample_spec
 -------------
 
diff --git a/docs/jwst/resample/main.rst b/docs/jwst/resample/main.rst
index e0008f0cd..08e0069f2 100644
--- a/docs/jwst/resample/main.rst
+++ b/docs/jwst/resample/main.rst
@@ -34,6 +34,60 @@ WCS of each input image and the WCS of the defined output product.
 This mapping function gets passed to cdriz to drive the actual
 drizzling to create the output product.
 
+
+Error Propagation
+-----------------
+
+The error associated with each resampled pixel can in principle be derived
+from the variance components associated with each input pixel, weighted by
+the square of the input user weights and the square of the overlap between
+the input and output pixels. In practice, the cdriz routine does not currently
+support propagating variance data alongside science images, so the output
+error cannot be precisely calculated.
+
+To approximate the error on a resampled pixel, the variance arrays associated
+with each input image are resampled individually, then combined with a weighted
+sum.  The process is:
+
+#. For each input image, take the square root of each of the read noise variance
+   array to make an error image.
+
+#. Drizzle the read noise error image onto the output WCS, with drizzle
+   parameters matching those used for the science data.
+
+#. Square the resampled read noise to make a variance array.
+
+   a. If the resampling `weight_type` is an inverse variance map (`ivm`), weight
+      the resampled variance by the square of its own inverse.
+
+   #. If the `weight_type` is the exposure time (`exptime`), weight the
+      resampled variance by the square of the exposure time for the image.
+
+#. Add the weighted, resampled read noise variance to a running sum across all
+   images.  Add the weights (unsquared) to a separate running sum across
+   all images.
+
+#. Perform the same steps for the Poisson noise variance and the flat variance.
+   For these components, the weight for the sum is either the resampled read
+   noise variance or else the exposure time.
+
+#. For each variance component (read noise, Poisson, and flat), divide the
+   weighted variance sum by the total weight, squared.
+
+After each variance component is resampled and summed, the final error
+array is computed as the square root of the sum of the three independent
+variance components.  This error image is stored in the ``err`` attribute
+in the output data model or the ``'ERR'`` extension of the FITS file.
+
+It is expected that the output errors computed in this way will
+generally overestimate the true error on the resampled data.  The magnitude
+of the overestimation depends on the details of the pixel weights
+and error images.  Note, however, that drizzling error images produces
+a much better estimate of the output error than directly drizzling
+the variance images, since the kernel overlap weights do not need to be
+squared for combining error values.
+
+
 Context Image
 -------------
 
@@ -93,6 +147,7 @@ context array ``con``, one can do something like this:
 For convenience, this functionality was implemented in the
 :py:func:`~jwst.resample.resample_utils.decode_context` function.
 
+
 Spectroscopic Data
 ------------------
 
diff --git a/jwst/extract_1d/extract1d.py b/jwst/extract_1d/extract1d.py
index 83400a023..979a3c00e 100644
--- a/jwst/extract_1d/extract1d.py
+++ b/jwst/extract_1d/extract1d.py
@@ -676,10 +676,17 @@ def _extract_colpix(image_data, x, j, limits):
     # compute number of data points:
     ns = image_data.shape[0] - 1
     ns12 = ns + 0.5
-
-    npts = sum(map(lambda x: min(ns, int(math.floor(x[1] + 0.5))) -
-                   max(0, int(math.floor(x[0] + 0.5))) + 1, intervals))
-    npts = max(npts, 1)
+    npts = 0
+    for interval in intervals:
+        if interval[0] == interval[1]:
+            # no data between limits
+            continue
+        maxval = min(ns, int(math.floor(interval[1] + 0.5)))
+        minval = max(0, int(math.floor(interval[0] + 0.5)))
+        if maxval - minval + 1 > 0:
+            npts += maxval - minval + 1
+    if npts == 0:
+        return [], [], []
 
     # pre-allocate data arrays:
     y = np.empty(npts, dtype=np.float64)
@@ -689,6 +696,8 @@ def _extract_colpix(image_data, x, j, limits):
     # populate data and weights:
     k = 0
     for i in intervals:
+        if i[0] == i[1]:
+            continue
         i1 = i[0] if i[0] >= -0.5 else -0.5
         i2 = i[1] if i[1] <= ns12 else ns12
 
@@ -697,6 +706,7 @@ def _extract_colpix(image_data, x, j, limits):
         ii2 = min(ns, int(math.floor(i2 + 0.5)))
 
         # special case: ii1 == ii2:
+        # take the value at the pixel
         if ii1 == ii2:
             v = image_data[ii1, x]
             val[k] = v
diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py
index 19a38489f..e33b5305a 100644
--- a/jwst/resample/resample.py
+++ b/jwst/resample/resample.py
@@ -1,5 +1,6 @@
 import logging
 import os
+import warnings
 
 import numpy as np
 from drizzle import util
@@ -373,20 +374,19 @@ class ResampleData:
             )
             del data, inwht
 
-        # Resample variances array in self.input_models to output_model
-        self.resample_variance_array("var_rnoise", output_model)
-        self.resample_variance_array("var_poisson", output_model)
-        self.resample_variance_array("var_flat", output_model)
-        output_model.err = np.sqrt(
-            np.nansum(
-                [
-                    output_model.var_rnoise,
-                    output_model.var_poisson,
-                    output_model.var_flat
-                ],
-                axis=0
-            )
-        )
+        # Resample variance arrays in self.input_models to output_model
+        self.resample_variance_arrays(output_model)
+        var_components = [
+            output_model.var_rnoise,
+            output_model.var_poisson,
+            output_model.var_flat
+        ]
+        output_model.err = np.sqrt(np.nansum(var_components,axis=0))
+
+        # nansum returns zero for input that is all NaN -
+        # set those values to NaN instead
+        all_nan = np.all(np.isnan(var_components), axis=0)
+        output_model.err[all_nan] = np.nan
 
         self.update_exposure_times(output_model)
         self.output_models.append(output_model)
@@ -396,85 +396,160 @@ class ResampleData:
 
         return self.output_models
 
-    def resample_variance_array(self, name, output_model):
-        """Resample variance arrays from self.input_models to the output_model
+    def resample_variance_arrays(self, output_model):
+        """Resample variance arrays from self.input_models to the output_model.
 
-        Resample the ``name`` variance array to the same name in output_model,
-        using a cumulative sum.
+        Variance images from each input model are resampled individually and
+        added to a weighted sum. If weight_type is 'ivm', the inverse of the
+        resampled read noise variance is used as the weight for all the variance
+        components. If weight_type is 'exptime', the exposure time is used.
 
-        This modifies output_model in-place.
+        The output_model is modified in place.
         """
-        output_wcs = output_model.meta.wcs
-        inverse_variance_sum = np.full_like(output_model.data, np.nan)
-
-        log.info(f"Resampling {name}")
+        log.info("Resampling variance components")
+        weighted_rn_var = np.full_like(output_model.data, np.nan)
+        weighted_pn_var = np.full_like(output_model.data, np.nan)
+        weighted_flat_var = np.full_like(output_model.data, np.nan)
+        total_weight_rn_var = np.zeros_like(output_model.data)
+        total_weight_pn_var = np.zeros_like(output_model.data)
+        total_weight_flat_var = np.zeros_like(output_model.data)
         for model in self.input_models:
-            variance = getattr(model, name)
-            if variance is None or variance.size == 0:
-                log.debug(
-                    f"No data for '{name}' for model "
-                    f"{repr(model.meta.filename)}. Skipping ..."
+            # Do the read noise variance first, so it can be
+            # used for weights if needed
+            rn_var = self._resample_one_variance_array(
+                "var_rnoise", model, output_model)
+
+            # Find valid weighting values in the variance
+            if rn_var is not None:
+                mask = (rn_var > 0) & np.isfinite(rn_var)
+            else:
+                mask = np.full_like(rn_var, False)
+
+            # Set the weight for the image from the weight type
+            weight = np.ones(output_model.data.shape)
+            if self.weight_type == "ivm" and rn_var is not None:
+                weight[mask] = rn_var[mask] ** -1
+            elif self.weight_type == "exptime":
+                if resample_utils.check_for_tmeasure(model):
+                    weight[:] = model.meta.exposure.measurement_time
+                else:
+                    weight[:] = model.meta.exposure.exposure_time
+
+            # Weight and add the readnoise variance
+            # Note: floating point overflow is an issue if variance weights
+            # are used - it can't be squared before multiplication
+            if rn_var is not None:
+                mask = (rn_var >= 0) & np.isfinite(rn_var) & (weight > 0)
+                weighted_rn_var[mask] = np.nansum(
+                    [weighted_rn_var[mask],
+                     rn_var[mask] * weight[mask] * weight[mask]],
+                    axis=0
                 )
-                continue
-
-            elif variance.shape != model.data.shape:
-                log.warning(
-                    f"Data shape mismatch for '{name}' for model "
-                    f"{repr(model.meta.filename)}. Skipping ..."
+                total_weight_rn_var[mask] += weight[mask]
+
+            # Now do poisson and flat variance, updating only valid new values
+            # (zero is a valid value; negative, inf, or NaN are not)
+            pn_var = self._resample_one_variance_array(
+                "var_poisson", model, output_model)
+            if pn_var is not None:
+                mask = (pn_var >= 0) & np.isfinite(pn_var) & (weight > 0)
+                weighted_pn_var[mask] = np.nansum(
+                    [weighted_pn_var[mask],
+                     pn_var[mask] * weight[mask] * weight[mask]],
+                    axis=0
                 )
-                continue
-
-            # Make input weight map of unity where there is science data
-            inwht = resample_utils.build_driz_weight(
-                model,
-                weight_type=None,
-                good_bits=self.good_bits
-            )
-
-            resampled_variance = np.zeros_like(output_model.data)
-            outwht = np.zeros_like(output_model.data)
-            outcon = np.zeros_like(output_model.con)
-
-            xmin, xmax, ymin, ymax = resample_utils._resample_range(
-                variance.shape,
-                model.meta.wcs.bounding_box
+                total_weight_pn_var[mask] += weight[mask]
+
+            flat_var = self._resample_one_variance_array(
+                "var_flat", model, output_model)
+            if flat_var is not None:
+                mask = (flat_var >= 0) & np.isfinite(flat_var) & (weight > 0)
+                weighted_flat_var[mask] = np.nansum(
+                    [weighted_flat_var[mask],
+                     flat_var[mask] * weight[mask] * weight[mask]],
+                    axis=0
+                )
+                total_weight_flat_var[mask] += weight[mask]
+
+        # We now have a sum of the weighted resampled variances.
+        # Divide by the total weights, squared, and set in the output model.
+        # Zero weight and missing values are NaN in the output.
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", "invalid value*", RuntimeWarning)
+            warnings.filterwarnings("ignore", "divide by zero*", RuntimeWarning)
+
+            output_variance = (weighted_rn_var
+                               / total_weight_rn_var / total_weight_rn_var)
+            setattr(output_model, "var_rnoise", output_variance)
+
+            output_variance = (weighted_pn_var
+                               / total_weight_pn_var / total_weight_pn_var)
+            setattr(output_model, "var_poisson", output_variance)
+
+            output_variance = (weighted_flat_var
+                               / total_weight_flat_var / total_weight_flat_var)
+            setattr(output_model, "var_flat", output_variance)
+
+    def _resample_one_variance_array(self, name, input_model, output_model):
+        """Resample one variance image from an input model.
+
+        The error image is passed to drizzle instead of the variance, to
+        better match kernel overlap and user weights to the data, in the
+        pixel averaging process. The drizzled error image is squared before
+        returning.
+        """
+        variance = getattr(input_model, name)
+        if variance is None or variance.size == 0:
+            log.debug(
+                f"No data for '{name}' for model "
+                f"{repr(input_model.meta.filename)}. Skipping ..."
             )
+            return
 
-            iscale = model.meta.iscale
-
-            # Resample the variance array. Fill "unpopulated" pixels with NaNs.
-            self.drizzle_arrays(
-                variance,
-                inwht,
-                model.meta.wcs,
-                output_wcs,
-                resampled_variance,
-                outwht,
-                outcon,
-                iscale=iscale**2,
-                pixfrac=self.pixfrac,
-                kernel=self.kernel,
-                fillval=np.nan,
-                xmin=xmin,
-                xmax=xmax,
-                ymin=ymin,
-                ymax=ymax
+        elif variance.shape != input_model.data.shape:
+            log.warning(
+                f"Data shape mismatch for '{name}' for model "
+                f"{repr(input_model.meta.filename)}. Skipping ..."
             )
+            return
 
-            # Add the inverse of the resampled variance to a running sum.
-            # Update only pixels (in the running sum) with valid new values:
-            mask = resampled_variance > 0
+        # Make input weight map
+        inwht = resample_utils.build_driz_weight(
+            input_model,
+            weight_type=self.weight_type,  # weights match science
+            good_bits=self.good_bits
+        )
 
-            inverse_variance_sum[mask] = np.nansum(
-                [inverse_variance_sum[mask], np.reciprocal(resampled_variance[mask])],
-                axis=0
-            )
+        resampled_error = np.zeros_like(output_model.data)
+        outwht = np.zeros_like(output_model.data)
+        outcon = np.zeros_like(output_model.con)
 
-        # We now have a sum of the inverse resampled variances.  We need the
-        # inverse of that to get back to units of variance.
-        output_variance = np.reciprocal(inverse_variance_sum)
+        xmin, xmax, ymin, ymax = resample_utils._resample_range(
+            variance.shape,
+            input_model.meta.wcs.bounding_box
+        )
 
-        setattr(output_model, name, output_variance)
+        iscale = input_model.meta.iscale
+
+        # Resample the error array. Fill "unpopulated" pixels with NaNs.
+        self.drizzle_arrays(
+            np.sqrt(variance),
+            inwht,
+            input_model.meta.wcs,
+            output_model.meta.wcs,
+            resampled_error,
+            outwht,
+            outcon,
+            iscale=iscale,
+            pixfrac=self.pixfrac,
+            kernel=self.kernel,
+            fillval=np.nan,
+            xmin=xmin,
+            xmax=xmax,
+            ymin=ymin,
+            ymax=ymax
+        )
+        return resampled_error ** 2
 
     def update_exposure_times(self, output_model):
         """Modify exposure time metadata in-place"""
@@ -485,7 +560,7 @@ class ResampleData:
         measurement_time_failures = []
         for exposure in self.input_models.models_grouped:
             total_exposure_time += exposure[0].meta.exposure.exposure_time
-            if not resample_utils._check_for_tmeasure(exposure[0]):
+            if not resample_utils.check_for_tmeasure(exposure[0]):
                 measurement_time_failures.append(1)
             else:
                 total_measurement_time += exposure[0].meta.exposure.measurement_time
diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py
index 08572108d..ed4447d8b 100644
--- a/jwst/resample/resample_utils.py
+++ b/jwst/resample/resample_utils.py
@@ -182,7 +182,7 @@ def build_driz_weight(model, weight_type=None, good_bits=None):
             inv_variance = 1.0
         result = inv_variance * dqmask
     elif weight_type == 'exptime':
-        if _check_for_tmeasure(model):
+        if check_for_tmeasure(model):
             exptime = model.meta.exposure.measurement_time
         else:
             exptime = model.meta.exposure.exposure_time
@@ -311,7 +311,7 @@ def _resample_range(data_shape, bbox=None):
     return xmin, xmax, ymin, ymax
 
 
-def _check_for_tmeasure(model):
+def check_for_tmeasure(model):
     '''
     Check if the measurement_time keyword is present in the datamodel
     for use in exptime weighting. If not, revert to using exposure_time.


Test patch:
diff --git a/jwst/extract_1d/tests/test_extract_src_flux.py b/jwst/extract_1d/tests/test_extract_src_flux.py
index e3952745f..cfab06f55 100644
--- a/jwst/extract_1d/tests/test_extract_src_flux.py
+++ b/jwst/extract_1d/tests/test_extract_src_flux.py
@@ -4,12 +4,13 @@ Test for extract_1d._extract_src_flux
 import math
 
 import numpy as np
+import pytest
 
 from jwst.extract_1d import extract1d
 
 
-def test_extract_src_flux():
-
+@pytest.fixture
+def inputs_constant():
     shape = (9, 5)
     image = np.arange(shape[0] * shape[1], dtype=np.float32).reshape(shape)
     var_poisson = image.copy()
@@ -24,6 +25,14 @@ def test_extract_src_flux():
     upper = np.zeros(shape[1], dtype=np.float64) + 7.   # middle of pixel 7
     srclim = [[lower, upper]]
 
+    return (image, var_poisson, var_rnoise, var_rflat,
+            x, j, lam, srclim, weights, bkgmodels)
+
+
+def test_extract_src_flux(inputs_constant):
+    (image, var_poisson, var_rnoise, var_rflat,
+     x, j, lam, srclim, weights, bkgmodels) = inputs_constant
+
     (total_flux, f_var_poisson, f_var_rnoise, f_var_flat,
      bkg_flux, b_var_poisson, b_var_rnoise, b_var_flat,
      tarea, twht) = extract1d._extract_src_flux(
@@ -61,3 +70,48 @@ def test_extract_src_flux():
     assert np.isnan(total_flux)
 
     assert tarea == 0.
+
+
+@pytest.mark.parametrize('test_type', ['all_empty', 'all_equal'])
+def test_extract_src_flux_empty_interval(inputs_constant, test_type):
+    (image, var_poisson, var_rnoise, var_rflat,
+     x, j, lam, srclim, weights, bkgmodels) = inputs_constant
+
+    if test_type == 'all_empty':
+        # no limits provided
+        srclim = []
+    else:
+        # empty extraction range: upper equals lower
+        srclim[0][1] = srclim[0][0].copy()
+
+    (total_flux, f_var_poisson, f_var_rnoise, f_var_flat,
+     bkg_flux, b_var_poisson, b_var_rnoise, b_var_flat,
+     tarea, twht) = extract1d._extract_src_flux(
+        image, var_poisson, var_rnoise, var_rflat,
+        x, j, lam, srclim, weights, bkgmodels)
+
+    # empty interval, so no flux returned
+    assert np.isnan(total_flux)
+    assert bkg_flux == 0.
+    assert tarea == 0.
+
+
+@pytest.mark.parametrize('offset', [-100, 100])
+def test_extract_src_flux_interval_out_of_range(inputs_constant, offset):
+    (image, var_poisson, var_rnoise, var_rflat,
+     x, j, lam, srclim, weights, bkgmodels) = inputs_constant
+
+    # extraction limits out of range
+    srclim[0][0] += offset
+    srclim[0][1] += offset
+
+    (total_flux, f_var_poisson, f_var_rnoise, f_var_flat,
+     bkg_flux, b_var_poisson, b_var_rnoise, b_var_flat,
+     tarea, twht) = extract1d._extract_src_flux(
+        image, var_poisson, var_rnoise, var_rflat,
+        x, j, lam, srclim, weights, bkgmodels)
+
+    # empty interval, so no flux returned
+    assert np.isnan(total_flux)
+    assert bkg_flux == 0.
+    assert tarea == 0.
diff --git a/jwst/extract_1d/tests/test_fit_background_model.py b/jwst/extract_1d/tests/test_fit_background_model.py
index 3758564f4..5aeb59695 100644
--- a/jwst/extract_1d/tests/test_fit_background_model.py
+++ b/jwst/extract_1d/tests/test_fit_background_model.py
@@ -2,6 +2,8 @@
 Test for extract_1d._fit_background_model
 """
 import math
+from copy import deepcopy
+
 import numpy as np
 import pytest
 
@@ -138,3 +140,89 @@ def test_handles_one_value(inputs_constant):
     assert math.isclose(b_var_flat_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
 
     assert npts == 0
+
+
+@pytest.mark.parametrize('test_type', ['all_empty', 'all_equal'])
+def test_handles_empty_interval(inputs_constant, test_type):
+    image, var_poisson, var_rnoise, var_rflat, x, j, bkglim, bkg_fit, bkg_order = inputs_constant
+
+    if test_type == 'all_empty':
+        # no limits provided
+        bkglim = []
+    else:
+        # empty extraction range: upper equals lower
+        bkglim[0][1] = bkglim[0][0].copy()
+
+    # No data available: background model is 0
+    (bkg_model, b_var_poisson_model, b_var_rnoise_model, b_var_flat_model, npts) = \
+        extract1d._fit_background_model(image, var_poisson, var_rnoise, var_rflat,
+                                        x, j, bkglim, bkg_fit, bkg_order)
+
+    assert math.isclose(bkg_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(bkg_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_poisson_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_poisson_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_rnoise_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_rnoise_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_flat_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_flat_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert npts == 0
+
+
+@pytest.mark.parametrize('offset', [-100, 100])
+def test_handles_interval_out_of_range(inputs_constant, offset):
+    image, var_poisson, var_rnoise, var_rflat, x, j, bkglim, bkg_fit, bkg_order = inputs_constant
+    bkglim[0][0] += offset
+    bkglim[0][1] += offset
+
+    # No data available: background model is 0
+    (bkg_model, b_var_poisson_model, b_var_rnoise_model, b_var_flat_model, npts) = \
+        extract1d._fit_background_model(image, var_poisson, var_rnoise, var_rflat,
+                                        x, j, bkglim, bkg_fit, bkg_order)
+
+    assert math.isclose(bkg_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(bkg_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_poisson_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_poisson_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_rnoise_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_rnoise_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_flat_model(0.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_flat_model(8.), 0.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert npts == 0
+
+
+def test_handles_one_empty_interval(inputs_constant):
+    image, var_poisson, var_rnoise, var_rflat, x, j, bkglim, bkg_fit, bkg_order = inputs_constant
+
+    # add an extra interval that is empty
+    bkglim.append(deepcopy(bkglim[0]))
+    bkglim[1][0] += 2.0
+    bkglim[1][1] = bkglim[1][0].copy()
+    print(bkglim)
+
+    # should ignore the second interval and return a valid answer for the first
+    (bkg_model, b_var_poisson_model, b_var_rnoise_model, b_var_flat_model, npts) = \
+        extract1d._fit_background_model(image, var_poisson, var_rnoise, var_rflat,
+                                        x, j, bkglim, bkg_fit, bkg_order)
+
+    assert math.isclose(bkg_model(0.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(bkg_model(8.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_poisson_model(0.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_poisson_model(8.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_rnoise_model(0.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_rnoise_model(8.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert math.isclose(b_var_flat_model(0.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+    assert math.isclose(b_var_flat_model(8.), 22.0, rel_tol=1.e-8, abs_tol=1.e-8)
+
+    assert npts == 2
diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py
index 08e1fd436..73f914181 100644
--- a/jwst/resample/tests/test_resample_step.py
+++ b/jwst/resample/tests/test_resample_step.py
@@ -457,8 +457,10 @@ def test_wcs_keywords(nircam_rate):
     assert result.meta.wcsinfo.vparity is None
 
 
-@pytest.mark.parametrize("n_images", [1, 2, 3, 9])
-def test_resample_variance(nircam_rate, n_images):
+@pytest.mark.parametrize("n_images,weight_type",
+                         [(1, 'ivm'), (2, 'ivm'), (3, 'ivm'), (9, 'ivm'),
+                          (1, 'exptime'), (2, 'exptime'), (3, 'exptime'), (9, 'exptime')])
+def test_resample_variance(nircam_rate, n_images, weight_type):
     """Test that resampled variance and error arrays are computed properly"""
     err = 0.02429
     var_rnoise = 0.00034
@@ -474,7 +476,7 @@ def test_resample_variance(nircam_rate, n_images):
     for n in range(n_images):
         c.append(im.copy())
 
-    result = ResampleStep.call(c, blendheaders=False)
+    result = ResampleStep.call(c, blendheaders=False, weight_type=weight_type)
 
     # Verify that the combined uncertainty goes as 1 / sqrt(N)
     assert_allclose(result.err[5:-5, 5:-5].mean(), err / np.sqrt(n_images), atol=1e-5)
@@ -493,7 +495,13 @@ def test_resample_undefined_variance(nircam_rate, shape):
     c = ModelContainer([im])
 
     with pytest.warns(RuntimeWarning, match="var_rnoise array not available"):
-        ResampleStep.call(c, blendheaders=False)
+        result = ResampleStep.call(c, blendheaders=False)
+
+    # no valid variance - output error and variance are all NaN
+    assert_allclose(result.err, np.nan)
+    assert_allclose(result.var_rnoise, np.nan)
+    assert_allclose(result.var_poisson, np.nan)
+    assert_allclose(result.var_flat, np.nan)
 
 
 @pytest.mark.parametrize('ratio', [0.7, 1.2])