fix input for background.traces, raise error in FlatTrace for negative trace

fix docstring

fix for masked values in FitTrace

.

.

.

..

,

.

.

Author: cshanahan1

CHANGES.rst

@@ -39,7 +39,12 @@ Bug Fixes
   peaks. Previously for Gaussian, the entire fit failed. [#205, #206]
 
 - Fixed input of `traces` in `Background`. Added a condition to 'FlatTrace' that
-  trace position must be a positive number. [#211]
+
+- Fix in FitTrace to set fully-masked column bin peaks to NaN. Previously, for
+  peak_method='max' these were set to 0.0, and for peak_method='centroid' they
+  were set to the number of rows in the image, biasing the final fit to all bin
+  peaks. [#206]
+
 
 Other changes
 ^^^^^^^^^^^^^

specreduce/background.py

@@ -32,9 +32,9 @@ class Background(_ImageParser):
     ----------
     image : `~astropy.nddata.NDData`-like or array-like
         image with 2-D spectral image data
-    traces : trace, int, float (single or list)
+    traces : List, `tracing.Trace`, int, float
         Individual or list of trace object(s) (or integers/floats to define
-        FlatTraces) to extract the background. If None, a FlatTrace at the
+        FlatTraces) to extract the background. If None, a `FlatTrace` at the
         center of the image (according to `disp_axis`) will be used.
     width : float
         width of extraction aperture in pixels
@@ -44,8 +44,22 @@ class Background(_ImageParser):
         pixels.
     disp_axis : int
         dispersion axis
+        [default: 1]
     crossdisp_axis : int
         cross-dispersion axis
+        [default: 0]
+    mask_treatment : string, optional
+        The method for handling masked or non-finite data. Choice of `filter`,
+        `omit`, or `zero-fill`. If `filter` is chosen, masked and non-finite
+        data will not contribute to the background statistic that is calculated
+        in each column along `disp_axis`. If `omit` is chosen, columns along
+        disp_axis with any masked/non-finite data values will be fully masked
+        (i.e, 2D mask is collapsed to 1D and applied). If `zero-fill` is chosen,
+        masked/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.
+        [default: ``filter``]
     """
     # required so numpy won't call __rsub__ on individual elements
     # https://stackoverflow.com/a/58409215
@@ -57,6 +71,8 @@ class Background(_ImageParser):
     statistic: str = 'average'
     disp_axis: int = 1
     crossdisp_axis: int = 0
+    mask_treatment : str = 'filter'
+    _valid_mask_treatment_methods = ('filter', 'omit', 'zero-fill')
 
     # TO-DO: update bkg_array with Spectrum1D alternative (is bkg_image enough?)
     bkg_array = deprecated_attribute('bkg_array', '1.3')
@@ -82,9 +98,29 @@ def __post_init__(self):
             dispersion axis
         crossdisp_axis : int
             cross-dispersion axis
+        mask_treatment : string
+            The method for handling masked or non-finite data. Choice of `filter`,
+            `omit`, or `zero-fill`. If `filter` is chosen, masked/non-finite data
+            will be filtered during the fit to each bin/column (along disp. axis) to
+            find the peak. If `omit` is chosen, columns along disp_axis with any
+            masked/non-finite data values will be fully masked (i.e, 2D mask is
+            collapsed to 1D and applied). If `zero-fill` is chosen, masked/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.
         """
+
+        # Parse image, including masked/nonfinite data handling based on
+        # choice of `mask_treatment`. Any uncaught nonfinte data values will be
+        # masked as well. Returns a Spectrum1D.
         self.image = self._parse_image(self.image)
 
+        # always work with masked array, even if there is no masked
+        # or nonfinite data, in case padding is needed. if not, mask will be
+        # dropped at the end and a regular array will be returned.
+        img = np.ma.masked_array(self.image.data, self.image.mask)
+
         if self.width < 0:
             raise ValueError("width must be positive")
         if self.width == 0:
@@ -95,9 +131,13 @@ def __post_init__(self):
 
         bkg_wimage = np.zeros_like(self.image.data, dtype=np.float64)
         for trace in self.traces:
+            # note: ArrayTrace can have masked values, but if it does a MaskedArray
+            # will be returned so this should be reflected in the window size here
+            # (i.e, np.nanmax is not required.)
             windows_max = trace.trace.data.max() + self.width/2
             windows_min = trace.trace.data.min() - self.width/2
-            if windows_max >= self.image.shape[self.crossdisp_axis]:
+
+            if windows_max > self.image.shape[self.crossdisp_axis]:
                 warnings.warn("background window extends beyond image boundaries " +
                               f"({windows_max} >= {self.image.shape[self.crossdisp_axis]})")
             if windows_min < 0:
@@ -115,27 +155,26 @@ def __post_init__(self):
             raise ValueError("background regions overlapped")
         if np.any(np.sum(bkg_wimage, axis=self.crossdisp_axis) == 0):
             raise ValueError("background window does not remain in bounds across entire dispersion axis")  # noqa
+        # check if image contained within background window is fully-nonfinite and raise an error
+        if np.all(img.mask[bkg_wimage > 0]):
+            raise ValueError("Image is fully masked within background window determined by `width`.")
 
         if self.statistic == 'median':
             # make it clear in the expose image that partial pixels are fully-weighted
             bkg_wimage[bkg_wimage > 0] = 1
 
         self.bkg_wimage = bkg_wimage
 
-        # mask user-highlighted and invalid values (if any) before taking stats
-        or_mask = (np.logical_or(~np.isfinite(self.image.data), self.image.mask)
-                   if self.image.mask is not None
-                   else ~np.isfinite(self.image.data))
-
         if self.statistic == 'average':
-            image_ma = np.ma.masked_array(self.image.data, mask=or_mask)
-            self._bkg_array = np.ma.average(image_ma,
+            self._bkg_array = np.ma.average(img,
                                             weights=self.bkg_wimage,
-                                            axis=self.crossdisp_axis).data
+                                            axis=self.crossdisp_axis)
+
         elif self.statistic == 'median':
-            med_mask = np.logical_or(self.bkg_wimage == 0, or_mask)
-            image_ma = np.ma.masked_array(self.image.data, mask=med_mask)
-            self._bkg_array = np.ma.median(image_ma, axis=self.crossdisp_axis).data
+            # combine where background weight image is 0 with image masked (which already
+            # accounts for non-finite data that wasn't already masked)
+            img.mask = np.logical_or(self.bkg_wimage == 0, self.image.mask)
+            self._bkg_array = np.ma.median(img, axis=self.crossdisp_axis)
         else:
             raise ValueError("statistic must be 'average' or 'median'")
 
@@ -204,7 +243,19 @@ def two_sided(cls, image, trace_object, separation, **kwargs):
             dispersion axis
         crossdisp_axis : int
             cross-dispersion axis
+        mask_treatment : string
+            The method for handling masked or non-finite data. Choice of `filter`,
+            `omit`, or `zero-fill`. If `filter` is chosen, masked/non-finite data
+            will be filtered during the fit to each bin/column (along disp. axis) to
+            find the peak. If `omit` is chosen, columns along disp_axis with any
+            masked/non-finite data values will be fully masked (i.e, 2D mask is
+            collapsed to 1D and applied). If `zero-fill` is chosen, masked/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.
         """
+
         image = _ImageParser._get_data_from_image(image) if image is not None else cls.image
         kwargs['traces'] = [trace_object-separation, trace_object+separation]
         return cls(image=image, **kwargs)
@@ -241,6 +292,17 @@ def one_sided(cls, image, trace_object, separation, **kwargs):
             dispersion axis
         crossdisp_axis : int
             cross-dispersion axis
+        mask_treatment : string
+            The method for handling masked or non-finite data. Choice of `filter`,
+            `omit`, or `zero-fill`. If `filter` is chosen, masked/non-finite data
+            will be filtered during the fit to each bin/column (along disp. axis) to
+            find the peak. If `omit` is chosen, columns along disp_axis with any
+            masked/non-finite data values will be fully masked (i.e, 2D mask is
+            collapsed to 1D and applied). If `zero-fill` is chosen, masked/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.
         """
         image = _ImageParser._get_data_from_image(image) if image is not None else cls.image
         kwargs['traces'] = [trace_object+separation]

specreduce/core.py

@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
+from copy import deepcopy
 import inspect
 from dataclasses import dataclass
 
@@ -53,14 +54,20 @@ def _parse_image(self, image, disp_axis=1):
 
         img = self._get_data_from_image(image)
 
-        # mask and uncertainty are set as None when they aren't specified upon
-        # creating a Spectrum1D object, so we must check whether these
-        # attributes are absent *and* whether they are present but set as None
-        if getattr(image, 'mask', None) is not None:
-            mask = image.mask
-        else:
-            mask = ~np.isfinite(img)
+        mask = getattr(image, 'mask', None)
+
+        # next, handle masked and nonfinite data in image.
+        # A mask will be created from any nonfinite image data, and combined
+        # with any additional 'mask' passed in. If image is being parsed within
+        # a specreduce operation that has 'mask_treatment' options, this will be
+        # handled as well. Note that input data may be modified if a fill value
+        # is chosen to handle masked data. The returned image will always have
+        # `image.mask` even if there are no nonfinte or masked values.
+        img, mask = self._mask_and_nonfinite_data_handling(image=img, mask=mask)
 
+        # mask (handled above) and uncertainty are set as None when they aren't
+        # specified upon creating a Spectrum1D object, so we must check whether
+        # these attributes are absent *and* whether they are present but set as None
         if getattr(image, 'uncertainty', None) is not None:
             uncertainty = image.uncertainty
         else:
@@ -71,8 +78,107 @@ def _parse_image(self, image, disp_axis=1):
         spectral_axis = getattr(image, 'spectral_axis',
                                 np.arange(img.shape[disp_axis]) * u.pix)
 
-        return Spectrum1D(img * unit, spectral_axis=spectral_axis,
-                          uncertainty=uncertainty, mask=mask)
+
+        img = Spectrum1D(img * unit, spectral_axis=spectral_axis,
+                         uncertainty=uncertainty, mask=mask)
+
+        return img
+
+    @staticmethod
+    def _get_data_from_image(image):
+        """Extract data array from various input types for `image`.
+           Retruns `np.ndarray` of image data."""
+
+        if isinstance(image, u.quantity.Quantity):
+            img = image.value
+        if isinstance(image, np.ndarray):
+            img = image
+        else:  # NDData, including CCDData and Spectrum1D
+            img = image.data
+        return img
+
+    def _mask_and_nonfinite_data_handling(self, image, mask):
+        """
+        This function handles the treatment of masked and nonfinite data,
+        including input validation.
+
+        All operations in Specreduce can take in a mask for the data as
+        part of the input NDData. Additionally, any non-finite values in the
+        data that aren't in the user-supplied mask will be combined bitwise
+        with the input mask.
+
+        There are three options currently implemented for the treatment
+        of masked and nonfinite data - filter, omit, and zero-fill.
+        Depending on the step, all or a subset of these three options are valid.
+
+        """
+
+        # valid options depend on Specreduce step, and are set accordingly there
+        # for steps that this isn't implemeted for yet, default to 'filter',
+        # which will return unmodified input image and mask
+        mask_treatment = getattr(self, 'mask_treatment', 'filter')
+
+        # make sure chosen option is valid. if _valid_mask_treatment_methods
+        # is not an attribue, proceed with 'filter' to return back inupt data
+        # and mask that is combined with nonfinite data.
+        if mask_treatment is not None:  # None in operations where masks aren't relevant (e.g FlatTrace)
+            valid_mask_treatment_methods = getattr(self, '_valid_mask_treatment_methods', ['filter'])
+            if mask_treatment not in valid_mask_treatment_methods:
+                raise ValueError(f"`mask_treatment` must be one of {valid_mask_treatment_methods}")
+
+        # make sure there is always a 'mask', even when all data is unmasked and finite.
+        if mask is not None:
+            mask = self.image.mask
+            # always mask any previously uncaught nonfinite values in image array
+            # combining these with the (optional) user-provided mask on `image.mask`
+            mask = np.logical_or(mask, ~np.isfinite(image))
+        else:
+            mask = ~np.isfinite(image)
+
+        # if mask option is the default 'filter' option, or None,
+        # nothing needs to be done. input mask (combined with nonfinite data)
+        # remains with data as-is.
+
+        if mask_treatment == 'zero-fill':
+            # make a copy of the input image since we will be modifying it
+            image = deepcopy(image)
+
+            # if mask_treatment is 'zero_fill', set masked values to zero in
+            # image data and drop image.mask. note that this is done after
+            # _combine_mask_with_nonfinite_from_data, so non-finite values in
+            # data (which are now in the mask) will also be set to zero.
+            # set masked values to zero
+            image[mask] = 0.
+
+            # masked array with no masked values, so accessing image.mask works
+            # but we don't need the actual mask anymore since data has been set to 0
+            mask = np.zeros(image.shape)
+
+        elif mask_treatment == 'omit':
+            # collapse 2d mask (after accounting for addl non-finite values in
+            # data) to a 1d mask, along dispersion axis, to fully mask columns
+            # that have any masked values.
+
+            # must have a crossdisp_axis specified to use 'omit' optoin
+            if hasattr(self, 'crossdisp_axis'):
+                crossdisp_axis = self.crossdisp_axis
+            if hasattr(self, '_crossdisp_axis'):
+                crossdisp_axis = self._crossdisp_axis
+
+            # create a 1d mask along crossdisp axis - if any column has a single nan,
+            # the entire column should be masked
+            reduced_mask = np.logical_or.reduce(mask,
+                                                axis=crossdisp_axis)
+
+            # back to a 2D mask
+            shape = (image.shape[0], 1) if crossdisp_axis == 0 else (1, image.shape[1])
+            mask = np.tile(reduced_mask, shape)
+
+        # check for case where entire image is masked.
+        if mask.all():
+            raise ValueError('Image is fully masked. Check for invalid values.')
+
+        return image, mask
 
     @staticmethod
     def _get_data_from_image(image):