Repository: spacetelescope/jwst

Issue:
Extraction background edge case bug
_Issue [JP-3601](https://jira.stsci.edu/browse/JP-3601) was created on JIRA by [Melanie Clarke](https://jira.stsci.edu/secure/ViewProfile.jspa?name=mclarke):_

While experimenting with setting polynomial background regions in extract_1d for BOTS data, I came across a crash when background upper and lower limits are both out of range for some part of the array.  This happens when setting a region above and below the curved spectral trace in a narrow subarray: there is data within the background limits for part of the array, but not the full array.  In this case, the limits are checked and both upper and lower limits are set to the edge of the array (0 at the bottom, ns at the top).

The crash comes from line 702 in jwst/extract_1d/extract1d.py:

IndexError: index 6 is out of bounds for axis 0 with size 6

I think the problem is that when counting intervals (npts, line 680), intervals with equal upper and lower limits are not counted, but when computing background levels, the value at that point is added to the val array (line 702).  I think the fix should be to skip intervals for which upper and lower limits are equal – it is not likely that that point contains a valid background level anyway.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8433?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 56.26%. 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 [(`db1747f`)](https://app.codecov.io/gh/spacetelescope/jwst/pull/8433?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> Report is 36 commits behind head on master.


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


```diff
@@             Coverage Diff             @@
##           master    #8433       +/-   ##
===========================================
- Coverage   75.31%   56.26%   -19.05%     
===========================================
  Files         474      390       -84     
  Lines       38965    38861      -104     
===========================================
- Hits        29345    21867     -7478     
- Misses       9620    16994     +7374     
```

| [Flag](https://app.codecov.io/gh/spacetelescope/jwst/pull/8433/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/8433/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/8433?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: In testing, the initial simple fix was a little too simple: it did not handle correctly the cases where part of a pixel falls within the specified region, or the case where an empty interval is specified, but the limits are within range of the data.  The updated version should now handle both of those cases correctly, as well as fixing the original crash.
melanieclarke: Thanks @nden - I'll see if I can add some unit tests to cover the relevant cases.
melanieclarke: @nden - unit tests added for a couple relevant cases in background and source extraction.

Looking at the regression tests, it looks like the test_nirspec_fs_spec3[s00024-s1600a1-x1d] failure is the only relevant change.  For this spectral extraction, the WCS is invalid: all wavelengths are set to NaN and the extraction region is out of range.   The old version defaulted the flux to 0.0 and the errors to NaN in this case; the new version defaults the flux to NaN and the errors to 0.0.  I think the new version is closer to the intended default values.

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 7c45d080d..ca6f6531e 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
 ----------
 
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


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