From f65a3b7c6fc2660ad30d6ddd99fd257064f8bfee Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Sun, 6 Sep 2026 20:27:02 +0100 Subject: [PATCH 1/2] - Fixed `find_arc_lines` to accept any of the Astropy's uncertainty types. --- CHANGES.rst | 10 +++++ specreduce/line_matching.py | 16 +++++++- specreduce/tests/test_line_matching.py | 54 +++++++++++++++++++++++++- 3 files changed, 77 insertions(+), 3 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 0464fbdc..8bc746a6 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -38,6 +38,16 @@ 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] + Other changes ^^^^^^^^^^^^^ diff --git a/specreduce/line_matching.py b/specreduce/line_matching.py index f8f1b807..423f7497 100644 --- a/specreduce/line_matching.py +++ b/specreduce/line_matching.py @@ -29,8 +29,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, the square root of + the absolute flux is used as an estimate. fwhm Estimated full-width half-maximum of the lines in pixels. @@ -55,9 +60,16 @@ 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 it as the square root of the flux. + # 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))) + 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"] diff --git a/specreduce/tests/test_line_matching.py b/specreduce/tests/test_line_matching.py index a96d0675..9c229fcc 100644 --- a/specreduce/tests/test_line_matching.py +++ b/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,58 @@ 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) + + +@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.remote_data @pytest.mark.filterwarnings("ignore:No observer defined on WCS") @pytest.mark.filterwarnings("ignore:Model is linear in parameters") From c9bc1b919879ca6bedffc8ec1060170668668d26 Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Sun, 6 Sep 2026 21:24:38 +0100 Subject: [PATCH 2/2] =?UTF-8?q?-=20Added=20a=20=C2=B4specreduce.utils.util?= =?UTF-8?q?s.measure=5Fnoise=C2=B4=20function=20to=20robustly=20estimate?= =?UTF-8?q?=20noise=20in=20a=202D=20spectrum.=20This=20is=20for=20now=20us?= =?UTF-8?q?ed=20only=20in=20`specreduce.line=5Fmatching.find=5Farc=5Flines?= =?UTF-8?q?`=20as=20a=20fallback=20noise=20estimate=20if=20the=20spectrum?= =?UTF-8?q?=20doesn't=20contain=20uncertainties,=20but=20may=20be=20useful?= =?UTF-8?q?=20elsewhere=20as=20well.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGES.rst | 13 ++++ specreduce/line_matching.py | 16 +++-- specreduce/tests/test_line_matching.py | 21 ++++++- specreduce/tests/test_utils.py | 56 ++++++++++++++++- specreduce/utils/utils.py | 83 +++++++++++++++++++++++++- 5 files changed, 179 insertions(+), 10 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 8bc746a6..3370f55f 100644 --- a/CHANGES.rst +++ b/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 ^^^^^^^^^^^ @@ -48,6 +54,13 @@ Bug Fixes 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 ^^^^^^^^^^^^^ diff --git a/specreduce/line_matching.py b/specreduce/line_matching.py index 423f7497..e774a60e 100644 --- a/specreduce/line_matching.py +++ b/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"] @@ -34,8 +36,8 @@ def find_arc_lines( 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, the square root of - the absolute flux is used as an estimate. + 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. @@ -61,12 +63,14 @@ def find_arc_lines( 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 it as the square root of the flux. - # 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 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) diff --git a/specreduce/tests/test_line_matching.py b/specreduce/tests/test_line_matching.py index 9c229fcc..37e64c93 100644 --- a/specreduce/tests/test_line_matching.py +++ b/specreduce/tests/test_line_matching.py @@ -114,7 +114,7 @@ def synthetic_arc(): 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) + return x * u.pix, flux * u.DN, np.full(x.size, sigma), np.array(centers) @pytest.mark.filterwarnings("ignore:The fit may be unsuccessful") @@ -132,7 +132,7 @@ def test_find_arc_lines_uncertainty_types(synthetic_arc, uncertainty_cls, transf 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 + spectral_axis, flux, sigma, _ = synthetic_arc reference = Spectrum( flux=flux, spectral_axis=spectral_axis, uncertainty=StdDevUncertainty(sigma) ) @@ -152,6 +152,23 @@ def test_find_arc_lines_uncertainty_types(synthetic_arc, uncertainty_cls, transf 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") diff --git a/specreduce/tests/test_utils.py b/specreduce/tests/test_utils.py index 43408603..7cbbed6e 100644 --- a/specreduce/tests/test_utils.py +++ b/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)) diff --git a/specreduce/utils/utils.py b/specreduce/utils/utils.py index 3886f90a..513074e0 100644 --- a/specreduce/utils/utils.py +++ b/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,