Merge pull request #125 from ojustino/bkg-overlap-fix

Found cause of background "overlap" error

CHANGES.rst

@@ -4,8 +4,9 @@
 Bug Fixes
 ^^^^^^^^^
 
-- Improved errors/warnings when background region extends beyond bounds of image. [#127]
-
+- Improved errors/warnings when background region extends beyond bounds of image [#127]
+- Fixed boxcar weighting bug that often resulted in peak pixels having weight
+  above 1 and erroneously triggered overlapping background errors [#125]
 
 1.1.0
 -----

specreduce/extract.py

@@ -19,24 +19,46 @@
 def _get_boxcar_weights(center, hwidth, npix):
     """
     Compute weights given an aperture center, half width,
-    and number of pixels
+    and number of pixels.
+
+    Based on `get_boxcar_weights()` from a JDAT Notebook by Karl Gordon:
+    https://github.com/spacetelescope/jdat_notebooks/blob/main/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb
+
+    Parameters
+    ----------
+    center : float, required
+        The index of the aperture's center pixel on the larger image's
+        cross-dispersion axis.
+
+    hwidth : float, required
+        Half of the aperture's width in the cross-dispersion direction.
+
+    npix : float, required
+        The number of pixels in the larger image's cross-dispersion
+        axis.
+
+    Returns
+    -------
+    weights : `~numpy.ndarray`
+        A 2D image with weights assigned to pixels that fall within the
+        defined aperture.
     """
-    weights = np.zeros((npix))
+    weights = np.zeros(npix)
+
+    # shift center from integer to pixel space, where pixel N is [N-0.5, N+0.5),
+    # not [N, N+1). a pixel's integer index corresponds to its middle, not edge
+    center += 0.5
 
-    # pixels with full weight
-    fullpixels = [max(0, int(center - hwidth + 1)),
-                  min(int(center + hwidth), npix)]
-    weights[fullpixels[0]:fullpixels[1]] = 1.0
+    # pixels given full weight because they sit entirely within the aperture
+    fullpixels = [max(0, int(np.ceil(center - hwidth))),
+                  min(int(np.floor(center + hwidth)), npix)]
+    weights[fullpixels[0]:fullpixels[1]] = 1
 
-    # pixels at the edges of the boxcar with partial weight
+    # pixels at the edges of the boxcar with partial weight, if any
     if fullpixels[0] > 0:
-        w = hwidth - (center - fullpixels[0] + 0.5)
-        if w >= 0:
-            weights[fullpixels[0] - 1] = w
-        else:
-            weights[fullpixels[0]] = 1. + w
+        weights[fullpixels[0] - 1] = hwidth - (center - fullpixels[0])
     if fullpixels[1] < npix:
-        weights[fullpixels[1]] = hwidth - (fullpixels[1] - center - 0.5)
+        weights[fullpixels[1]] = hwidth - (fullpixels[1] - center)
 
     return weights
 
@@ -46,23 +68,26 @@ def _ap_weight_image(trace, width, disp_axis, crossdisp_axis, image_shape):
     """
     Create a weight image that defines the desired extraction aperture.
 
+    Based on `ap_weight_images()` from a JDAT Notebook by Karl Gordon:
+    https://github.com/spacetelescope/jdat_notebooks/blob/main/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb
+
     Parameters
     ----------
-    trace : Trace
+    trace : `~specreduce.tracing.Trace`, required
         trace object
-    width : float
+    width : float, required
         width of extraction aperture in pixels
-    disp_axis : int
+    disp_axis : int, required
         dispersion axis
-    crossdisp_axis : int
+    crossdisp_axis : int, required
         cross-dispersion axis
-    image_shape : tuple with 2 elements
+    image_shape : tuple with 2 elements, required
         size (shape) of image
 
     Returns
     -------
-    wimage : 2D image
-        weight image defining the aperture
+    wimage : `~numpy.ndarray`
+        a 2D weight image defining the aperture
     """
     wimage = np.zeros(image_shape)
     hwidth = 0.5 * width

specreduce/tests/test_extract.py

@@ -56,7 +56,7 @@ def test_boxcar_extraction():
 
     trace.set_position(14.3)
     spectrum = boxcar(width=4.7)
-    assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.0))
+    assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.15))
 
 
 def test_boxcar_outside_image_condition():

specreduce/tests/test_tracing.py

@@ -100,8 +100,8 @@ def test_kosmos_trace():
     tg = KosmosTrace(img, peak_method='gaussian')
     tc = KosmosTrace(img, peak_method='centroid')
     tm = KosmosTrace(img, peak_method='max')
-    # traces should all be close to 100 (note this passes with the current set values, but if
-    # changing the seed, noise, etc, then this test might need to be updated)
+    # traces should all be close to 100
+    # (values may need to be updated on changes to seed, noise, etc.)
     assert np.max(abs(tg.trace-100)) < sigma_pix
     assert np.max(abs(tc.trace-100)) < 3 * sigma_pix
     assert np.max(abs(tm.trace-100)) < 6 * sigma_pix
@@ -116,6 +116,14 @@ def test_kosmos_trace():
     img_win_nans = img.copy()
     img_win_nans[guess - window:guess + window] = np.nan
 
+    # ensure a low bin number is rejected
+    with pytest.raises(ValueError, match='bins must be >= 4'):
+        KosmosTrace(img, bins=3)
+
+    # ensure number of bins greater than number of dispersion pixels is rejected
+    with pytest.raises(ValueError, match=r'bins must be <*'):
+        KosmosTrace(img, bins=ncols)
+
     # error on trace of otherwise valid image with all-nan window around guess
     try:
         KosmosTrace(img_win_nans, guess=guess, window=window)