Merge pull request #245 from hpparvi/weighted_boxcar_extraction

Boxcar extraction with non-finite pixels

Author: hpparvi

CHANGES.rst

@@ -18,9 +18,16 @@ Other changes
 
 New Features
 ^^^^^^^^^^^^
+
 - Added the ``mask_treatment`` parameter to Background, Trace, and Boxcar Extract
-  operations to handle non-finite data and boolean masks. Choice of ``filter``,
-  ``omit``, or ``zero-fill``. [#216]
+  operations to handle non-finite data and boolean masks. Available options are
+  ``filter``, ``omit``, or ``zero-fill``, with ``exclude`` additionally available
+  for BoxcarExtract. [#216]
+
+- Modified BoxcarExtract to ignore non-finite pixels when ``mask_treatment`` is set
+  to ``exclude``; otherwise, non-finite values are propagated. Boxcar extraction is
+  now carried out as a weighed sum over the window. When no non-finite values are
+  present, the extracted spectra remain unchanged from the previous behaviour.
 
 API Changes
 ^^^^^^^^^^^

docs/extraction_quickstart.rst

@@ -88,10 +88,17 @@ then this will be used. Otherwise, the ``variance`` parameter must be set.::
   extract = specreduce.extract.HorneExtract(image-bg, trace, variance=var_array)
 
 An optional mask array for the image may be supplied to HorneExtract as well. 
-This follows the same convention and can either be attacted to ``image`` if it
-is and ``astropy.NDData`` object, or supplied as a keyword argrument. Note that
-any wavelengths columns containing any masked values will be omitted from the
-extraction.
+This follows the same convention and can either be attached to ``image`` if it
+is an ``astropy.NDData`` object, or supplied as a keyword argument.
+
+The extraction methods automatically detect non-finite pixels in the input
+image and combine them with the user-supplied mask to prevent them from biasing the
+extraction. In the boxcar extraction, the treatment of these pixels is controlled by
+the ``mask_treatment`` option. When set to ``exclude`` (the default), non-finite
+pixels within the extraction window are excluded from the extraction, and the extracted
+flux is scaled according to the effective number of unmasked pixels. When using other
+options (``filter`` or ``omit``), the non-finite values may be propagated or treated
+differently as documented in the API.
 
 The previous examples in this section show how to initialize the BoxcarExtract
 or HorneExtract objects with their required parameters. To extract the 1D
@@ -111,8 +118,8 @@ or, for example to override the original ``trace_object``::
 Spatial profile options
 -----------------------
 The Horne algorithm provides two options for fitting the spatial profile to the
-cross_dispersion direction of the source: a Gaussian fit (default),
-or an empirical 'interpolated_profile' option.
+cross dispersion direction of the source: a Gaussian fit (default),
+or an empirical ``interpolated_profile`` option.
 
 If the default Gaussian option is used, an optional background model may be
 supplied as well (default is a 2D Polynomial) to account

specreduce/extract.py

@@ -150,15 +150,15 @@ class BoxcarExtract(SpecreduceOperation):
         cross-dispersion axis
     mask_treatment
         The method for handling masked or non-finite data. Choice of ``filter``,
-        ``omit``, or ``zero-fill``. If `filter` is chosen, the mask is ignored
+        ``omit``, or ``exclude``. If ``filter`` is chosen, the mask is ignored
         and the non-finite data will passed to the extraction as is. If ``omit``
         is chosen, columns along disp_axis with any masked or non-finite data
         values will be fully masked (i.e, 2D mask is collapsed to 1D and applied).
-        If ``zero-fill`` is chosen, masked and non-finite data will be replaced
-        with 0.0 in the input image, and the mask will then be dropped.
-        For all three options, the input mask (optional on input NDData object)
-        will be combined with a mask generated from any non-finite values in the
-        image data.
+        If ``exclude`` is chosen, masked and non-finite data will be excluded
+        from the extraction and the spectrum extraction is carried out as a
+        weighted sum. For all three options, the input mask (optional on input
+        NDData object) will be combined with a mask generated from any non-finite
+        values in the image data.
 
     Returns
     -------
@@ -171,8 +171,9 @@ class BoxcarExtract(SpecreduceOperation):
     disp_axis: int = 1
     crossdisp_axis: int = 0
     # TODO: should disp_axis and crossdisp_axis be defined in the Trace object?
-    mask_treatment: Literal['filter', 'omit', 'zero-fill'] = 'filter'
-    _valid_mask_treatment_methods = ('filter', 'omit', 'zero-fill')
+    # TODO: should the 'filter' option be changed to 'propagate'
+    mask_treatment: Literal['filter', 'omit', 'exclude'] = 'exclude'
+    _valid_mask_treatment_methods = ('filter', 'omit', 'exclude')
 
     @property
     def spectrum(self):
@@ -182,61 +183,63 @@ def __call__(self, image: NDData | None = None,
                  trace: Trace | None = None,
                  width: float | None = None,
                  disp_axis: int | None = None,
-                 crossdisp_axis: int | None = None):
+                 crossdisp_axis: int | None = None) -> Spectrum1D:
         """
         Extract the 1D spectrum using the boxcar method.
 
         Parameters
         ----------
-        image : `~astropy.nddata.NDData`-like or array-like, required
-            image with 2-D spectral image data
-        trace : Trace, required
-            trace object
-        width : float, optional
-            width of extraction aperture in pixels [default: 5]
-        disp_axis : int, optional
-            dispersion axis [default: 1]
-        crossdisp_axis : int, optional
-            cross-dispersion axis [default: 0]
+        image
+            The image with 2-D spectral image data
+        trace
+            The trace object
+        width
+            The width of extraction aperture in pixels
+        disp_axis
+            The dispersion axis
+        crossdisp_axis
+            The cross-dispersion axis
 
         Returns
         -------
-        spec : `~specutils.Spectrum1D`
+        spec
             The extracted 1d spectrum with flux expressed in the same
             units as the input image, or u.DN, and pixel units
         """
         image = image if image is not None else self.image
-        trace = trace if trace is not None else self.trace_object
-        width = width if width is not None else self.width
-        disp_axis = disp_axis if disp_axis is not None else self.disp_axis
-        crossdisp_axis = crossdisp_axis if crossdisp_axis is not None else self.crossdisp_axis
+        trace = trace or self.trace_object
+        width = width or self.width
+        disp_axis = disp_axis or self.disp_axis
+        cdisp_axis = crossdisp_axis or self.crossdisp_axis
+        mask_mapping = {'filter': 'filter', 'exclude': 'filter', 'omit': 'omit'}
+
+        if width <= 0:
+            raise ValueError("The window width must be positive")
 
         # Parse image, including masked/nonfinite data handling based on
-        # choice of `mask_treatment`, which for BoxcarExtract can be filter, zero-fill, or
-        # omit. non-finite data will be masked, always. Returns a Spectrum1D.
+        # choice of `mask_treatment`, which for BoxcarExtract can be filter, omit, or
+        # exclude. non-finite data will be masked, always. Returns a Spectrum1D.
         self.image = self._parse_image(image, disp_axis=disp_axis,
-                                       mask_treatment=self.mask_treatment)
-
-        # # _parse_image returns a Spectrum1D. convert this to a masked array
-        # # for ease of calculations here (even if there is no masked data).
-        # img = np.ma.masked_array(self.image.data, self.image.mask)
+                                       mask_treatment=mask_mapping[self.mask_treatment])
+
+        # Spectrum extraction
+        # ===================
+        # Assign no weight to non-finite pixels outside the window. Non-finite pixels inside
+        # the window will be propagated to the sum if mask treatment is either ``filter`` or
+        # ``omit`` or excluded if the chosen mask treatment option is ``exclude``. In the
+        # latter case, the flux is calculated as the average of the non-masked pixels inside
+        # the window multiplied by the window width.
+        window_weights = _ap_weight_image(trace, width, disp_axis, cdisp_axis, self.image.shape)
+
+        if self.mask_treatment == 'exclude':
+            image_cleaned = np.where(~self.image.mask, self.image.data*window_weights, 0.0)
+            weights = np.where(~self.image.mask, window_weights, 0.0)
+            spectrum = image_cleaned.sum(axis=cdisp_axis) / weights.sum(axis=cdisp_axis) * width
+        else:
+            image_windowed = np.where(window_weights, self.image.data*window_weights, 0.0)
+            spectrum = np.sum(image_windowed, axis=cdisp_axis)
 
-        if width <= 0:
-            raise ValueError("width must be positive")
-
-        # weight image to use for extraction
-        wimg = _ap_weight_image(trace,
-                                width,
-                                disp_axis,
-                                crossdisp_axis,
-                                self.image.shape)
-
-        # extract, assigning no weight to non-finite pixels outside the window
-        # (non-finite pixels inside the window will still make it into the sum)
-        image_windowed = np.where(wimg, self.image.data*wimg, 0)
-        ext1d = np.sum(image_windowed, axis=crossdisp_axis)
-        return Spectrum1D(ext1d * self.image.unit,
-                          spectral_axis=self.image.spectral_axis)
+        return Spectrum1D(spectrum * self.image.unit, spectral_axis=self.image.spectral_axis)
 
 
 @dataclass

specreduce/tests/test_extract.py

@@ -83,14 +83,35 @@ def test_boxcar_extraction(mk_test_img):
     assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.15))
 
 
+def test_boxcar_nonfinite_handling(mk_test_img):
+    image = mk_test_img
+    image.data[14, 2] = np.nan
+    image.data[14, 4] = np.inf
+
+    trace = FlatTrace(image, 15.0)
+    boxcar = BoxcarExtract(image, trace, width=6, mask_treatment='filter')
+    spectrum = boxcar()
+    target = np.full_like(spectrum.flux.value, 90.)
+    target[2] = np.nan
+    target[4] = np.inf
+    np.testing.assert_equal(spectrum.flux.value, target)
+
+    trace = FlatTrace(image, 15.0)
+    boxcar = BoxcarExtract(image, trace, width=6, mask_treatment='exclude')
+    spectrum = boxcar()
+    target = np.full_like(spectrum.flux.value, 90.)
+    target[[2, 4]] = 91.2
+    np.testing.assert_allclose(spectrum.flux.value, target)
+
+
 def test_boxcar_outside_image_condition(mk_test_img):
     #
     # Trace is such that extraction aperture lays partially outside the image
     #
     image = mk_test_img
 
     trace = FlatTrace(image, 3.0)
-    boxcar = BoxcarExtract(image, trace)
+    boxcar = BoxcarExtract(image, trace, mask_treatment='filter')
 
     spectrum = boxcar(width=10.)
     assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 32.0))