set bin peaks to NaN in FitTrace when fully masked bin is encountered

CHANGES.rst

@@ -34,6 +34,11 @@ Bug Fixes
   for the Trace. Warning messages for ``peak_method`` == ``max`` and ``centroid``
   are also now reflective of what the bin peak is being set to. [#205]
 
+- 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/tests/test_tracing.py

@@ -183,7 +183,41 @@ def test_window_fit_trace(self):
 
     @pytest.mark.filterwarnings("ignore:The fit may be unsuccessful")
     @pytest.mark.filterwarnings("ignore:Model is linear in parameters")
-    def test_fit_trace_all_nan_columns(self):
+    @pytest.mark.filterwarnings("ignore:All pixels in bins")
+    def test_fit_trace_all_nan_cols(self):
+
+        # make sure that the actual trace that is fit is correct when
+        # all-masked bin peaks are set to NaN
+        img = self.mk_img(nrows=10, ncols=11)
+
+        img[:, 7] = np.nan
+        img[:, 4] = np.nan
+        img[:, 0] = np.nan
+
+        # peak_method = 'max'
+        truth = [1.6346154, 2.2371795, 2.8397436, 3.4423077, 4.0448718,
+                 4.6474359, 5.25, 5.8525641, 6.4551282, 7.0576923,
+                 7.6602564]
+        max_trace = FitTrace(img, peak_method='max')
+        np.testing.assert_allclose(truth, max_trace.trace)
+
+        # peak_method = 'gaussian'
+        truth = [1.947455, 2.383634, 2.8198131, 3.2559921, 3.6921712,
+                 4.1283502, 4.5645293, 5.0007083, 5.4368874, 5.8730665,
+                 6.3092455]
+        max_trace = FitTrace(img, peak_method='gaussian')
+        np.testing.assert_allclose(truth, max_trace.trace)
+
+        # peak_method = 'centroid'
+        truth = [2.5318835, 2.782069, 3.0322546, 3.2824402, 3.5326257,
+                 3.7828113, 4.0329969, 4.2831824, 4.533368, 4.7835536,
+                 5.0337391]
+        max_trace = FitTrace(img, peak_method='centroid')
+        np.testing.assert_allclose(truth, max_trace.trace)
+
+    @pytest.mark.filterwarnings("ignore:The fit may be unsuccessful")
+    @pytest.mark.filterwarnings("ignore:Model is linear in parameters")
+    def test_warn_msg_fit_trace_all_nan_cols(self):
 
         img = self.mk_img()
 
@@ -195,18 +229,15 @@ def test_fit_trace_all_nan_columns(self):
         mask[:, 30] = 1
         nddat = NDData(data=img, mask=mask, unit=u.DN)
 
-        match_str = 'All pixels in bins 20, 30, 100 are fully masked. '
+        match_str = 'All pixels in bins 20, 30, 100 are fully masked. Setting bin peaks to NaN.'
 
-        with pytest.warns(UserWarning, match=match_str +
-                          'Setting bin peaks to zero.'):
+        with pytest.warns(UserWarning, match=match_str):
             FitTrace(nddat, peak_method='max')
 
-        with pytest.warns(UserWarning, match=match_str +
-                          'Setting bin peaks to largest bin index \\(200\\)'):
+        with pytest.warns(UserWarning, match=match_str):
             FitTrace(nddat, peak_method='centroid')
 
-        with pytest.warns(UserWarning, match=match_str +
-                          'Setting bin peaks to nan.'):
+        with pytest.warns(UserWarning, match=match_str):
             FitTrace(nddat, peak_method='gaussian')
 
         # and when many bins are masked, that the message is consolidated
@@ -215,5 +246,5 @@ def test_fit_trace_all_nan_columns(self):
         nddat = NDData(data=img, mask=mask, unit=u.DN)
         with pytest.warns(UserWarning, match='All pixels in bins '
                           '0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ..., 20 are '
-                          'fully masked. Setting bin peaks to zero.'):
+                          'fully masked. Setting bin peaks to NaN.'):
             FitTrace(nddat)

specreduce/tracing.py

@@ -317,14 +317,14 @@ def _fit_trace(self, img):
             # or just a single, unbinned column if no bins
             z_i = img[ilum2, x_bins[i]:x_bins[i+1]].sum(axis=self._disp_axis)
 
-            if self.peak_method == 'gaussian':
-                # if binned column is fully masked for peak_method='gaussian',
-                # the fit value for this bin should be nan, then continue to next
+            # if this bin is fully masked, set bin peak to NaN so it can be
+            # filtered in the final fit to all bin peaks for the trace
+            if z_i.mask.all():
+                warn_bins.append(i)
+                y_bins[i] = np.nan
+                continue
 
-                if z_i.mask.all():
-                    warn_bins.append(i)
-                    y_bins[i] = np.nan
-                    continue
+            if self.peak_method == 'gaussian':
 
                 peak_y_i = ilum2[z_i.argmax()]
 
@@ -353,10 +353,7 @@ def _fit_trace(self, img):
                     y_bins[i] = popt_i.mean_0.value
                     popt_tot = popt_i
 
-            if z_i.mask.all():  # all-masked bins when peak_method is 'centroid' or 'max'
-                warn_bins.append(i)
-
-            if self.peak_method == 'centroid':
+            elif self.peak_method == 'centroid':
                 z_i_cumsum = np.cumsum(z_i)
                 # find the interpolated index where the cumulative array reaches
                 # half the total cumulative values
@@ -380,14 +377,10 @@ def _fit_trace(self, img):
             if len(warn_bins) > 20:
                 warn_bins = warn_bins[0: 10] + ['...'] + [warn_bins[-1]]
 
-            # warning message printed out depends on `peak_method`
-            warn_msgs = {'gaussian': 'nan', 'max': 'zero',
-                         'centroid': f'largest bin index ({str(img.shape[0])})'}
-
             warnings.warn(f"All pixels in {'bins' if len(warn_bins) else 'bin'} "
                           f"{', '.join([str(x) for x in warn_bins])}"
                           " are fully masked. Setting bin"
-                          f" peak{'s' if len(warn_bins) else ''} to {warn_msgs[self.peak_method]}.")
+                          f" peak{'s' if len(warn_bins) else ''} to NaN.")
 
         # recenter bin positions
         x_bins = (x_bins[:-1] + x_bins[1:]) / 2