From 3aa5acb04ac6fd49f303572f5f4368e894ec1ddc Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Thu, 3 Sep 2026 20:20:47 +0100 Subject: [PATCH 1/4] Fixed uncertainty propagation in WavelengthSolution1D.resample(). --- CHANGES.rst | 7 +++++++ specreduce/tests/test_wavesol1d.py | 30 +++++++++++++++++++++++++++++- specreduce/wavesol1d.py | 4 ++-- 3 files changed, 38 insertions(+), 3 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 0464fbdc..cf5e3916 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -30,6 +30,13 @@ New Features (``AWAV-GRA``) and vacuum (``WAVE-GRI``) spectral axis types in FITS WCS export. [#316] +Bug Fixes +^^^^^^^^^ + +- Fixed the variance propagation in ``WavelengthSolution1D.resample()``, which + weighted each pixel's variance by its wavelength width instead of the square of the + width, underestimating the uncertainty of the resampled flux. [#XXX] + API Changes ^^^^^^^^^^^ diff --git a/specreduce/tests/test_wavesol1d.py b/specreduce/tests/test_wavesol1d.py index 358ac33c..93f947d9 100644 --- a/specreduce/tests/test_wavesol1d.py +++ b/specreduce/tests/test_wavesol1d.py @@ -5,7 +5,7 @@ from astropy.io import fits from astropy.modeling import models, fitting from astropy.modeling.polynomial import Polynomial1D -from astropy.nddata import StdDevUncertainty +from astropy.nddata import StdDevUncertainty, VarianceUncertainty from gwcs import coordinate_frames, wcs from astropy.utils.exceptions import AstropyUserWarning from astropy.wcs import WCS as astropy_WCS @@ -419,3 +419,31 @@ def test_from_asdf_raises_on_unsupported_transform(tmp_path): asdf.AsdfFile({"wavelength_solution": {"gwcs": w, "wave_air": False}}).write_to(path) with pytest.raises(ValueError, match="shift followed by a polynomial"): WavelengthSolution1D.from_asdf(path) + + +def test_resample_propagates_variance(mk_ws_with_transform): + ws = mk_ws_with_transform + npix = pix_bounds[1] + sigma = 3.0 + spectrum = Spectrum( + flux=np.ones(npix) * u.count, + spectral_axis=np.arange(npix) * u.pix, + uncertainty=StdDevUncertainty(np.full(npix, sigma)), + ) + result = ws.resample(spectrum, nbins=10) + var = result.uncertainty.represent_as(VarianceUncertainty).array + dl = np.diff(result.spectral_axis.value)[0] + dldx = np.diff(ws.p2w(np.arange(npix + 1) - 0.5)) + + # The output is sum(w * f * dldx) / dl, so the variance of an interior bin is + # sum(w**2 * sigma**2 * dldx**2) / dl**2 over the contributing pixels. + ibin = 5 + edges = result.spectral_axis.value[ibin] + np.array([-0.5, 0.5]) * dl + x_edges = ws.wav_to_pix(edges) + 0.5 + i1, i2 = np.floor(x_edges).astype(int) + w = np.zeros(npix) + w[i1 + 1 : i2] = 1.0 + w[i1] = 1.0 - (x_edges[0] - i1) + w[i2] = x_edges[1] - i2 + expected = (w**2 * sigma**2 * dldx**2).sum() / dl**2 + np.testing.assert_allclose(var[ibin], expected, rtol=1e-6) diff --git a/specreduce/wavesol1d.py b/specreduce/wavesol1d.py index dd680165..51e03b4c 100644 --- a/specreduce/wavesol1d.py +++ b/specreduce/wavesol1d.py @@ -712,11 +712,11 @@ def resample( sl = slice(i1, i2 + 1) w = weights[sl] flux_wl[i] = (w * flux[sl] * dldx[sl]).sum() - ucty_wl[i] = (w**2 * ucty[sl] * dldx[sl]).sum() + ucty_wl[i] = (w**2 * ucty[sl] * dldx[sl] ** 2).sum() else: fracw = bin_edges_pix[i + 1] - bin_edges_pix[i] flux_wl[i] = fracw * flux[i1] * dldx[i1] - ucty_wl[i] = fracw**2 * ucty[i1] * dldx[i1] + ucty_wl[i] = fracw**2 * ucty[i1] * dldx[i1] ** 2 bin_widths_wav = np.diff(bin_edges_wav) flux_wl = flux_wl / bin_widths_wav * spectrum.flux.unit / self.unit From 72c12c17cdfea6a04afe18359edb8fc1b0460005 Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Sat, 5 Sep 2026 14:06:23 +0100 Subject: [PATCH 2/4] - Fixed a bug in WavelengthSolution1D.resample, where the method was multiplying the per-pixel flux by the pixel wavelength width before dividing by the bin width. The output was labelled as a flux density but was numerically off by the local dispersion and did not conserve the integrated flux. I really thought I fixed this months ago already... --- CHANGES.rst | 10 +++++-- specreduce/tests/test_wavesol1d.py | 47 +++++++++++++++++++++++++----- specreduce/wavesol1d.py | 23 ++++++++------- 3 files changed, 60 insertions(+), 20 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index cf5e3916..322b2d41 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -33,9 +33,13 @@ New Features Bug Fixes ^^^^^^^^^ -- Fixed the variance propagation in ``WavelengthSolution1D.resample()``, which - weighted each pixel's variance by its wavelength width instead of the square of the - width, underestimating the uncertainty of the resampled flux. [#XXX] +- Fixed ``WavelengthSolution1D.resample()``, which multiplied the per-pixel flux by + the wavelength width of each pixel before dividing by the bin width. The output was + labelled as a flux density per wavelength unit but was numerically off by the local + dispersion, so a flat spectrum in counts per pixel stayed flat instead of falling where + the dispersion grows, and the integrated flux was not conserved. The flux density and + its uncertainty are now computed from the fractional pixel overlaps alone, with the + variance carrying the squared weights. [#XXX] API Changes ^^^^^^^^^^^ diff --git a/specreduce/tests/test_wavesol1d.py b/specreduce/tests/test_wavesol1d.py index 93f947d9..377ca9c2 100644 --- a/specreduce/tests/test_wavesol1d.py +++ b/specreduce/tests/test_wavesol1d.py @@ -114,10 +114,12 @@ def test_resample(mk_spectrum, mk_ws_with_transform, mk_ws_without_transform): assert len(resampled.flux) == 50 assert resampled.flux.unit == u.count / u.angstrom - pix_edges = np.arange(spectrum.spectral_axis.size + 1) - 0.5 - f0 = (spectrum.flux.value * np.diff(ws._p2w(pix_edges))).sum() + # the input is in counts per pixel, so integrating the output density over the + # wavelength bins must return the total input counts, up to the accuracy of the + # interpolated wavelength-to-pixel inverse at the outermost bin edges + f0 = spectrum.flux.value.sum() f1 = (resampled.flux.value * np.diff(resampled.spectral_axis.value)[0]).sum() - np.testing.assert_approx_equal(f0, f1, 5) + np.testing.assert_allclose(f0, f1, rtol=1e-4) resampled = ws.resample(spectrum, wlbounds=wav_bounds) resampled = ws.resample(spectrum, bin_edges=np.linspace(*wav_bounds, num=50)) @@ -433,10 +435,9 @@ def test_resample_propagates_variance(mk_ws_with_transform): result = ws.resample(spectrum, nbins=10) var = result.uncertainty.represent_as(VarianceUncertainty).array dl = np.diff(result.spectral_axis.value)[0] - dldx = np.diff(ws.p2w(np.arange(npix + 1) - 0.5)) - # The output is sum(w * f * dldx) / dl, so the variance of an interior bin is - # sum(w**2 * sigma**2 * dldx**2) / dl**2 over the contributing pixels. + # The output is sum(w * f) / dl, so the variance of an interior bin is + # sum(w**2 * sigma**2) / dl**2 over the contributing pixels. ibin = 5 edges = result.spectral_axis.value[ibin] + np.array([-0.5, 0.5]) * dl x_edges = ws.wav_to_pix(edges) + 0.5 @@ -445,5 +446,37 @@ def test_resample_propagates_variance(mk_ws_with_transform): w[i1 + 1 : i2] = 1.0 w[i1] = 1.0 - (x_edges[0] - i1) w[i2] = x_edges[1] - i2 - expected = (w**2 * sigma**2 * dldx**2).sum() / dl**2 + expected = (w**2 * sigma**2).sum() / dl**2 np.testing.assert_allclose(var[ibin], expected, rtol=1e-6) + + +def test_resample_returns_flux_density_per_wavelength(): + """ + A flat spectrum in counts per pixel becomes counts per wavelength unit divided by + the local dispersion, and the integral over the bins conserves the total counts. + """ + npix, flux, sigma = 100, 100.0, 4.0 + spectrum = Spectrum( + flux=np.full(npix, flux) * u.count, + spectral_axis=np.arange(npix) * u.pix, + uncertainty=StdDevUncertainty(np.full(npix, sigma)), + ) + + # linear dispersion: with nbins == npix every bin is exactly one pixel wide + dispersion = 2.4 + ws = WavelengthSolution1D( + models.Shift(0) | models.Polynomial1D(1, c0=5000, c1=dispersion), (0, npix), u.angstrom + ) + result = ws.resample(spectrum) + assert result.flux.unit == u.count / u.angstrom + np.testing.assert_allclose(result.flux.value, flux / dispersion, rtol=1e-9) + np.testing.assert_allclose(result.uncertainty.array, sigma / dispersion, rtol=1e-9) + + # non-linear dispersion: the density follows 1 / dldx and the counts are conserved + p2w = models.Shift(0) | models.Polynomial1D(2, c0=5000, c1=1.0, c2=0.02) + ws = WavelengthSolution1D(p2w, (0, npix), u.angstrom) + result = ws.resample(spectrum, nbins=npix) + bin_widths = np.diff(np.linspace(*p2w(np.array([-0.5, npix - 0.5])), npix + 1)) + np.testing.assert_allclose((result.flux.value * bin_widths).sum(), flux * npix, rtol=1e-4) + dldx_at_centers = _diff_poly1d(p2w[1])(ws.wav_to_pix(result.spectral_axis.value)) + np.testing.assert_allclose(result.flux.value, flux / dldx_at_centers, rtol=2e-2) diff --git a/specreduce/wavesol1d.py b/specreduce/wavesol1d.py index 51e03b4c..b15d9f99 100644 --- a/specreduce/wavesol1d.py +++ b/specreduce/wavesol1d.py @@ -636,9 +636,12 @@ def resample( ) -> Spectrum: """Bin the given pixel-space 1D spectrum to a wavelength space conserving the flux. - This method bins a pixel-space spectrum to a wavelength space using the computed - pixel-to-wavelength and wavelength-to-pixel transformations and their derivatives with - respect to the spectral axis. The binning is exact and conserves the integrated flux. + The input flux is taken to be integrated per pixel (e.g. counts per pixel, as + produced by the extraction methods). Each wavelength bin receives the flux of the + pixels it overlaps, weighted by the overlapping fraction of each pixel, and the total + is divided by the bin width, so the output is a flux density per wavelength unit. + The binning is exact and conserves the integrated flux. The variance is propagated + with the squared weights, assuming independent pixel noise. Parameters ---------- @@ -662,7 +665,9 @@ def resample( Returns ------- - 1D spectrum binned to the specified wavelength bins. + 1D spectrum binned to the specified wavelength bins, with the flux in units of + the input flux unit per wavelength unit and the uncertainty in the same + uncertainty class as the input. """ if nbins is not None and nbins < 0: raise ValueError("Number of bins must be non-zero and positive.") @@ -700,8 +705,6 @@ def resample( ucty_wl = np.zeros(nbins) weights = np.zeros(npix) - dldx = np.diff(self.p2w(np.arange(pixels[0], pixels[-1] + 2) - 0.5)) - for i in range(nbins): i1, i2 = bin_edge_ix[i : i + 2] weights[:] = 0.0 @@ -711,12 +714,12 @@ def resample( weights[i2] = bin_edge_w[i + 1] sl = slice(i1, i2 + 1) w = weights[sl] - flux_wl[i] = (w * flux[sl] * dldx[sl]).sum() - ucty_wl[i] = (w**2 * ucty[sl] * dldx[sl] ** 2).sum() + flux_wl[i] = (w * flux[sl]).sum() + ucty_wl[i] = (w**2 * ucty[sl]).sum() else: fracw = bin_edges_pix[i + 1] - bin_edges_pix[i] - flux_wl[i] = fracw * flux[i1] * dldx[i1] - ucty_wl[i] = fracw**2 * ucty[i1] * dldx[i1] ** 2 + flux_wl[i] = fracw * flux[i1] + ucty_wl[i] = fracw**2 * ucty[i1] bin_widths_wav = np.diff(bin_edges_wav) flux_wl = flux_wl / bin_widths_wav * spectrum.flux.unit / self.unit From a2fb17669b3bcd4f8c2d9b909a009e6864bcd4b6 Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Sat, 5 Sep 2026 16:06:03 +0100 Subject: [PATCH 3/4] - Fixed `WavelengthSolution1D.resample` for spectra whose spectral axis does not start at pixel zero, and for solutions where the wavelength decreases with pixel number (the output is now always ascending in wavelength). - Changed `WavelengthSolution1D.resample` to propagate the input mask, copy the input metadata, and return no uncertainty instead of a fabricated zero-valued one when the input has none. --- CHANGES.rst | 11 +++- specreduce/tests/test_wavesol1d.py | 84 +++++++++++++++++++++++++++- specreduce/wavesol1d.py | 88 ++++++++++++++++++++---------- 3 files changed, 149 insertions(+), 34 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 322b2d41..463e70f8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -37,9 +37,14 @@ Bug Fixes the wavelength width of each pixel before dividing by the bin width. The output was labelled as a flux density per wavelength unit but was numerically off by the local dispersion, so a flat spectrum in counts per pixel stayed flat instead of falling where - the dispersion grows, and the integrated flux was not conserved. The flux density and - its uncertainty are now computed from the fractional pixel overlaps alone, with the - variance carrying the squared weights. [#XXX] + the dispersion grows, and the integrated flux was not conserved. [#XXX] + +- ``WavelengthSolution1D.resample()`` now bins spectra whose spectral axis does not start + at pixel zero correctly, supports solutions where the wavelength decreases with pixel + number (the output is always ascending in wavelength), propagates the input mask by + flagging every bin that received a contribution from a masked pixel, copies the input + metadata, returns no uncertainty instead of a fabricated zero-valued one when the + input has none, and rejects ``nbins=0`` as its error message already promised. [#XXX] API Changes ^^^^^^^^^^^ diff --git a/specreduce/tests/test_wavesol1d.py b/specreduce/tests/test_wavesol1d.py index 377ca9c2..43b8d682 100644 --- a/specreduce/tests/test_wavesol1d.py +++ b/specreduce/tests/test_wavesol1d.py @@ -127,7 +127,7 @@ def test_resample(mk_spectrum, mk_ws_with_transform, mk_ws_without_transform): # Resample a spectrum without uncertainty spectrum.uncertainty = None resampled = ws.resample(spectrum, nbins=50) - assert resampled.uncertainty is not None + assert resampled.uncertainty is None ws = mk_ws_without_transform with pytest.raises(ValueError, match="Wavelength solution not set."): @@ -136,6 +136,8 @@ def test_resample(mk_spectrum, mk_ws_with_transform, mk_ws_without_transform): ws = mk_ws_with_transform with pytest.raises(ValueError, match="Number of bins must be non-zero and positive"): ws.resample(mk_spectrum, nbins=-5) + with pytest.raises(ValueError, match="Number of bins must be non-zero and positive"): + ws.resample(mk_spectrum, nbins=0) def test_pix_to_wav(mk_ws_with_transform): @@ -450,6 +452,86 @@ def test_resample_propagates_variance(mk_ws_with_transform): np.testing.assert_allclose(var[ibin], expected, rtol=1e-6) +def _flat_spectrum(npix, flux=100.0, sigma=4.0, first_pixel=0): + return Spectrum( + flux=np.full(npix, flux) * u.count, + spectral_axis=(first_pixel + np.arange(npix)) * u.pix, + uncertainty=StdDevUncertainty(np.full(npix, sigma)), + ) + + +def test_resample_handles_offset_spectral_axis(): + """A spectrum whose spectral axis does not start at pixel zero must be binned correctly.""" + npix, first, dispersion = 100, 300, 2.4 + ws = WavelengthSolution1D( + models.Shift(0) | models.Polynomial1D(1, c0=5000, c1=dispersion), + (0, first + npix), + u.angstrom, + ) + flux = 100.0 + np.arange(npix) # a slope, so misindexing is visible + spectrum = _flat_spectrum(npix, first_pixel=first) + spectrum.flux[:] = flux * u.count + + result = ws.resample(spectrum) + np.testing.assert_allclose(result.flux.value, flux / dispersion, rtol=1e-9) + np.testing.assert_allclose(result.spectral_axis.value[0], ws.p2w(first), rtol=1e-12) + np.testing.assert_allclose(result.uncertainty.array, 4.0 / dispersion, rtol=1e-9) + + +def test_resample_handles_decreasing_dispersion(): + """Wavelength decreasing with pixel number gives an ascending, flux-conserving result.""" + npix, dispersion = 100, 2.4 + ws = WavelengthSolution1D( + models.Shift(0) | models.Polynomial1D(1, c0=9000, c1=-dispersion), (0, npix), u.angstrom + ) + spectrum = _flat_spectrum(npix) + spectrum.flux[:] = (100.0 + np.arange(npix)) * u.count + + result = ws.resample(spectrum) + assert np.all(np.diff(result.spectral_axis.value) > 0) + # pixel 0 is the reddest, so it lands in the last bin + np.testing.assert_allclose(result.flux.value, spectrum.flux.value[::-1] / dispersion, rtol=1e-9) + np.testing.assert_allclose(result.uncertainty.array, 4.0 / dispersion, rtol=1e-9) + + result = ws.resample(spectrum, nbins=npix // 4) + dlam = np.diff(result.spectral_axis.value)[0] + np.testing.assert_allclose( + (result.flux.value * dlam).sum(), spectrum.flux.value.sum(), rtol=1e-9 + ) + + +def test_resample_propagates_mask_and_meta(): + npix, dispersion = 100, 2.4 + ws = WavelengthSolution1D( + models.Shift(0) | models.Polynomial1D(1, c0=5000, c1=dispersion), (0, npix), u.angstrom + ) + spectrum = _flat_spectrum(npix) + assert ws.resample(spectrum).mask is None + + mask = np.zeros(npix, dtype=bool) + mask[40] = True + meta = {"OBJECT": "target", "HISTORY": ["extracted"]} + spectrum = Spectrum( + flux=spectrum.flux, spectral_axis=spectrum.spectral_axis, mask=mask, meta=meta + ) + + # one bin per pixel: exactly the masked pixel's bin is flagged + result = ws.resample(spectrum) + np.testing.assert_array_equal(result.mask, mask) + + # bins half a pixel out of phase with the pixels: the two overlapping bins are flagged + edges = ws.p2w(np.arange(npix + 1)) # pixel centres act as bin edges + result = ws.resample(spectrum, bin_edges=edges) + expected = np.zeros(npix, dtype=bool) + expected[39:41] = True + np.testing.assert_array_equal(result.mask, expected) + + assert result.meta == meta + assert result.meta is not meta + result.meta["HISTORY"].append("resampled") + assert meta["HISTORY"] == ["extracted"] + + def test_resample_returns_flux_density_per_wavelength(): """ A flat spectrum in counts per pixel becomes counts per wavelength unit divided by diff --git a/specreduce/wavesol1d.py b/specreduce/wavesol1d.py index b15d9f99..5a418ed1 100644 --- a/specreduce/wavesol1d.py +++ b/specreduce/wavesol1d.py @@ -1,8 +1,8 @@ import warnings +from copy import deepcopy from functools import cached_property from pathlib import Path from typing import Callable -from copy import deepcopy import asdf import astropy.units as u @@ -641,7 +641,15 @@ def resample( pixels it overlaps, weighted by the overlapping fraction of each pixel, and the total is divided by the bin width, so the output is a flux density per wavelength unit. The binning is exact and conserves the integrated flux. The variance is propagated - with the squared weights, assuming independent pixel noise. + with the squared weights, assuming independent pixel noise. The spectral axis of the + input must be in pixels with unit spacing but need not start at zero, and the + wavelength may increase or decrease with pixel number; the output is always + ascending in wavelength. + + Bins that are not aligned with the pixels share pixels with their neighbours, so + the resampled uncertainties, while correct for each bin, are correlated between + neighbouring bins. Sums or fits that treat the bins as independent underestimate + the uncertainty; resampling to a grid finer than the pixels makes this worse. Parameters ---------- @@ -666,10 +674,12 @@ def resample( Returns ------- 1D spectrum binned to the specified wavelength bins, with the flux in units of - the input flux unit per wavelength unit and the uncertainty in the same - uncertainty class as the input. + the input flux unit per wavelength unit. The uncertainty, if the input has one, + is returned in the same uncertainty class as the input; the mask, if the input + has one, flags every bin that received a contribution from a masked pixel; and + the metadata is a copy of the input metadata. """ - if nbins is not None and nbins < 0: + if nbins is not None and nbins <= 0: raise ValueError("Number of bins must be non-zero and positive.") if self._p2w is None: @@ -677,51 +687,69 @@ def resample( flux = spectrum.flux.value pixels = spectrum.spectral_axis.value + npix = flux.size - if spectrum.uncertainty is not None: + has_uncertainty = spectrum.uncertainty is not None + if has_uncertainty: ucty = spectrum.uncertainty.represent_as(VarianceUncertainty).array ucty_type = type(spectrum.uncertainty) - else: - ucty = np.zeros_like(flux) - ucty_type = VarianceUncertainty - npix = flux.size - nbins = npix if nbins is None else nbins - if wlbounds is None: - l1, l2 = self.p2w(pixels[[0, -1]] + np.array([-0.5, 0.5])) - else: - l1, l2 = wlbounds + has_mask = spectrum.mask is not None + mask = np.asarray(spectrum.mask, dtype=bool) if has_mask else None + nbins = npix if nbins is None else nbins if bin_edges is not None: - bin_edges_wav = np.asarray(bin_edges) + bin_edges_wav = np.sort(np.asarray(bin_edges, dtype=float)) nbins = bin_edges_wav.size - 1 else: + if wlbounds is None: + l1, l2 = sorted(self.p2w(pixels[[0, -1]] + np.array([-0.5, 0.5]))) + else: + l1, l2 = sorted(wlbounds) bin_edges_wav = np.linspace(l1, l2, num=nbins + 1) - - bin_edges_pix = np.clip(self.w2p(bin_edges_wav) + 0.5, 0, npix - 1e-12) - bin_edge_ix = np.floor(bin_edges_pix).astype(int) - bin_edge_w = bin_edges_pix - bin_edge_ix bin_centers_wav = 0.5 * (bin_edges_wav[:-1] + bin_edges_wav[1:]) + + # Bin edges in array-index space, where pixel j spans [j, j + 1). The spectral axis + # need not start at zero, and the wavelength may decrease with pixel number, in + # which case the left and right edges of a bin swap places in pixel space. + x = np.clip(self.w2p(bin_edges_wav) + 0.5 - pixels[0], 0, npix - 1e-12) + x_left, x_right = np.minimum(x[:-1], x[1:]), np.maximum(x[:-1], x[1:]) + i_left, i_right = np.floor(x_left).astype(int), np.floor(x_right).astype(int) + flux_wl = np.zeros(nbins) - ucty_wl = np.zeros(nbins) + ucty_wl = np.zeros(nbins) if has_uncertainty else None + mask_wl = np.zeros(nbins, dtype=bool) if has_mask else None weights = np.zeros(npix) for i in range(nbins): - i1, i2 = bin_edge_ix[i : i + 2] - weights[:] = 0.0 + i1, i2 = i_left[i], i_right[i] if i1 != i2: + weights[:] = 0.0 weights[i1 + 1 : i2] = 1.0 - weights[i1] = 1 - bin_edge_w[i] - weights[i2] = bin_edge_w[i + 1] + weights[i1] = 1.0 - (x_left[i] - i1) + weights[i2] = x_right[i] - i2 sl = slice(i1, i2 + 1) w = weights[sl] flux_wl[i] = (w * flux[sl]).sum() - ucty_wl[i] = (w**2 * ucty[sl]).sum() + if has_uncertainty: + ucty_wl[i] = (w**2 * ucty[sl]).sum() + if has_mask: + mask_wl[i] = (mask[sl] & (w > 0)).any() else: - fracw = bin_edges_pix[i + 1] - bin_edges_pix[i] + fracw = x_right[i] - x_left[i] flux_wl[i] = fracw * flux[i1] - ucty_wl[i] = fracw**2 * ucty[i1] + if has_uncertainty: + ucty_wl[i] = fracw**2 * ucty[i1] + if has_mask: + mask_wl[i] = mask[i1] & (fracw > 0) bin_widths_wav = np.diff(bin_edges_wav) flux_wl = flux_wl / bin_widths_wav * spectrum.flux.unit / self.unit - ucty_wl = VarianceUncertainty(ucty_wl / bin_widths_wav**2).represent_as(ucty_type) - return Spectrum(flux_wl, bin_centers_wav * self.unit, uncertainty=ucty_wl) + if has_uncertainty: + ucty_wl = VarianceUncertainty(ucty_wl / bin_widths_wav**2).represent_as(ucty_type) + return Spectrum( + flux_wl, + bin_centers_wav * self.unit, + uncertainty=ucty_wl, + mask=mask_wl, + meta=deepcopy(spectrum.meta), + ) From 4c56a82fbc11be0b443e5e728c5c84e1f6ce2c4d Mon Sep 17 00:00:00 2001 From: Hannu Parviainen Date: Sat, 5 Sep 2026 16:16:44 +0100 Subject: [PATCH 4/4] - Changed `WavelengthSolution1D.resample` to store the bin edges actually used on the output spectral axis, so spectral_axis.bin_edges is exact for non-uniform bin grids instead of being inferred from the bin centres. --- CHANGES.rst | 4 ++++ specreduce/tests/test_wavesol1d.py | 22 ++++++++++++++++++++++ specreduce/wavesol1d.py | 12 +++++++----- 3 files changed, 33 insertions(+), 5 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 463e70f8..cadf2ac9 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -46,6 +46,10 @@ Bug Fixes metadata, returns no uncertainty instead of a fabricated zero-valued one when the input has none, and rejects ``nbins=0`` as its error message already promised. [#XXX] +- The spectrum returned by ``WavelengthSolution1D.resample()`` now carries the bin edges + actually used on its spectral axis, so ``spectral_axis.bin_edges`` is exact for + non-uniform ``bin_edges`` grids instead of being inferred from the bin centres. [#XXX] + API Changes ^^^^^^^^^^^ diff --git a/specreduce/tests/test_wavesol1d.py b/specreduce/tests/test_wavesol1d.py index 43b8d682..3b31b61d 100644 --- a/specreduce/tests/test_wavesol1d.py +++ b/specreduce/tests/test_wavesol1d.py @@ -532,6 +532,28 @@ def test_resample_propagates_mask_and_meta(): assert meta["HISTORY"] == ["extracted"] +def test_resample_keeps_explicit_bin_edges(): + """Non-uniform bin edges must be carried on the output spectral axis as given.""" + npix, dispersion = 100, 2.4 + ws = WavelengthSolution1D( + models.Shift(0) | models.Polynomial1D(1, c0=5000, c1=dispersion), (0, npix), u.angstrom + ) + spectrum = _flat_spectrum(npix) + lo, hi = ws.p2w(np.array([-0.5, npix - 0.5])) + edges = lo + (hi - lo) * np.linspace(0, 1, 9) ** 2 # widths growing along the axis + + result = ws.resample(spectrum, bin_edges=edges) + np.testing.assert_allclose(result.spectral_axis.bin_edges.value, edges, rtol=1e-12) + np.testing.assert_allclose( + result.spectral_axis.value, 0.5 * (edges[:-1] + edges[1:]), rtol=1e-12 + ) + # integrating over the true bin widths conserves the counts + widths = np.diff(result.spectral_axis.bin_edges.value) + np.testing.assert_allclose( + (result.flux.value * widths).sum(), spectrum.flux.value.sum(), rtol=1e-9 + ) + + def test_resample_returns_flux_density_per_wavelength(): """ A flat spectrum in counts per pixel becomes counts per wavelength unit divided by diff --git a/specreduce/wavesol1d.py b/specreduce/wavesol1d.py index 5a418ed1..e9c33b44 100644 --- a/specreduce/wavesol1d.py +++ b/specreduce/wavesol1d.py @@ -18,7 +18,7 @@ from scipy import optimize from scipy.interpolate import interp1d -from specutils import Spectrum +from specutils import Spectrum, SpectralAxis __all__ = ["WavelengthSolution1D"] @@ -640,7 +640,9 @@ def resample( produced by the extraction methods). Each wavelength bin receives the flux of the pixels it overlaps, weighted by the overlapping fraction of each pixel, and the total is divided by the bin width, so the output is a flux density per wavelength unit. - The binning is exact and conserves the integrated flux. The variance is propagated + The bin edges are stored on the output spectral axis, so the integrated flux is + recovered exactly as ``(flux * np.diff(spectral_axis.bin_edges)).sum()`` for any + bin grid. The binning is exact and conserves the integrated flux. The variance is propagated with the squared weights, assuming independent pixel noise. The spectral axis of the input must be in pixels with unit spacing but need not start at zero, and the wavelength may increase or decrease with pixel number; the output is always @@ -674,7 +676,8 @@ def resample( Returns ------- 1D spectrum binned to the specified wavelength bins, with the flux in units of - the input flux unit per wavelength unit. The uncertainty, if the input has one, + the input flux unit per wavelength unit and a spectral axis that carries the bin + edges actually used. The uncertainty, if the input has one, is returned in the same uncertainty class as the input; the mask, if the input has one, flags every bin that received a contribution from a masked pixel; and the metadata is a copy of the input metadata. @@ -706,7 +709,6 @@ def resample( else: l1, l2 = sorted(wlbounds) bin_edges_wav = np.linspace(l1, l2, num=nbins + 1) - bin_centers_wav = 0.5 * (bin_edges_wav[:-1] + bin_edges_wav[1:]) # Bin edges in array-index space, where pixel j spans [j, j + 1). The spectral axis # need not start at zero, and the wavelength may decrease with pixel number, in @@ -748,7 +750,7 @@ def resample( ucty_wl = VarianceUncertainty(ucty_wl / bin_widths_wav**2).represent_as(ucty_type) return Spectrum( flux_wl, - bin_centers_wav * self.unit, + SpectralAxis(bin_edges_wav * self.unit, bin_specification="edges"), uncertainty=ucty_wl, mask=mask_wl, meta=deepcopy(spectrum.meta),