Fix calculation of FWHM when the input data has non-finite values

CHANGES.rst

@@ -19,6 +19,7 @@ Other Changes and Additions
 Bug Fixes
 ^^^^^^^^^
 + Fixed dependence on non-release version of astrowidgets for overwrite capability on output images. [#108]
++ Fixed computation of FWHM when fitting to data that includes NaNs. [#164]
 
 1.3.9 (2023-06-16)
 ------------------

stellarphot/photometry/source_detection.py

@@ -36,10 +36,16 @@ def _fit_2dgaussian(data):
     gfit : `astropy.modeling.Model`
         The best-fit 2D Gaussian model.
     """
-    props = data_properties(data - np.min(data))
+    mask = ~np.isfinite(data)
+
+    # If there are non-finite pixels they need to be masked out
+    # or many of the returned properties will be NaN.
+    props = data_properties(data - np.min(data[~mask]), mask=mask)
 
     init_const = 0.  # subtracted data minimum above
-    init_amplitude = np.ptp(data)
+    # ptp = peak-to-peak, i.e. max - min, need to also exclude non-finite
+    # values here.
+    init_amplitude = np.ptp(data[~mask])
 
     g_init = (Const2D(init_const)
                   + Gaussian2D(amplitude=init_amplitude,
@@ -262,8 +268,8 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
         x, y = compute_fwhm(ccd, sources, fwhm_estimate=fwhm,
                             x_column='xcentroid', y_column='ycentroid',
                             sky_per_pix_avg=sky_per_pix_avg)
-        sources['fwhm_x'] = x 
-        sources['fwhm_y'] = y 
+        sources['fwhm_x'] = x
+        sources['fwhm_y'] = y
         sources['width'] = (x + y) / 2 # Average of x and y FWHM
 
         # Flag bogus fwhm values returned from fitting (no objects
@@ -297,5 +303,3 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5,
 
     sl_data = SourceListData(input_data=sources, colname_map=colnamemap)
     return sl_data
-
-

stellarphot/photometry/tests/test_detection.py

@@ -25,6 +25,28 @@ def test_compute_fwhm(units):
     fwhm_x, fwhm_y = compute_fwhm(fake_image.image, sources,
                                   x_column='x_mean', y_column='y_mean')
 
+    expected_fwhm = np.array(sources['x_stddev'] * gaussian_sigma_to_fwhm)
+    assert np.allclose(fwhm_x, expected_fwhm, rtol=1e-2)
+
+
+def test_compute_fwhm_with_NaNs():
+    # Regression test for https://github.com/feder-observatory/stellarphot/issues/161
+    # We should be able to find FWHM for a source even with NaNs in the image.
+    fake_image = FakeImage()
+    sources = fake_image.sources
+    x, y = sources['x_mean'].astype(int)[0], sources['y_mean'].astype(int)[0]
+    image = fake_image.image.copy()
+
+    # Add a NaN to the image at the location of the first source. Note the
+    # usual row/column swap when going to x/y coordinates.
+    image[y, x] = np.nan
+
+    fwhm_x, fwhm_y = compute_fwhm(image, sources,
+                                  x_column='x_mean', y_column='y_mean', fit=True)
+
+    expected_fwhm = np.array(sources['x_stddev'] * gaussian_sigma_to_fwhm)
+    assert np.allclose(fwhm_x, expected_fwhm, rtol=1e-2)
+
 
 def test_detect_source_number_location():
     """
@@ -76,4 +98,4 @@ def test_detect_source_with_padding():
                                      sky_per_pix_avg=sky_per_pix, padding=95)
 
     # Did we drop one source because it was too close to the edge?
-    assert len(sources) - 1 == len(found_sources)
+    assert len(sources) - 1 == len(found_sources)