find_arc_lines uncertainty handling fix

CHANGES.rst

@@ -30,6 +30,12 @@ New Features
   (``AWAV-GRA``) and vacuum (``WAVE-GRI``) spectral axis types in FITS WCS
   export. [#316]
 
+- Added ``specreduce.utils.measure_noise``, a robust estimator of the per-pixel
+  noise standard deviation of a 1D spectrum or of each row of a 2D spectral
+  image. It measures the sigma-clipped ``mad_std`` of the second difference of
+  the flux along the dispersion axis, so it is insensitive to a smooth continuum
+  or residual background and to emission lines. [#XXX]
+
 API Changes
 ^^^^^^^^^^^
 
@@ -38,6 +44,23 @@ API Changes
   on instantiation, and its removal has been rescheduled from v2.0 to v1.11.
   Use ``specreduce.wavecal1d.WavelengthCalibration1D`` instead. [#316]
 
+Bug Fixes
+^^^^^^^^^
+
+- ``line_matching.find_arc_lines`` now accepts spectra with any of the Astropy
+  uncertainty types (``StdDevUncertainty``, ``VarianceUncertainty``, or
+  ``InverseVariance``). Non-standard-deviation uncertainties are converted to
+  ``StdDevUncertainty`` on a copy of the spectrum before the line finding, which
+  previously failed with a unit conversion error for variance-type uncertainties.
+  [#XXX]
+
+- ``line_matching.find_arc_lines`` now estimates the noise of a spectrum without
+  an uncertainty using ``measure_noise`` instead of the square root of the absolute
+  flux. The old fallback amounted to a fixed detection threshold of
+  ``noise_factor**2`` flux units regardless of the actual noise, which missed
+  all lines in faint spectra and reported noise spikes as lines in noisy or
+  background-subtracted ones. [#XXX]
+
 Other changes
 ^^^^^^^^^^^^^
 

specreduce/line_matching.py

@@ -14,6 +14,8 @@
 
 from specutils import Spectrum
 
+from specreduce.utils.utils import measure_noise
+
 __all__ = ["find_arc_lines", "match_lines_wcs"]
 
 
@@ -29,8 +31,13 @@ def find_arc_lines(
 
     Parameters
     ----------
-    spectrum : The extracted arc spectrum to search for lines. It should be background-subtracted
-        and must have an "uncertainty" attribute.
+    spectrum
+        The extracted arc spectrum to search for lines. It should be background-subtracted.
+        The uncertainty can be any of the Astropy uncertainty types
+        (`~astropy.nddata.StdDevUncertainty`, `~astropy.nddata.VarianceUncertainty`, or
+        `~astropy.nddata.InverseVariance`); it is converted to a standard deviation
+        before the line finding. If the spectrum has no uncertainty, a constant per-pixel
+        noise is estimated from the data with `~specreduce.utils.utils.measure_noise`.
 
     fwhm
         Estimated full-width half-maximum of the lines in pixels.
@@ -55,9 +62,18 @@ def find_arc_lines(
     if fwhm.unit != spectrum.spectral_axis.unit:
         raise ValueError("fwhm must have the same units as spectrum.spectral_axis.")
 
+    # The line finding and fitting are always done using standard deviation uncertainties.
+    # If the spectrum has no uncertainty, estimate a constant per-pixel noise from the
+    # scatter in the data itself. If it has a variance or inverse variance uncertainty,
+    # convert it to a standard deviation. Either way, work on a copy so that the input
+    # spectrum is left untouched.
     if spectrum.uncertainty is None:
         spectrum = deepcopy(spectrum)
-        spectrum.uncertainty = StdDevUncertainty(np.sqrt(np.abs(spectrum.flux.value)))
+        noise = measure_noise(spectrum.flux.value)
+        spectrum.uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, noise))
+    elif not isinstance(spectrum.uncertainty, StdDevUncertainty):
+        spectrum = deepcopy(spectrum)
+        spectrum.uncertainty = spectrum.uncertainty.represent_as(StdDevUncertainty)
 
     detected_lines = find_lines_threshold(spectrum, noise_factor=noise_factor)
     detected_lines = detected_lines[detected_lines["line_type"] == "emission"]

specreduce/tests/test_line_matching.py

@@ -4,7 +4,7 @@
 
 from astropy.wcs import WCS
 from astropy.modeling import models
-from astropy.nddata import StdDevUncertainty
+from astropy.nddata import StdDevUncertainty, VarianceUncertainty, InverseVariance
 from specutils.fitting import fit_generic_continuum
 
 from specreduce.calibration_data import load_pypeit_calibration_lines
@@ -100,6 +100,75 @@ def test_find_arc_lines(mk_test_data):
     assert len(lines) > 1
 
 
+@pytest.fixture
+def synthetic_arc():
+    """
+    A small synthetic arc spectrum with Gaussian emission lines at known pixel positions
+    and a constant standard deviation of 3 DN per pixel.
+    """
+    rng = np.random.default_rng(42)
+    x = np.arange(500.0)
+    centers = [50.0, 120.0, 210.0, 330.0, 410.0]
+    flux = np.zeros_like(x)
+    for c, a in zip(centers, [200.0, 80.0, 500.0, 60.0, 150.0]):
+        flux += a * np.exp(-0.5 * ((x - c) / 2.1) ** 2)
+    sigma = 3.0
+    flux += rng.normal(0.0, sigma, x.size)
+    return x * u.pix, flux * u.DN, np.full(x.size, sigma), np.array(centers)
+
+
+@pytest.mark.filterwarnings("ignore:The fit may be unsuccessful")
+@pytest.mark.filterwarnings("ignore:Spectrum is not below the threshold")
+@pytest.mark.parametrize(
+    "uncertainty_cls, transform",
+    [
+        (StdDevUncertainty, lambda s: s),
+        (VarianceUncertainty, lambda s: s**2),
+        (InverseVariance, lambda s: 1.0 / s**2),
+    ],
+)
+def test_find_arc_lines_uncertainty_types(synthetic_arc, uncertainty_cls, transform):
+    """
+    find_arc_lines must accept any of the three astropy uncertainty types and produce
+    the same lines as it does for an equivalent StdDevUncertainty.
+    """
+    spectral_axis, flux, sigma, _ = synthetic_arc
+    reference = Spectrum(
+        flux=flux, spectral_axis=spectral_axis, uncertainty=StdDevUncertainty(sigma)
+    )
+    spectrum = Spectrum(
+        flux=flux, spectral_axis=spectral_axis, uncertainty=uncertainty_cls(transform(sigma))
+    )
+
+    expected = find_arc_lines(reference, fwhm=5, window=3, noise_factor=5)
+    lines = find_arc_lines(spectrum, fwhm=5, window=3, noise_factor=5)
+
+    assert len(expected) >= 5
+    assert len(lines) == len(expected)
+    np.testing.assert_allclose(lines["centroid"].value, expected["centroid"].value)
+    np.testing.assert_allclose(lines["fwhm"].value, expected["fwhm"].value)
+    np.testing.assert_allclose(lines["amplitude"].value, expected["amplitude"].value)
+    # The input spectrum must not be modified in place.
+    assert isinstance(spectrum.uncertainty, uncertainty_cls)
+
+
+@pytest.mark.filterwarnings("ignore:The fit may be unsuccessful")
+@pytest.mark.filterwarnings("ignore:Spectrum is not below the threshold")
+@pytest.mark.parametrize("scale", [1.0 / 30.0, 10.0 / 3.0], ids=["faint", "noisy"])
+def test_find_arc_lines_estimates_noise_without_uncertainty(synthetic_arc, scale):
+    """
+    Without an uncertainty, find_arc_lines must estimate the noise from the data so that
+    the detection threshold follows the actual noise level, independent of the flux scale.
+    """
+    spectral_axis, flux, _, centers = synthetic_arc
+    spectrum = Spectrum(flux=flux * scale, spectral_axis=spectral_axis)
+
+    lines = find_arc_lines(spectrum, fwhm=5, window=3, noise_factor=5)
+
+    assert len(lines) == len(centers)
+    np.testing.assert_allclose(np.sort(lines["centroid"].value), centers, atol=0.5)
+
+
 @pytest.mark.remote_data
 @pytest.mark.filterwarnings("ignore:No observer defined on WCS")
 @pytest.mark.filterwarnings("ignore:Model is linear in parameters")

specreduce/tests/test_utils.py

@@ -6,7 +6,7 @@
 
 from specutils import Spectrum
 from specreduce.tracing import FitTrace
-from specreduce.utils.utils import measure_cross_dispersion_profile
+from specreduce.utils.utils import measure_cross_dispersion_profile, measure_noise
 
 
 def mk_gaussian_img(nrows=20, ncols=16, mean=10, stddev=4):
@@ -169,3 +169,57 @@ def test_errors_warnings(self):
                                              'or None to use all '
                                              'cross-dispersion pixels.'):
             measure_cross_dispersion_profile(img, width='.')
+
+
+class TestMeasureNoise:
+    """Tests for the robust second-difference noise estimator."""
+
+    @staticmethod
+    def _arc(rng, x, centers, amplitudes, fwhm=4.0):
+        sigma = fwhm / 2.3548
+        flux = np.zeros_like(x)
+        for c, a in zip(centers, amplitudes):
+            flux += a * np.exp(-0.5 * ((x - c) / sigma) ** 2)
+        return flux
+
+    def test_white_noise(self):
+        rng = np.random.default_rng(1)
+        flux = rng.normal(0.0, 3.0, 2000)
+        assert measure_noise(flux) == pytest.approx(3.0, rel=0.05)
+
+    def test_lines_and_residual_background(self):
+        rng = np.random.default_rng(2)
+        x = np.arange(2000.0)
+        flux = self._arc(rng, x, rng.uniform(0, 2000, 40), 10 ** rng.uniform(1, 3.3, 40))
+        flux += 120.0 * (x / 2000.0 - 0.5) ** 2
+        flux += rng.normal(0.0, 3.0, x.size)
+        noise = measure_noise(flux)
+        assert 0.9 * 3.0 < noise < 1.3 * 3.0
+
+    def test_2d_per_row_along_dispersion_axis(self):
+        rng = np.random.default_rng(3)
+        sigmas = np.array([1.0, 2.0, 4.0, 8.0])
+        flux = rng.normal(0.0, sigmas[:, None], (4, 1000))
+        noise = measure_noise(flux, axis=-1)
+        assert noise.shape == (4,)
+        np.testing.assert_allclose(noise, sigmas, rtol=0.1)
+        noise_t = measure_noise(flux.T, axis=0)
+        np.testing.assert_allclose(noise_t, noise)
+
+    def test_quantity_input_keeps_unit(self):
+        rng = np.random.default_rng(4)
+        flux = rng.normal(0.0, 3.0, 2000) * u.DN
+        noise = measure_noise(flux)
+        assert isinstance(noise, u.Quantity)
+        assert noise.unit == u.DN
+        assert noise.value == pytest.approx(3.0, rel=0.05)
+
+    def test_nan_pixels_are_ignored(self):
+        rng = np.random.default_rng(5)
+        flux = rng.normal(0.0, 3.0, 2000)
+        flux[100:110] = np.nan
+        assert measure_noise(flux) == pytest.approx(3.0, rel=0.05)
+
+    def test_too_short_input_raises(self):
+        with pytest.raises(ValueError, match="at least 5"):
+            measure_noise(np.zeros(4))

specreduce/utils/utils.py

@@ -1,10 +1,91 @@
 import numpy as np
+from astropy import units as u
+from astropy.stats import mad_std, sigma_clipped_stats
 
 from specreduce.core import parse_image
 from specreduce.tracing import Trace, FlatTrace
 from specreduce.extract import _ap_weight_image, _align_along_trace
 
-__all__ = ['measure_cross_dispersion_profile', '_align_along_trace']
+__all__ = ['measure_cross_dispersion_profile', 'measure_noise', '_align_along_trace']
+
+
+def measure_noise(
+    data: np.ndarray | u.Quantity,
+    axis: int = -1,
+    sigma: float = 3.0,
+    maxiters: int | None = 10,
+) -> float | np.ndarray | u.Quantity:
+    """
+    Estimate the per-pixel noise standard deviation of a spectrum from the data itself.
+
+    The estimate is the sigma-clipped median absolute deviation of the second difference
+    of the flux along the dispersion axis, ``2 f[i] - f[i-2] - f[i+2]``, scaled to the
+    standard deviation of a single pixel. Differencing removes any smooth continuum or
+    residual background, so the estimate does not depend on the spectrum being
+    background-subtracted, while the iterative sigma clipping removes the pixels
+    dominated by emission or absorption lines before the scatter is measured.
+
+    The second difference of white noise has a variance of six times the per-pixel
+    variance, so the clipped ``mad_std`` of the differences is divided by the square
+    root of six. This is the same differencing used by the DER_SNR algorithm
+    (Stoehr et al. 2008), with the plain median replaced by a sigma-clipped robust
+    standard deviation to reduce the bias from dense line lists.
+
+    Parameters
+    ----------
+    data
+        The flux array. Can be 1D or N-dimensional; for a 2D spectral image the noise
+        is estimated separately along ``axis`` for each row (or column). Non-finite
+        values and masked elements of a masked array are ignored.
+    axis
+        The dispersion axis along which the differences are taken.
+    sigma
+        The clipping threshold in units of the robust standard deviation.
+    maxiters
+        The maximum number of clipping iterations, or `None` to iterate until
+        convergence.
+
+    Returns
+    -------
+    float, ndarray, or Quantity
+        The estimated noise standard deviation. A scalar for 1D input, otherwise an
+        array with ``axis`` removed. If ``data`` is a `~astropy.units.Quantity`, the
+        result carries the same unit.
+
+    Notes
+    -----
+    The estimator assumes that the noise is uncorrelated between pixels two apart.
+    For spectra that have been smoothed or resampled onto a finer grid, the
+    differencing suppresses part of the correlated noise and the result is biased
+    low by up to a few tens of percent.
+    """
+    unit = None
+    if isinstance(data, u.Quantity):
+        unit = data.unit
+        data = data.value
+    if np.ma.isMaskedArray(data):
+        data = data.astype(float).filled(np.nan)
+    values = np.asarray(data, dtype=float)
+
+    if values.ndim == 0 or values.shape[axis] < 5:
+        raise ValueError("measure_noise requires at least 5 pixels along the dispersion axis.")
+
+    values = np.moveaxis(values, axis, -1)
+    diff2 = 2.0 * values[..., 2:-2] - values[..., :-4] - values[..., 4:]
+
+    # Non-finite differences (from NaN or masked pixels) are excluded through an explicit
+    # mask. The placeholder value is never used, but must be finite to keep astropy from
+    # warning about invalid input.
+    invalid = ~np.isfinite(diff2)
+    diff2 = np.where(invalid, 0.0, diff2)
+    _, _, clipped_std = sigma_clipped_stats(
+        diff2, mask=invalid, sigma=sigma, maxiters=maxiters, stdfunc=mad_std, axis=-1
+    )
+    noise = clipped_std / np.sqrt(6.0)
+
+    if unit is not None:
+        noise = noise * unit
+    return noise
 
 
 def measure_cross_dispersion_profile(image, trace=None, crossdisp_axis=0,