From ea5e6a282659ffa6822bdd360d1eef9777a98af0 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 22:50:44 -0600 Subject: [PATCH 01/18] synth_data refactor: add SynthImage builder core with background + outputs --- specreduce/tests/test_synth_data.py | 197 ++++------------------------ specreduce/utils/synth_data.py | 189 +++++++++++++++++++++++++- 2 files changed, 208 insertions(+), 178 deletions(-) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index b8bcb870..9dd9b9e0 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -1,187 +1,40 @@ +import numpy as np import pytest from astropy import units as u from astropy.modeling import models from astropy.nddata import CCDData -from astropy.wcs import WCS -from specreduce.utils.synth_data import make_2d_trace_image, make_2d_arc_image, make_2d_spec_image +from specreduce.utils.synth_data import SynthImage -def test_make_2d_trace_image(): - ccdim = make_2d_trace_image( - nx=3000, - ny=1000, - background=5, - trace_center=None, - trace_order=3, - trace_coeffs={'c0': 0, 'c1': 50, 'c2': 100}, - profile=models.Gaussian1D(amplitude=100, stddev=10) - ) - assert ccdim.data.shape == (1000, 3000) - assert isinstance(ccdim, CCDData) +def test_empty_image_is_zeros(): + arr = SynthImage(nx=50, ny=20).to_array() + assert arr.shape == (20, 50) + assert np.all(arr == 0) -@pytest.mark.remote_data -@pytest.mark.filterwarnings("ignore:No observer defined on WCS") -def test_make_2d_arc_image_defaults(): - ccdim = make_2d_arc_image() - assert isinstance(ccdim, CCDData) +def test_add_background_constant(): + arr = SynthImage(nx=50, ny=20).add_background(7).to_array() + assert arr.shape == (20, 50) + assert np.allclose(arr, 7) -@pytest.mark.remote_data -@pytest.mark.filterwarnings("ignore:No observer defined on WCS") -def test_make_2d_arc_pass_wcs(): - nx = 3000 - ny = 1000 - wave_unit = u.Angstrom - extent = [3000, 6000] +def test_add_background_stacks(): + arr = SynthImage(nx=10, ny=10).add_background(3).add_background(4).to_array() + assert np.allclose(arr, 7) - # test passing a valid WCS with dispersion along X - wcs = WCS(naxis=2) - wcs.wcs.ctype[0] = 'WAVE' - wcs.wcs.ctype[1] = 'PIXEL' - wcs.wcs.cunit[0] = wave_unit - wcs.wcs.cunit[1] = u.pixel - wcs.wcs.crval[0] = extent[0] - wcs.wcs.cdelt[0] = (extent[1] - extent[0]) / nx - wcs.wcs.crval[1] = 0 - wcs.wcs.cdelt[1] = 1 - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=None, - wave_unit=None, - wcs=wcs - ) - assert ccdim.data.shape == (1000, 3000) - assert isinstance(ccdim, CCDData) +def test_immutability(): + base = SynthImage(nx=10, ny=10).add_background(5) + derived = base.add_background(2) + assert np.allclose(base.to_array(), 5) + assert np.allclose(derived.to_array(), 7) + assert base is not derived - # test passing a tilt model - tilt_model = models.Chebyshev1D(degree=2, c0=50, c1=0, c2=100) - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=None, - wave_unit=None, - wcs=wcs, - tilt_func=tilt_model - ) - assert ccdim.data.shape == (1000, 3000) - assert isinstance(ccdim, CCDData) - # make sure WCS without spectral axis gets rejected - wcs.wcs.ctype[0] = 'PIXEL' - assert wcs.spectral.naxis == 0 - with pytest.raises(ValueError, match='Provided WCS must have a spectral axis'): - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=None, - wave_unit=None, - wcs=wcs - ) - - # test passing valid WCS with dispersion along Y while using air wavelengths - wcs = WCS(naxis=2) - wcs.wcs.ctype[1] = 'AWAV' - wcs.wcs.ctype[0] = 'PIXEL' - wcs.wcs.cunit[1] = wave_unit - wcs.wcs.cunit[0] = u.pixel - wcs.wcs.crval[1] = extent[0] - wcs.wcs.cdelt[1] = (extent[1] - extent[0]) / nx - wcs.wcs.crval[0] = 0 - wcs.wcs.cdelt[0] = 1 - - ccdim = make_2d_arc_image( - nx=ny, - ny=nx, - extent=None, - wave_unit=None, - wave_air=True, - wcs=wcs, - tilt_func=tilt_model - ) - assert ccdim.data.shape == (3000, 1000) - assert isinstance(ccdim, CCDData) - - # make sure no WCS and no extent gets rejected - with pytest.raises(ValueError, match='Must specify either a wavelength extent or a WCS'): - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=None, - wave_unit=None, - wcs=None - ) - - # make sure if extent is provided, it has the right length - with pytest.raises(ValueError, match='Wavelength extent must be of length 2'): - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=[1, 2, 3], - wave_unit=None, - wcs=None - ) - - # make sure a 1D WCS gets rejected - wcs = WCS(naxis=1) - wcs.wcs.ctype[0] = 'WAVE' - wcs.wcs.cunit[0] = wave_unit - wcs.wcs.crval[0] = extent[0] - wcs.wcs.cdelt[0] = (extent[1] - extent[0]) / nx - - with pytest.raises(ValueError, match='WCS must have NAXIS=2 for a 2D image'): - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=None, - wave_unit=None, - wcs=wcs - ) - - # make sure a WCS with no spectral axis gets rejected - wcs = WCS(naxis=2) - wcs.wcs.ctype[1] = 'PIXEL' - wcs.wcs.ctype[0] = 'PIXEL' - wcs.wcs.cunit[1] = u.pixel - wcs.wcs.cunit[0] = u.pixel - wcs.wcs.crval[1] = extent[0] - wcs.wcs.cdelt[1] = (extent[1] - extent[0]) / nx - wcs.wcs.crval[0] = 0 - wcs.wcs.cdelt[0] = 1 - - with pytest.raises(ValueError, match='Provided WCS must have a spectral axis'): - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=None, - wave_unit=None, - wcs=wcs - ) - - # make sure invalid wave_unit is caught - with pytest.raises(ValueError, match='Wavelength unit must be a length unit'): - ccdim = make_2d_arc_image( - nx=nx, - ny=ny, - extent=[100, 300], - wave_unit=u.pixel - ) - - # make sure a non-polynomial tilt_func gets rejected - with pytest.raises( - ValueError, - match='The only tilt functions currently supported are 1D polynomials' - ): - ccdim = make_2d_arc_image( - tilt_func=models.Gaussian1D - ) - - -@pytest.mark.remote_data -@pytest.mark.filterwarnings("ignore:No observer defined on WCS") -def test_make_2d_spec_image_defaults(): - ccdim = make_2d_spec_image() - assert isinstance(ccdim, CCDData) +def test_to_ccddata_no_wcs(): + ccd = SynthImage(nx=30, ny=15).add_background(5).to_ccddata() + assert isinstance(ccd, CCDData) + assert ccd.unit == u.count + assert ccd.data.shape == (15, 30) + assert ccd.wcs is None diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index a270db93..6dddb313 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -1,6 +1,6 @@ # Licensed under a 3-clause BSD style license - see ../../licenses/LICENSE.rst import warnings -from typing import Sequence +from dataclasses import dataclass, field import numpy as np from astropy import units as u @@ -9,14 +9,191 @@ from astropy.stats import gaussian_fwhm_to_sigma from astropy.wcs import WCS +from specutils import Spectrum + from specreduce.calibration_data import load_pypeit_calibration_lines __all__ = [ - 'make_2d_trace_image', - 'make_2d_arc_image', - 'make_2d_spec_image' + "SynthImage", + "make_2d_trace_image", + "make_2d_arc_image", + "make_2d_spec_image", ] +_ALLOWED_TILT = ( + models.Legendre1D, + models.Chebyshev1D, + models.Polynomial1D, + models.Hermite1D, +) + + +@dataclass(frozen=True) +class _RenderContext: + """Shared geometry passed to every layer's ``render`` method.""" + nx: int + ny: int + xx: np.ndarray + yy: np.ndarray + wcs: WCS | None + disp_axis: int + + +@dataclass(frozen=True) +class BackgroundLayer: + """A constant additive background level in counts.""" + level: float + + def render(self, ctx: _RenderContext) -> np.ndarray: + return np.full((ctx.ny, ctx.nx), float(self.level)) + + +class SynthImage: + """ + Immutable, composable builder for synthetic 2D spectroscopic images. + + Build an image by chaining ``add_*`` methods, then render it with one of the + ``to_*`` terminal methods. Each ``add_*`` returns a *new* ``SynthImage``; the + original is never mutated, so a base configuration can be safely branched. + + Parameters + ---------- + nx + Size of the image along the X (dispersion) axis. + ny + Size of the image along the Y (spatial) axis. + wcs + Optional 2D WCS with a single spectral axis. If not provided and arc + layers are present, a linear ``WAVE``/``PIXEL`` WCS is built from + ``extent`` and ``wave_unit``. + extent + Beginning and end wavelengths used to build a default WCS when ``wcs`` + is not supplied and arc layers are present. + wave_unit + Wavelength unit for the default WCS. + seed + Seed for the random number generator used by the noise layers. If + ``None``, noise is non-deterministic. + """ + + def __init__( + self, + nx: int = 3000, + ny: int = 1000, + wcs: WCS | None = None, + extent=(3500, 7000), + wave_unit: u.Unit = u.Angstrom, + seed: int | None = None, + ): + self.nx = nx + self.ny = ny + self._wcs = wcs + self._extent = extent + self._wave_unit = wave_unit + self._seed = seed + self._layers = () + self._poisson = False + self._read_noise = None + + def _clone(self, **changes) -> "SynthImage": + new = SynthImage.__new__(SynthImage) + new.nx = self.nx + new.ny = self.ny + new._wcs = self._wcs + new._extent = self._extent + new._wave_unit = self._wave_unit + new._seed = self._seed + new._layers = self._layers + new._poisson = self._poisson + new._read_noise = self._read_noise + for key, value in changes.items(): + setattr(new, key, value) + return new + + def add_background(self, level: float) -> "SynthImage": + """Add a constant background level (counts).""" + return self._clone(_layers=self._layers + (BackgroundLayer(level),)) + + def _resolve_wcs(self): + has_arc = any(isinstance(layer, ArcLayer) for layer in self._layers) + if self._wcs is not None: + wcs = self._wcs + if wcs.spectral.naxis != 1: + raise ValueError("Provided WCS must have a spectral axis.") + if wcs.naxis != 2: + raise ValueError("WCS must have NAXIS=2 for a 2D image.") + elif has_arc: + if self._extent is None: + raise ValueError("Must specify either a wavelength extent or a WCS.") + if len(self._extent) != 2: + raise ValueError("Wavelength extent must be of length 2.") + if u.get_physical_type(self._wave_unit) != "length": + raise ValueError("Wavelength unit must be a length unit.") + wcs = WCS(naxis=2) + wcs.wcs.ctype[0] = "WAVE" + wcs.wcs.ctype[1] = "PIXEL" + wcs.wcs.cunit[0] = self._wave_unit + wcs.wcs.cunit[1] = u.pixel + wcs.wcs.crval[0] = self._extent[0] + wcs.wcs.cdelt[0] = (self._extent[1] - self._extent[0]) / self.nx + wcs.wcs.crval[1] = 0 + wcs.wcs.cdelt[1] = 1 + else: + wcs = None + + if wcs is None: + disp_axis = 1 + else: + is_spectral = [a["coordinate_type"] == "spectral" for a in wcs.get_axis_types()] + disp_axis = 0 if is_spectral[0] else 1 + return wcs, disp_axis + + def _render(self): + wcs, disp_axis = self._resolve_wcs() + x = np.arange(self.nx) + y = np.arange(self.ny) + xx, yy = np.meshgrid(x, y) + ctx = _RenderContext(self.nx, self.ny, xx, yy, wcs, disp_axis) + + signal = np.zeros((self.ny, self.nx)) + for layer in self._layers: + signal = signal + layer.render(ctx) + + rng = np.random.default_rng(self._seed) + if self._poisson: + from photutils.datasets import apply_poisson_noise + signal = apply_poisson_noise(signal, seed=rng) + if self._read_noise is not None: + signal = signal + rng.normal(0.0, self._read_noise, size=signal.shape) + + return signal, wcs + + def to_array(self) -> np.ndarray: + """Render and return the image as a plain ``numpy.ndarray`` (counts).""" + return self._render()[0] + + def to_ccddata(self) -> CCDData: + """Render and return the image as a `~astropy.nddata.CCDData`.""" + data, wcs = self._render() + return CCDData(data, unit=u.count, wcs=wcs) + + def to_spectrum(self) -> Spectrum: + """Render and return the image as a `~specutils.Spectrum`.""" + data, wcs = self._render() + if wcs is not None: + return Spectrum(flux=data * u.count, wcs=wcs) + return Spectrum(flux=data * u.count, spectral_axis_index=data.ndim - 1) + + +@dataclass(frozen=True) +class SourceLayer: + pass # implemented in Task 2 + + +@dataclass(frozen=True) +class ArcLayer: + pass # implemented in Task 3 + def make_2d_trace_image( nx: int = 3000, @@ -92,7 +269,7 @@ def make_2d_arc_image( nx: int = 3000, ny: int = 1000, wcs: WCS | None = None, - extent: Sequence[int | float] = (3500, 7000), + extent=(3500, 7000), wave_unit: u.Unit = u.Angstrom, wave_air: bool = False, background: int | float = 5, @@ -345,7 +522,7 @@ def make_2d_spec_image( nx: int = 3000, ny: int = 1000, wcs: WCS | None = None, - extent: Sequence[int | float] = (6500, 9500), + extent=(6500, 9500), wave_unit: u.Unit = u.Angstrom, wave_air: bool = False, background: int | float = 5, From e0a3b9e7f77f792896936972867c3fc0eed6754c Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 22:56:41 -0600 Subject: [PATCH 02/18] synth_data refactor: add SourceLayer and stackable add_source --- specreduce/tests/test_synth_data.py | 24 ++++++++++++++++++++ specreduce/utils/synth_data.py | 34 ++++++++++++++++++++++++++++- 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index 9dd9b9e0..c5698b48 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -38,3 +38,27 @@ def test_to_ccddata_no_wcs(): assert ccd.unit == u.count assert ccd.data.shape == (15, 30) assert ccd.wcs is None + + +def test_add_source_adds_flux_at_trace(): + base = SynthImage(nx=200, ny=100).add_background(5) + src = base.add_source(profile=models.Gaussian1D(amplitude=100, stddev=10)) + base_arr = base.to_array() + src_arr = src.to_array() + # source adds flux somewhere; total flux strictly increases + assert src_arr.sum() > base_arr.sum() + assert src_arr.max() > base_arr.max() + + +def test_add_source_default_profile(): + arr = SynthImage(nx=200, ny=100).add_source().to_array() + assert arr.shape == (100, 200) + assert arr.max() > 0 + + +def test_add_source_stackable(): + one = SynthImage(nx=200, ny=100).add_source( + profile=models.Gaussian1D(amplitude=100, stddev=10) + ) + two = one.add_source(profile=models.Gaussian1D(amplitude=100, stddev=10)) + assert two.to_array().sum() > one.to_array().sum() diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 6dddb313..c5a2baa2 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -114,6 +114,19 @@ def add_background(self, level: float) -> "SynthImage": """Add a constant background level (counts).""" return self._clone(_layers=self._layers + (BackgroundLayer(level),)) + def add_source( + self, + profile: Model | None = None, + trace_center: float | None = None, + trace_order: int = 3, + trace_coeffs: dict | None = None, + ) -> "SynthImage": + """Add a continuum source with a Chebyshev-traced spatial profile.""" + if profile is None: + profile = models.Moffat1D(amplitude=10, alpha=0.1) + layer = SourceLayer(profile, trace_center, trace_order, trace_coeffs) + return self._clone(_layers=self._layers + (layer,)) + def _resolve_wcs(self): has_arc = any(isinstance(layer, ArcLayer) for layer in self._layers) if self._wcs is not None: @@ -187,7 +200,26 @@ def to_spectrum(self) -> Spectrum: @dataclass(frozen=True) class SourceLayer: - pass # implemented in Task 2 + """A continuum source whose spatial profile follows a Chebyshev trace. + + The dispersion axis is the X (column) axis, matching the historical + ``make_2d_trace_image`` behaviour (the trace ignores any WCS). + """ + profile: Model + trace_center: float | None = None + trace_order: int = 3 + trace_coeffs: dict | None = None + + def render(self, ctx: _RenderContext) -> np.ndarray: + trace_center = ctx.ny / 2 if self.trace_center is None else self.trace_center + trace_coeffs = ( + {"c0": 0, "c1": 50, "c2": 100} + if self.trace_coeffs is None + else self.trace_coeffs + ) + trace_mod = models.Chebyshev1D(degree=self.trace_order, **trace_coeffs) + trace = ctx.yy - trace_center + trace_mod(ctx.xx / ctx.nx) + return self.profile(trace) @dataclass(frozen=True) From 31da4ad7aadfc9b4beb2de07f01db0aef9a7a323 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 23:16:03 -0600 Subject: [PATCH 03/18] synth_data refactor: add ArcLayer, add_arcs, and WCS resolution --- specreduce/tests/test_synth_data.py | 91 +++++++++++++++++++++++++++++ specreduce/utils/synth_data.py | 62 +++++++++++++++++++- 2 files changed, 152 insertions(+), 1 deletion(-) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index c5698b48..a9dffc83 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -3,6 +3,7 @@ from astropy import units as u from astropy.modeling import models from astropy.nddata import CCDData +from astropy.wcs import WCS from specreduce.utils.synth_data import SynthImage @@ -62,3 +63,93 @@ def test_add_source_stackable(): ) two = one.add_source(profile=models.Gaussian1D(amplitude=100, stddev=10)) assert two.to_array().sum() > one.to_array().sum() + + +def _linear_wcs(nx, extent=(3000, 6000), wave_unit=u.Angstrom): + wcs = WCS(naxis=2) + wcs.wcs.ctype[0] = "WAVE" + wcs.wcs.ctype[1] = "PIXEL" + wcs.wcs.cunit[0] = wave_unit + wcs.wcs.cunit[1] = u.pixel + wcs.wcs.crval[0] = extent[0] + wcs.wcs.cdelt[0] = (extent[1] - extent[0]) / nx + wcs.wcs.crval[1] = 0 + wcs.wcs.cdelt[1] = 1 + return wcs + + +# --- local (no network) validation tests --- + +def test_arc_requires_wcs_or_extent(): + img = SynthImage(nx=100, ny=50, extent=None).add_arcs() + with pytest.raises(ValueError, match="Must specify either a wavelength extent or a WCS"): + img.to_array() + + +def test_arc_extent_wrong_length(): + img = SynthImage(nx=100, ny=50, extent=[1, 2, 3]).add_arcs() + with pytest.raises(ValueError, match="Wavelength extent must be of length 2"): + img.to_array() + + +def test_arc_bad_wave_unit(): + img = SynthImage(nx=100, ny=50, extent=[100, 300], wave_unit=u.pixel).add_arcs() + with pytest.raises(ValueError, match="Wavelength unit must be a length unit"): + img.to_array() + + +def test_arc_bad_tilt_func(): + wcs = _linear_wcs(100) + img = SynthImage(nx=100, ny=50, wcs=wcs).add_arcs(tilt_func=models.Gaussian1D) + with pytest.raises( + ValueError, match="The only tilt functions currently supported are 1D polynomials" + ): + img.to_array() + + +def test_provided_wcs_must_have_spectral_axis(): + wcs = WCS(naxis=2) + wcs.wcs.ctype[0] = "PIXEL" + wcs.wcs.ctype[1] = "PIXEL" + img = SynthImage(nx=100, ny=50, wcs=wcs).add_arcs() + with pytest.raises(ValueError, match="Provided WCS must have a spectral axis"): + img.to_array() + + +def test_provided_wcs_must_be_2d(): + wcs = WCS(naxis=1) + wcs.wcs.ctype[0] = "WAVE" + wcs.wcs.cunit[0] = u.Angstrom + img = SynthImage(nx=100, ny=50, wcs=wcs).add_arcs() + with pytest.raises(ValueError, match="WCS must have NAXIS=2 for a 2D image"): + img.to_array() + + +# --- network tests (load pypeit line lists) --- + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_arc_image_with_extent_builds_wcs(): + ccd = SynthImage(nx=300, ny=100, extent=(3500, 7000)).add_arcs(["HeI"]).to_ccddata() + assert isinstance(ccd, CCDData) + assert ccd.data.shape == (100, 300) + assert ccd.wcs is not None + assert ccd.data.max() > 0 + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_arc_image_supplied_wcs_and_tilt(): + wcs = _linear_wcs(300) + tilt = models.Chebyshev1D(degree=2, c0=50, c1=0, c2=100) + ccd = SynthImage(nx=300, ny=100, wcs=wcs).add_arcs(["HeI"], tilt_func=tilt).to_ccddata() + assert ccd.data.shape == (100, 300) + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_arcs_stackable(): + wcs = _linear_wcs(300) + one = SynthImage(nx=300, ny=100, wcs=wcs).add_arcs(["HeI"]) + two = one.add_arcs(["NeI"]) + assert two.to_array().sum() > one.to_array().sum() diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index c5a2baa2..47ccdd58 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -127,6 +127,24 @@ def add_source( layer = SourceLayer(profile, trace_center, trace_order, trace_coeffs) return self._clone(_layers=self._layers + (layer,)) + def add_arcs( + self, + linelists=("HeI",), + line_fwhm: float = 5.0, + amplitude_scale: float = 1.0, + wave_air: bool = False, + tilt_func: Model | None = None, + ) -> "SynthImage": + """Add emission lines from one or more pypeit calibration line lists.""" + if tilt_func is None: + tilt_func = models.Legendre1D(degree=0) + if isinstance(linelists, str): + linelists = (linelists,) + layer = ArcLayer( + tuple(linelists), line_fwhm, amplitude_scale, wave_air, tilt_func + ) + return self._clone(_layers=self._layers + (layer,)) + def _resolve_wcs(self): has_arc = any(isinstance(layer, ArcLayer) for layer in self._layers) if self._wcs is not None: @@ -224,7 +242,49 @@ def render(self, ctx: _RenderContext) -> np.ndarray: @dataclass(frozen=True) class ArcLayer: - pass # implemented in Task 3 + """Emission lines from one or more pypeit calibration line lists. + + Requires a resolvable WCS (supplied on ``SynthImage`` or built from its + ``extent``). ``tilt_func`` applies a cross-dispersion tilt to simulate + curved lines. + """ + linelists: tuple + line_fwhm: float = 5.0 + amplitude_scale: float = 1.0 + wave_air: bool = False + tilt_func: Model = field(default_factory=lambda: models.Legendre1D(degree=0)) + + def render(self, ctx: _RenderContext) -> np.ndarray: + xx, yy = ctx.xx, ctx.yy + if self.tilt_func is not None: + if not isinstance(self.tilt_func, _ALLOWED_TILT): + raise ValueError( + "The only tilt functions currently supported are 1D polynomials " + "from astropy.models." + ) + if ctx.disp_axis == 0: + xx = xx + self.tilt_func((yy - ctx.ny / 2) / ctx.ny) + else: + yy = yy + self.tilt_func((xx - ctx.nx / 2) / ctx.nx) + + z = np.zeros((ctx.ny, ctx.nx)) + linelist = load_pypeit_calibration_lines(list(self.linelists), wave_air=self.wave_air) + if linelist is not None: + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", message="No observer defined on WCS.*") + line_disp_positions = ctx.wcs.spectral.world_to_pixel(linelist["wavelength"]) + line_sigma = gaussian_fwhm_to_sigma * self.line_fwhm + for line_pos, ampl in zip(line_disp_positions, linelist["amplitude"]): + line_mod = models.Gaussian1D( + amplitude=ampl * self.amplitude_scale, + mean=line_pos, + stddev=line_sigma, + ) + if ctx.disp_axis == 0: + z += line_mod(xx) + else: + z += line_mod(yy) + return z def make_2d_trace_image( From 1137d6e73e802b213c42a1e1fcd22c323addce29 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 23:19:41 -0600 Subject: [PATCH 04/18] synth_data refactor: add add_skylines wrapper for OH airglow lines --- specreduce/tests/test_synth_data.py | 9 +++++++++ specreduce/utils/synth_data.py | 4 ++++ 2 files changed, 13 insertions(+) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index a9dffc83..c00ed4da 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -153,3 +153,12 @@ def test_arcs_stackable(): one = SynthImage(nx=300, ny=100, wcs=wcs).add_arcs(["HeI"]) two = one.add_arcs(["NeI"]) assert two.to_array().sum() > one.to_array().sum() + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_add_skylines_matches_add_arcs(): + wcs = _linear_wcs(300, extent=(6500, 9500)) + sky = SynthImage(nx=300, ny=100, wcs=wcs, seed=1).add_skylines("OH_GMOS") + arc = SynthImage(nx=300, ny=100, wcs=wcs, seed=1).add_arcs("OH_GMOS") + assert np.array_equal(sky.to_array(), arc.to_array()) diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 47ccdd58..3a123f82 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -145,6 +145,10 @@ def add_arcs( ) return self._clone(_layers=self._layers + (layer,)) + def add_skylines(self, linelists="OH_GMOS", **kwargs) -> "SynthImage": + """Add night-sky airglow emission lines (OH lists), wrapping ``add_arcs``.""" + return self.add_arcs(linelists, **kwargs) + def _resolve_wcs(self): has_arc = any(isinstance(layer, ArcLayer) for layer in self._layers) if self._wcs is not None: From 97e0cd8951b0f07bdc8a3e8f37a5a0c5b406b676 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 23:23:46 -0600 Subject: [PATCH 05/18] synth_data refactor: add seeded Poisson and read noise (add_rdnoise alias) --- specreduce/tests/test_synth_data.py | 39 +++++++++++++++++++++++++++++ specreduce/utils/synth_data.py | 10 ++++++++ 2 files changed, 49 insertions(+) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index c00ed4da..2a7727c2 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -162,3 +162,42 @@ def test_add_skylines_matches_add_arcs(): sky = SynthImage(nx=300, ny=100, wcs=wcs, seed=1).add_skylines("OH_GMOS") arc = SynthImage(nx=300, ny=100, wcs=wcs, seed=1).add_arcs("OH_GMOS") assert np.array_equal(sky.to_array(), arc.to_array()) + + +def test_read_noise_changes_image_and_is_reproducible(): + base = SynthImage(nx=40, ny=40, seed=42).add_background(100) + noisy = base.add_read_noise(5) + a = noisy.to_array() + b = noisy.to_array() + assert np.array_equal(a, b) # same seed -> reproducible + assert not np.allclose(a, 100) # read noise actually applied + assert abs(a.std() - 5) < 1.0 # ~ sigma + + +def test_poisson_noise_reproducible_with_seed(): + img = SynthImage(nx=40, ny=40, seed=7).add_background(100).add_poisson_noise() + assert np.array_equal(img.to_array(), img.to_array()) + + +def test_different_seed_differs(): + a = SynthImage(nx=40, ny=40, seed=1).add_background(100).add_poisson_noise().to_array() + b = SynthImage(nx=40, ny=40, seed=2).add_background(100).add_poisson_noise().to_array() + assert not np.array_equal(a, b) + + +def test_unseeded_runs(): + arr = SynthImage(nx=20, ny=20).add_background(100).add_poisson_noise().to_array() + assert arr.shape == (20, 20) + + +def test_rdnoise_alias_equivalent(): + a = SynthImage(nx=20, ny=20, seed=3).add_background(50).add_read_noise(4).to_array() + b = SynthImage(nx=20, ny=20, seed=3).add_background(50).add_rdnoise(4).to_array() + assert np.array_equal(a, b) + + +def test_poisson_then_read_noise_both_applied(): + # both noise stages applied; result differs from poisson-only with same seed + poisson_only = SynthImage(nx=40, ny=40, seed=5).add_background(100).add_poisson_noise() + both = poisson_only.add_read_noise(5) + assert not np.array_equal(poisson_only.to_array(), both.to_array()) diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 3a123f82..92fde3b9 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -149,6 +149,16 @@ def add_skylines(self, linelists="OH_GMOS", **kwargs) -> "SynthImage": """Add night-sky airglow emission lines (OH lists), wrapping ``add_arcs``.""" return self.add_arcs(linelists, **kwargs) + def add_poisson_noise(self) -> "SynthImage": + """Apply Poisson noise to the rendered signal (requires photutils).""" + return self._clone(_poisson=True) + + def add_read_noise(self, sigma: float) -> "SynthImage": + """Add Gaussian read noise of standard deviation ``sigma`` (counts).""" + return self._clone(_read_noise=sigma) + + add_rdnoise = add_read_noise + def _resolve_wcs(self): has_arc = any(isinstance(layer, ArcLayer) for layer in self._layers) if self._wcs is not None: From c5835b1110aee296dd88abfdcbb7f1921f4eedf0 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 23:33:05 -0600 Subject: [PATCH 06/18] synth_data refactor: cover to_spectrum output --- specreduce/tests/test_synth_data.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index 2a7727c2..a5b44ed8 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -5,6 +5,8 @@ from astropy.nddata import CCDData from astropy.wcs import WCS +from specutils import Spectrum + from specreduce.utils.synth_data import SynthImage @@ -201,3 +203,20 @@ def test_poisson_then_read_noise_both_applied(): poisson_only = SynthImage(nx=40, ny=40, seed=5).add_background(100).add_poisson_noise() both = poisson_only.add_read_noise(5) assert not np.array_equal(poisson_only.to_array(), both.to_array()) + + +def test_to_spectrum_no_wcs(): + sp = SynthImage(nx=60, ny=20).add_background(5).to_spectrum() + assert isinstance(sp, Spectrum) + assert sp.flux.shape == (20, 60) + assert sp.flux.unit == u.count + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_to_spectrum_with_wcs(): + wcs = _linear_wcs(300) + sp = SynthImage(nx=300, ny=100, wcs=wcs).add_arcs(["HeI"]).to_spectrum() + assert isinstance(sp, Spectrum) + assert sp.flux.shape == (100, 300) + assert sp.spectral_axis is not None From 57b21fbfa36bf8261d70e328f0fbda756b92a991 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sun, 31 May 2026 23:39:45 -0600 Subject: [PATCH 07/18] refactor synth_data: migrate synth_data callers to SynthImage --- specreduce/tests/conftest.py | 21 ++++++++++++--------- specreduce/tests/test_line_matching.py | 21 ++++++++++++--------- specreduce/tests/test_tracing.py | 10 ++++++++-- 3 files changed, 32 insertions(+), 20 deletions(-) diff --git a/specreduce/tests/conftest.py b/specreduce/tests/conftest.py index ec0b0e04..63983f60 100644 --- a/specreduce/tests/conftest.py +++ b/specreduce/tests/conftest.py @@ -5,7 +5,7 @@ from astropy.wcs import WCS from specreduce.tilt_correction import TiltCorrection -from specreduce.utils.synth_data import make_2d_arc_image +from specreduce.utils.synth_data import SynthImage # Arc frame creation code taken from Tim Pickering's example notebook @@ -31,14 +31,17 @@ def mk_arc_frames(): arcs = [] for ll in (["HeI", "NeI", "XeI"], ["ArI"]): - arc = make_2d_arc_image( - nx=512, - ny=128, - linelists=ll, - wcs=blue_channel_wcs, - line_fwhm=3, - tilt_func=tilt_mod, - amplitude_scale=1e-2, + arc = ( + SynthImage(nx=512, ny=128, wcs=blue_channel_wcs) + .add_background(5) + .add_arcs( + linelists=ll, + line_fwhm=3, + tilt_func=tilt_mod, + amplitude_scale=1e-2, + ) + .add_poisson_noise() + .to_ccddata() ) arc.wcs = None arc.uncertainty = StdDevUncertainty(np.full_like(arc.data, 5)) diff --git a/specreduce/tests/test_line_matching.py b/specreduce/tests/test_line_matching.py index c70af135..a96d0675 100644 --- a/specreduce/tests/test_line_matching.py +++ b/specreduce/tests/test_line_matching.py @@ -12,7 +12,7 @@ from specreduce.extract import BoxcarExtract from specreduce.line_matching import match_lines_wcs, find_arc_lines from specreduce.tracing import FlatTrace -from specreduce.utils.synth_data import make_2d_arc_image +from specreduce.utils.synth_data import SynthImage @pytest.fixture @@ -53,14 +53,17 @@ def mk_test_data(): linear_wcs = WCS(header=linear_header) tilt_mod = models.Legendre1D(degree=2, c0=50, c1=0, c2=100) - match_im = make_2d_arc_image( - nx=1400, - ny=1024, - linelists=['HeI', 'NeI'], - wcs=linear_wcs, - line_fwhm=5, - tilt_func=tilt_mod, - amplitude_scale=5e-4 + match_im = ( + SynthImage(nx=1400, ny=1024, wcs=linear_wcs) + .add_background(5) + .add_arcs( + linelists=['HeI', 'NeI'], + line_fwhm=5, + tilt_func=tilt_mod, + amplitude_scale=5e-4, + ) + .add_poisson_noise() + .to_ccddata() ) arclist = load_pypeit_calibration_lines(['HeI', 'NeI'])['wavelength'] diff --git a/specreduce/tests/test_tracing.py b/specreduce/tests/test_tracing.py index 3b85a003..917d7332 100644 --- a/specreduce/tests/test_tracing.py +++ b/specreduce/tests/test_tracing.py @@ -3,10 +3,16 @@ from astropy.modeling import fitting, models from astropy.nddata import NDData import astropy.units as u -from specreduce.utils.synth_data import make_2d_trace_image +from specreduce.utils.synth_data import SynthImage from specreduce.tracing import Trace, FlatTrace, ArrayTrace, FitTrace -IM = make_2d_trace_image() +IM = ( + SynthImage(nx=3000, ny=1000) + .add_background(5) + .add_source(profile=models.Moffat1D(amplitude=10, alpha=0.1)) + .add_poisson_noise() + .to_ccddata() +) def mk_img(nrows=200, ncols=160, nan_slices=None, add_noise=True): From d83bf11f7cc2ee519d8cf3dfb18e7555b5efbe22 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 00:00:19 -0600 Subject: [PATCH 08/18] synth_data refactor: deprecate make_2d_* shims in favor of SynthImage --- specreduce/tests/test_synth_data.py | 42 ++ specreduce/utils/synth_data.py | 586 +++++++++------------------- 2 files changed, 223 insertions(+), 405 deletions(-) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index a5b44ed8..a657ae0f 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -3,11 +3,17 @@ from astropy import units as u from astropy.modeling import models from astropy.nddata import CCDData +from astropy.utils.exceptions import AstropyDeprecationWarning from astropy.wcs import WCS from specutils import Spectrum from specreduce.utils.synth_data import SynthImage +from specreduce.utils.synth_data import ( + make_2d_trace_image, + make_2d_arc_image, + make_2d_spec_image, +) def test_empty_image_is_zeros(): @@ -220,3 +226,39 @@ def test_to_spectrum_with_wcs(): assert isinstance(sp, Spectrum) assert sp.flux.shape == (100, 300) assert sp.spectral_axis is not None + + +def test_make_2d_trace_image_deprecated_and_matches(): + with pytest.warns(AstropyDeprecationWarning): + ccd = make_2d_trace_image(add_noise=False) + expected = ( + SynthImage(nx=3000, ny=1000) + .add_background(5) + .add_source(profile=models.Moffat1D(amplitude=10, alpha=0.1)) + .to_array() + ) + assert isinstance(ccd, CCDData) + assert np.array_equal(ccd.data, expected) + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_make_2d_arc_image_deprecated_and_matches(): + with pytest.warns(AstropyDeprecationWarning): + ccd = make_2d_arc_image(nx=300, ny=100, add_noise=False) + expected = ( + SynthImage(nx=300, ny=100, extent=(3500, 7000)) + .add_background(5) + .add_arcs(("HeI",)) + .to_array() + ) + assert np.array_equal(ccd.data, expected) + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_make_2d_spec_image_deprecated(): + with pytest.warns(AstropyDeprecationWarning): + ccd = make_2d_spec_image(nx=300, ny=100, add_noise=False) + assert isinstance(ccd, CCDData) + assert ccd.data.shape == (100, 300) diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 92fde3b9..a157e8c5 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -7,6 +7,7 @@ from astropy.modeling import models, Model from astropy.nddata import CCDData from astropy.stats import gaussian_fwhm_to_sigma +from astropy.utils.decorators import deprecated from astropy.wcs import WCS from specutils import Spectrum @@ -74,6 +75,111 @@ class SynthImage: seed Seed for the random number generator used by the noise layers. If ``None``, noise is non-deterministic. + + Examples + -------- + This is an example of modeling a spectrograph whose output is curved in the + cross-dispersion direction: + + .. plot:: + :include-source: + + import matplotlib.pyplot as plt + import numpy as np + from astropy.modeling import models + import astropy.units as u + from specreduce.utils.synth_data import SynthImage + + model_deg2 = models.Legendre1D(degree=2, c0=50, c1=0, c2=100) + im = ( + SynthImage(nx=3000, ny=1000) + .add_background(5) + .add_arcs(['HeI', 'ArI', 'ArII'], line_fwhm=3, tilt_func=model_deg2) + .add_poisson_noise() + .to_ccddata() + ) + fig = plt.figure(figsize=(10, 6)) + plt.imshow(im) + + The FITS WCS standard implements ideal world coordinate functions based on the physics + of simple dispersers. This is described in detail by Paper III, + https://www.aanda.org/articles/aa/pdf/2006/05/aa3818-05.pdf. This can be used to model a + non-linear dispersion relation based on the properties of a spectrograph. This example + recreates Figure 5 in that paper using a spectrograph with a 450 lines/mm volume phase + holographic grism. Standard gratings only use the first three ``PV`` terms: + + .. plot:: + :include-source: + + import numpy as np + import matplotlib.pyplot as plt + from astropy.wcs import WCS + import astropy.units as u + from specreduce.utils.synth_data import SynthImage + + non_linear_header = { + 'CTYPE1': 'AWAV-GRA', # Grating dispersion function with air wavelengths + 'CUNIT1': 'Angstrom', # Dispersion units + 'CRPIX1': 719.8, # Reference pixel [pix] + 'CRVAL1': 7245.2, # Reference value [Angstrom] + 'CDELT1': 2.956, # Linear dispersion [Angstrom/pix] + 'PV1_0': 4.5e5, # Grating density [1/m] + 'PV1_1': 1, # Diffraction order + 'PV1_2': 27.0, # Incident angle [deg] + 'PV1_3': 1.765, # Reference refraction + 'PV1_4': -1.077e6, # Refraction derivative [1/m] + 'CTYPE2': 'PIXEL', # Spatial detector coordinates + 'CUNIT2': 'pix', # Spatial units + 'CRPIX2': 1, # Reference pixel + 'CRVAL2': 0, # Reference value + 'CDELT2': 1 # Spatial units per pixel + } + + linear_header = { + 'CTYPE1': 'AWAV', # Grating dispersion function with air wavelengths + 'CUNIT1': 'Angstrom', # Dispersion units + 'CRPIX1': 719.8, # Reference pixel [pix] + 'CRVAL1': 7245.2, # Reference value [Angstrom] + 'CDELT1': 2.956, # Linear dispersion [Angstrom/pix] + 'CTYPE2': 'PIXEL', # Spatial detector coordinates + 'CUNIT2': 'pix', # Spatial units + 'CRPIX2': 1, # Reference pixel + 'CRVAL2': 0, # Reference value + 'CDELT2': 1 # Spatial units per pixel + } + + non_linear_wcs = WCS(non_linear_header) + linear_wcs = WCS(linear_header) + + # this re-creates Paper III, Figure 5 + pix_array = 200 + np.arange(1400) + nlin = non_linear_wcs.spectral.pixel_to_world(pix_array) + lin = linear_wcs.spectral.pixel_to_world(pix_array) + resid = (nlin - lin).to(u.Angstrom) + plt.plot(pix_array, resid) + plt.xlabel("Pixel") + plt.ylabel("Correction (Angstrom)") + plt.show() + + nlin_im = ( + SynthImage(nx=600, ny=512, wcs=non_linear_wcs) + .add_background(5) + .add_arcs(['HeI', 'NeI'], line_fwhm=3, wave_air=True) + .add_poisson_noise() + .to_ccddata() + ) + lin_im = ( + SynthImage(nx=600, ny=512, wcs=linear_wcs) + .add_background(5) + .add_arcs(['HeI', 'NeI'], line_fwhm=3, wave_air=True) + .add_poisson_noise() + .to_ccddata() + ) + + # subtracting the linear simulation from the non-linear one shows how the + # positions of lines diverge between the two cases + plt.imshow(nlin_im.data - lin_im.data) + plt.show() """ def __init__( @@ -301,76 +407,41 @@ def render(self, ctx: _RenderContext) -> np.ndarray: return z +@deprecated("1.10", alternative="SynthImage") def make_2d_trace_image( nx: int = 3000, ny: int = 1000, background: int | float = 5, trace_center: int | float | None = None, trace_order: int = 3, - trace_coeffs: None | dict[str, int | float] = None, - profile: Model = models.Moffat1D(amplitude=10, alpha=0.1), - add_noise: bool = True + trace_coeffs: dict | None = None, + profile: Model | None = None, + add_noise: bool = True, ) -> CCDData: - """ - Create synthetic 2D spectroscopic image with a single source. The spatial (y-axis) position - of the source along the dispersion (x-axis) direction is modeled using a Chebyshev polynomial. - The flux units are counts and the noise is modeled as Poisson. - - Parameters - ---------- - nx - Size of image in X axis which is assumed to be the dispersion axis - - ny - Size of image in Y axis which is assumed to be the spatial axis - - background - Level of constant background in counts - - trace_center - Zeropoint of the trace. If None, then use center of Y (spatial) axis. - - trace_order - Order of the Chebyshev polynomial used to model the source's trace - - trace_coeffs - Dict containing the Chebyshev polynomial coefficients to use in the trace model + """Deprecated. Use :class:`SynthImage` instead. - profile - Model to use for the source's spatial profile - - add_noise - If True, add Poisson noise to the image - - Returns - ------- - ccd_im - `~astropy.nddata.CCDData` instance containing synthetic 2D spectroscopic image + Equivalent to ``SynthImage(nx, ny).add_background(background) + .add_source(...)`` followed by ``.add_poisson_noise()`` when ``add_noise``, + then ``.to_ccddata()``. """ - if trace_coeffs is None: - trace_coeffs = {'c0': 0, 'c1': 50, 'c2': 100} - x = np.arange(nx) - y = np.arange(ny) - xx, yy = np.meshgrid(x, y) - - if trace_center is None: - trace_center = ny / 2 - - trace_mod = models.Chebyshev1D(degree=trace_order, **trace_coeffs) - trace = yy - trace_center + trace_mod(xx/nx) - z = background + profile(trace) - + if profile is None: + profile = models.Moffat1D(amplitude=10, alpha=0.1) + img = ( + SynthImage(nx=nx, ny=ny) + .add_background(background) + .add_source( + profile=profile, + trace_center=trace_center, + trace_order=trace_order, + trace_coeffs=trace_coeffs, + ) + ) if add_noise: - from photutils.datasets import apply_poisson_noise - trace_image = apply_poisson_noise(z) - else: - trace_image = z - - ccd_im = CCDData(trace_image, unit=u.count) - - return ccd_im + img = img.add_poisson_noise() + return img.to_ccddata() +@deprecated("1.10", alternative="SynthImage") def make_2d_arc_image( nx: int = 3000, ny: int = 1000, @@ -379,251 +450,32 @@ def make_2d_arc_image( wave_unit: u.Unit = u.Angstrom, wave_air: bool = False, background: int | float = 5, - line_fwhm: float = 5., - linelists: list[str] = ('HeI',), - amplitude_scale: float = 1., - tilt_func: Model = models.Legendre1D(degree=0), - add_noise: bool = True + line_fwhm: float = 5.0, + linelists=("HeI",), + amplitude_scale: float = 1.0, + tilt_func: Model | None = None, + add_noise: bool = True, ) -> CCDData: - """ - Create synthetic 2D spectroscopic image of reference emission lines, e.g. a calibration - arc lamp. Currently, linelists from ``pypeit`` are supported and are selected by string or - list of strings that is passed to `~specreduce.calibration_data.load_pypeit_calibration_lines`. - If a ``wcs`` is not provided, one is created using ``extent`` and ``wave_unit`` with - dispersion along the X axis. - - Parameters - ---------- - nx - Size of image in X axis which is assumed to be the dispersion axis - - ny - Size of image in Y axis which is assumed to be the spatial axis - - wcs - 2D WCS to apply to the image. Must have a spectral axis defined along with - appropriate spectral wavelength units. - - extent - If ``wcs`` is not provided, this defines the beginning and end wavelengths - of the dispersion axis. - - wave_unit - If ``wcs`` is not provided, this defines the wavelength units of the - dispersion axis. - - wave_air - If True, convert the vacuum wavelengths used by ``pypeit`` to air wavelengths. - - background - Level of constant background in counts - - line_fwhm - Gaussian FWHM of the emission lines in pixels - - linelists - Specification for linelists to load from ``pypeit`` - - amplitude_scale - Scale factor to apply to amplitudes provided in the linelists - - tilt_func - The tilt function to apply along the cross-dispersion axis to simulate - tilted or curved emission lines. - - add_noise - If True, add Poisson noise to the image; requires ``photutils`` to be installed. - - Returns - ------- - ccd_im - `~astropy.nddata.CCDData` instance containing synthetic 2D spectroscopic image - - Examples - -------- - This is an example of modeling a spectrograph whose output is curved in the - cross-dispersion direction: - - .. plot:: - :include-source: - - import matplotlib.pyplot as plt - import numpy as np - from astropy.modeling import models - import astropy.units as u - from specreduce.utils.synth_data import make_2d_arc_image - - model_deg2 = models.Legendre1D(degree=2, c0=50, c1=0, c2=100) - im = make_2d_arc_image( - linelists=['HeI', 'ArI', 'ArII'], - line_fwhm=3, - tilt_func=model_deg2 - ) - fig = plt.figure(figsize=(10, 6)) - plt.imshow(im) - - The FITS WCS standard implements ideal world coordinate functions based on the physics - of simple dispersers. This is described in detail by Paper III, - https://www.aanda.org/articles/aa/pdf/2006/05/aa3818-05.pdf. This can be used to model a - non-linear dispersion relation based on the properties of a spectrograph. This example - recreates Figure 5 in that paper using a spectrograph with a 450 lines/mm volume phase - holographic grism. Standard gratings only use the first three ``PV`` terms: - - .. plot:: - :include-source: - - import numpy as np - import matplotlib.pyplot as plt - from astropy.wcs import WCS - import astropy.units as u - from specreduce.utils.synth_data import make_2d_arc_image - - non_linear_header = { - 'CTYPE1': 'AWAV-GRA', # Grating dispersion function with air wavelengths - 'CUNIT1': 'Angstrom', # Dispersion units - 'CRPIX1': 719.8, # Reference pixel [pix] - 'CRVAL1': 7245.2, # Reference value [Angstrom] - 'CDELT1': 2.956, # Linear dispersion [Angstrom/pix] - 'PV1_0': 4.5e5, # Grating density [1/m] - 'PV1_1': 1, # Diffraction order - 'PV1_2': 27.0, # Incident angle [deg] - 'PV1_3': 1.765, # Reference refraction - 'PV1_4': -1.077e6, # Refraction derivative [1/m] - 'CTYPE2': 'PIXEL', # Spatial detector coordinates - 'CUNIT2': 'pix', # Spatial units - 'CRPIX2': 1, # Reference pixel - 'CRVAL2': 0, # Reference value - 'CDELT2': 1 # Spatial units per pixel - } - - linear_header = { - 'CTYPE1': 'AWAV', # Grating dispersion function with air wavelengths - 'CUNIT1': 'Angstrom', # Dispersion units - 'CRPIX1': 719.8, # Reference pixel [pix] - 'CRVAL1': 7245.2, # Reference value [Angstrom] - 'CDELT1': 2.956, # Linear dispersion [Angstrom/pix] - 'CTYPE2': 'PIXEL', # Spatial detector coordinates - 'CUNIT2': 'pix', # Spatial units - 'CRPIX2': 1, # Reference pixel - 'CRVAL2': 0, # Reference value - 'CDELT2': 1 # Spatial units per pixel - } - - non_linear_wcs = WCS(non_linear_header) - linear_wcs = WCS(linear_header) - - # this re-creates Paper III, Figure 5 - pix_array = 200 + np.arange(1400) - nlin = non_linear_wcs.spectral.pixel_to_world(pix_array) - lin = linear_wcs.spectral.pixel_to_world(pix_array) - resid = (nlin - lin).to(u.Angstrom) - plt.plot(pix_array, resid) - plt.xlabel("Pixel") - plt.ylabel("Correction (Angstrom)") - plt.show() - - nlin_im = make_2d_arc_image( - nx=600, - ny=512, - linelists=['HeI', 'NeI'], - line_fwhm=3, - wave_air=True, - wcs=non_linear_wcs + """Deprecated. Use :class:`SynthImage` with ``.add_arcs(...)`` instead.""" + if tilt_func is None: + tilt_func = models.Legendre1D(degree=0) + img = ( + SynthImage(nx=nx, ny=ny, wcs=wcs, extent=extent, wave_unit=wave_unit) + .add_background(background) + .add_arcs( + linelists=linelists, + line_fwhm=line_fwhm, + amplitude_scale=amplitude_scale, + wave_air=wave_air, + tilt_func=tilt_func, ) - lin_im = make_2d_arc_image( - nx=600, - ny=512, - linelists=['HeI', 'NeI'], - line_fwhm=3, - wave_air=True, - wcs=linear_wcs - ) - - # subtracting the linear simulation from the non-linear one shows how the - # positions of lines diverge between the two cases - plt.imshow(nlin_im.data - lin_im.data) - plt.show() - - """ - if wcs is None: - if extent is None: - raise ValueError("Must specify either a wavelength extent or a WCS.") - if len(extent) != 2: - raise ValueError("Wavelength extent must be of length 2.") - if u.get_physical_type(wave_unit) != 'length': - raise ValueError("Wavelength unit must be a length unit.") - wcs = WCS(naxis=2) - wcs.wcs.ctype[0] = 'WAVE' - wcs.wcs.ctype[1] = 'PIXEL' - wcs.wcs.cunit[0] = wave_unit - wcs.wcs.cunit[1] = u.pixel - wcs.wcs.crval[0] = extent[0] - wcs.wcs.cdelt[0] = (extent[1] - extent[0]) / nx - wcs.wcs.crval[1] = 0 - wcs.wcs.cdelt[1] = 1 - else: - if wcs.spectral.naxis != 1: - raise ValueError("Provided WCS must have a spectral axis.") - if wcs.naxis != 2: - raise ValueError("WCS must have NAXIS=2 for a 2D image.") - - x = np.arange(nx) - y = np.arange(ny) - xx, yy = np.meshgrid(x, y) - - is_spectral = [a['coordinate_type'] == "spectral" for a in wcs.get_axis_types()] - if is_spectral[0]: - disp_axis = 0 - else: - disp_axis = 1 - - if tilt_func is not None: - if not isinstance( - tilt_func, - (models.Legendre1D, models.Chebyshev1D, models.Polynomial1D, models.Hermite1D) - ): - raise ValueError( - "The only tilt functions currently supported are 1D polynomials " - "from astropy.models." - ) - - if disp_axis == 0: - xx = xx + tilt_func((yy - ny/2)/ny) - else: - yy = yy + tilt_func((xx - nx/2)/nx) - - z = background + np.zeros((ny, nx)) - - linelist = load_pypeit_calibration_lines(linelists, wave_air=wave_air) - - if linelist is not None: - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", message="No observer defined on WCS.*") - line_disp_positions = wcs.spectral.world_to_pixel(linelist['wavelength']) - - line_sigma = gaussian_fwhm_to_sigma * line_fwhm - for line_pos, ampl in zip(line_disp_positions, linelist['amplitude']): - line_mod = models.Gaussian1D( - amplitude=ampl * amplitude_scale, - mean=line_pos, - stddev=line_sigma - ) - if disp_axis == 0: - z += line_mod(xx) - else: - z += line_mod(yy) - + ) if add_noise: - from photutils.datasets import apply_poisson_noise - arc_image = apply_poisson_noise(z) - else: - arc_image = z - - ccd_im = CCDData(arc_image, unit=u.count, wcs=wcs) - - return ccd_im + img = img.add_poisson_noise() + return img.to_ccddata() +@deprecated("1.10", alternative="SynthImage") def make_2d_spec_image( nx: int = 3000, ny: int = 1000, @@ -632,115 +484,39 @@ def make_2d_spec_image( wave_unit: u.Unit = u.Angstrom, wave_air: bool = False, background: int | float = 5, - line_fwhm: float = 5., - linelists: list[str] = ('OH_GMOS',), - amplitude_scale: float = 1., - tilt_func: Model = models.Legendre1D(degree=0), + line_fwhm: float = 5.0, + linelists=("OH_GMOS",), + amplitude_scale: float = 1.0, + tilt_func: Model | None = None, trace_center: int | float | None = None, trace_order: int = 3, - trace_coeffs: None | dict[str, int | float] = None, - source_profile: Model = models.Moffat1D(amplitude=10, alpha=0.1), - add_noise: bool = True + trace_coeffs: dict | None = None, + source_profile: Model | None = None, + add_noise: bool = True, ) -> CCDData: - """ - Make a synthetic 2D spectrum image containing both emission lines and - a trace of a continuum source. - - Parameters - ---------- - nx - Number of pixels in the dispersion direction. - - ny - Number of pixels in the spatial direction. - - wcs - 2D WCS to apply to the image. Must have a spectral axis defined along with - appropriate spectral wavelength units. - - extent - If ``wcs`` is not provided, this defines the beginning and end wavelengths - of the dispersion axis. - - wave_unit - If ``wcs`` is not provided, this defines the wavelength units of the - dispersion axis. - - wave_air - If True, convert the vacuum wavelengths used by ``pypeit`` to air wavelengths. - - background - Constant background level in counts. - - line_fwhm - Gaussian FWHM of the emission lines in pixels - - linelists - Specification for linelists to load from ``pypeit`` - - amplitude_scale - Scale factor to apply to amplitudes provided in the linelists - - tilt_func - The tilt function to apply along the cross-dispersion axis to simulate - tilted or curved emission lines. - - trace_center - Zeropoint of the trace. If None, then use center of Y (spatial) axis. - - trace_order - Order of the Chebyshev polynomial used to model the source's trace - - trace_coeffs - Dict containing the Chebyshev polynomial coefficients to use in the trace model - - source_profile - Model to use for the source's spatial profile - - add_noise - If True, add Poisson noise to the image; requires ``photutils`` to be installed. - - - Returns - ------- - ccd_im - `~astropy.nddata.CCDData` instance containing synthetic 2D spectroscopic image - """ - if trace_coeffs is None: - trace_coeffs = {'c0': 0, 'c1': 50, 'c2': 100} - - arc_image = make_2d_arc_image( - nx=nx, - ny=ny, - wcs=wcs, - extent=extent, - wave_unit=wave_unit, - wave_air=wave_air, - background=0, - line_fwhm=line_fwhm, - linelists=linelists, - amplitude_scale=amplitude_scale, - tilt_func=tilt_func, - add_noise=False - ) - - trace_image = make_2d_trace_image( - nx=nx, - ny=ny, - background=0, - trace_center=trace_center, - trace_order=trace_order, - trace_coeffs=trace_coeffs, - profile=source_profile, - add_noise=False + """Deprecated. Use :class:`SynthImage` with ``.add_arcs(...)`` and + ``.add_source(...)`` instead.""" + if tilt_func is None: + tilt_func = models.Legendre1D(degree=0) + if source_profile is None: + source_profile = models.Moffat1D(amplitude=10, alpha=0.1) + img = ( + SynthImage(nx=nx, ny=ny, wcs=wcs, extent=extent, wave_unit=wave_unit) + .add_background(background) + .add_arcs( + linelists=linelists, + line_fwhm=line_fwhm, + amplitude_scale=amplitude_scale, + wave_air=wave_air, + tilt_func=tilt_func, + ) + .add_source( + profile=source_profile, + trace_center=trace_center, + trace_order=trace_order, + trace_coeffs=trace_coeffs, + ) ) - - spec_image = arc_image.data + trace_image.data + background - if add_noise: - from photutils.datasets import apply_poisson_noise - spec_image = apply_poisson_noise(spec_image) - - ccd_im = CCDData(spec_image, unit=u.count, wcs=arc_image.wcs) - - return ccd_im + img = img.add_poisson_noise() + return img.to_ccddata() From d27ca052f17a152df49d8bf79232239bcf0f884e Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 00:03:38 -0600 Subject: [PATCH 09/18] migrate quickstart synth_data calls to SynthImage --- docs/getting_started/quickstart.py | 50 +++++++++++++++--------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/docs/getting_started/quickstart.py b/docs/getting_started/quickstart.py index d3adbd62..8ce1a150 100644 --- a/docs/getting_started/quickstart.py +++ b/docs/getting_started/quickstart.py @@ -9,7 +9,7 @@ import astropy.units as u from photutils.datasets import apply_poisson_noise -from specreduce.utils.synth_data import make_2d_arc_image, make_2d_trace_image +from specreduce.utils.synth_data import SynthImage def make_2d_spec_image( @@ -96,29 +96,26 @@ def make_2d_spec_image( if trace_coeffs is None: trace_coeffs = {"c0": 0, "c1": 50, "c2": 100} - arc_image = make_2d_arc_image( - nx=nx, - ny=ny, - wcs=wcs, - extent=extent, - wave_unit=wave_unit, - wave_air=wave_air, - background=0, - line_fwhm=line_fwhm, - linelists=linelists, - tilt_func=tilt_func, - add_noise=False, + arc_image = ( + SynthImage(nx=nx, ny=ny, wcs=wcs, extent=extent, wave_unit=wave_unit) + .add_arcs( + linelists=linelists, + line_fwhm=line_fwhm, + wave_air=wave_air, + tilt_func=tilt_func, + ) + .to_ccddata() ) - trace_image = make_2d_trace_image( - nx=nx, - ny=ny, - background=0, - trace_center=trace_center, - trace_order=trace_order, - trace_coeffs=trace_coeffs, - profile=source_profile, - add_noise=False, + trace_image = ( + SynthImage(nx=nx, ny=ny) + .add_source( + profile=source_profile, + trace_center=trace_center, + trace_order=trace_order, + trace_coeffs=trace_coeffs, + ) + .to_ccddata() ) wl = wcs.spectral.pixel_to_world(np.arange(nx)).to(u.nm).value @@ -176,11 +173,14 @@ def make_science_and_arcs(ndisp: int = 1000, ncross: int = 300): source_profile=models.Moffat1D(amplitude=1, alpha=0.3), ) - arcargs = dict(wcs=wcs, line_fwhm=3, background=0, add_noise=False) arcs = [] for linelist in ["HeI", "NeI"]: - arc = make_2d_arc_image(ndisp, ncross, linelists=[linelist], **arcargs) - arc.data = apply_poisson_noise(200*(arc.data / arc.data.max()) + 10) - 10 + arc = ( + SynthImage(nx=ndisp, ny=ncross, wcs=wcs) + .add_arcs(linelists=[linelist], line_fwhm=3) + .to_ccddata() + ) + arc.data = apply_poisson_noise(200 * (arc.data / arc.data.max()) + 10) - 10 arcs.append(arc) return science, arcs From 38d064490b66670e1dde3387617bab729dadd385 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 00:20:23 -0600 Subject: [PATCH 10/18] synth_data refactor: update CLAUDE.md, docs, and tighten deprecation test --- CLAUDE.md | 22 +++++++++++++++------- docs/getting_started/quickstart.ipynb | 13 ++----------- docs/getting_started/quickstart.py | 2 +- specreduce/tests/test_synth_data.py | 10 +++++++++- 4 files changed, 27 insertions(+), 20 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 39d2eb82..9f2d2323 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -52,7 +52,7 @@ The reduction workflow follows a modular pipeline: **Trace → Background → Ex ### Core Modules -- **core.py**: Base class `SpecreduceOperation` and `_ImageParser` for image format coercion. Handles Spectrum1D, CCDData, NDData, Quantity, and ndarray inputs. Defines `MaskingOption` enum and masking strategies. +- **core.py**: Base class `SpecreduceOperation` and the `parse_image()` function for image format coercion (replaced the former `_ImageParser` class in v1.9). Handles Spectrum (specutils v2), CCDData, NDData, Quantity, and ndarray inputs. `MaskingOption` is a `typing.Literal` of the seven supported masking strategies. - **tracing.py**: Trace determination classes - `Trace`: Base trace (center of image) @@ -68,27 +68,35 @@ The reduction workflow follows a modular pipeline: **Trace → Background → Ex - **wavecal1d.py**: Current wavelength calibration implementation using `WavelengthCalibration1D`. Supports automated line matching, template matching, and produces GWCS-based WCS objects. +- **wavesol1d.py**: `WavelengthSolution1D` — manages the pixel↔wavelength polynomial mapping underlying the calibration. + +- **line_matching.py**: Helpers for matching detected lines to reference line lists. + +- **tilt_correction.py**: `TiltCorrection` for correcting 2D spectral tilt/curvature. + +- **tilt_solution.py**: `TiltSolution` encapsulating the polynomial tilt transformation. + +- **table_utils.py**: Shared helpers for the QTable structures used across modules. + - **wavelength_calibration.py**: Legacy wavelength calibration (deprecated in v1.7.0, removal in v2.0) - **fluxcal.py**: `FluxCalibration` class for flux calibration with magnitude-to-flux conversion and airmass extinction correction - **calibration_data.py**: Spectrophotometric standards and line lists -- **utils/synth_data.py**: Synthetic spectroscopic data generation for testing - -- **compat.py**: Compatibility layer for specutils v1.x and v2.x +- **utils/synth_data.py**: Synthetic spectroscopic data generation for testing. `SynthImage` is an immutable, chainable builder (`add_background`/`add_source`/`add_arcs`/`add_skylines`/`add_poisson_noise`/`add_read_noise`, then `to_array`/`to_ccddata`/`to_spectrum`). The old `make_2d_trace_image`/`make_2d_arc_image`/`make_2d_spec_image` functions are deprecated shims around it. ### Key Design Patterns -1. **Image Format Flexibility**: All operations accept multiple input formats via `_ImageParser` -2. **Masking Strategies**: Seven masking options defined in `MaskingOption` enum +1. **Image Format Flexibility**: All operations accept multiple input formats via `parse_image()` +2. **Masking Strategies**: Seven masking options in the `MaskingOption` Literal (`apply`, `ignore`, `propagate`, `zero_fill`, `nan_fill`, `apply_mask_only`, `apply_nan_only`) 3. **Astropy Integration**: Uses Astropy models, units, and conventions throughout 4. **GWCS Support**: Wavelength calibration produces proper WCS objects ## Dependencies - Python ≥3.11 -- Core: numpy≥1.24, astropy≥5.3, scipy≥1.10, specutils≥1.9.1, matplotlib≥3.10, gwcs +- Core: numpy≥1.24, astropy≥6.0, scipy≥1.14, specutils≥2.0, matplotlib≥3.10, gwcs - Optional: photutils≥1.0 (stellar profile fitting), synphot (synthetic photometry) ## Code Style diff --git a/docs/getting_started/quickstart.ipynb b/docs/getting_started/quickstart.ipynb index b1d10a43..8f2819ae 100644 --- a/docs/getting_started/quickstart.ipynb +++ b/docs/getting_started/quickstart.ipynb @@ -53,16 +53,7 @@ "cell_type": "markdown", "id": "1ba4ac1f-0742-4a0b-81dc-6197f82a9c17", "metadata": {}, - "source": [ - "## Data preparation\n", - "\n", - "First, we read in our science and arc spectra. We use synthetic data\n", - "created by the `quickstart.make_science_and_arcs` utility function, which itself calls\n", - "`specreduce.utils.synth_data.make_2d_arc_image` and `specreduce.utils.synth_data.make_2d_spec_image`.\n", - "The science frame (`sci`) and the two arc frames (`arc_he` and `arc_ne`)\n", - "are returned as `astropy.nddata.CCDData` objects with `astropy.nddata.StdDevUncertainty`\n", - "uncertainties.\n" - ] + "source": "## Data preparation\n\nFirst, we read in our science and arc spectra. We use synthetic data\ncreated by the `quickstart.make_science_and_arcs` utility function, which itself uses\n`specreduce.utils.synth_data.SynthImage` to build the science and arc frames.\nThe science frame (`sci`) and the two arc frames (`arc_he` and `arc_ne`)\nare returned as `astropy.nddata.CCDData` objects with `astropy.nddata.StdDevUncertainty`\nuncertainties." }, { "cell_type": "code", @@ -751,4 +742,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/docs/getting_started/quickstart.py b/docs/getting_started/quickstart.py index 8ce1a150..8ab3933b 100644 --- a/docs/getting_started/quickstart.py +++ b/docs/getting_started/quickstart.py @@ -21,7 +21,7 @@ def make_2d_spec_image( wave_air: bool = False, background: int | float = 5, line_fwhm: float = 5.0, - linelists: list[str] = ("OH_GMOS"), + linelists: list[str] = ("OH_GMOS",), airglow_amplitude: float = 1.0, spectrum_amplitude: float = 1.0, tilt_func: Model = models.Legendre1D(degree=0), diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index a657ae0f..adcb3be1 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -257,8 +257,16 @@ def test_make_2d_arc_image_deprecated_and_matches(): @pytest.mark.remote_data @pytest.mark.filterwarnings("ignore:No observer defined on WCS") -def test_make_2d_spec_image_deprecated(): +def test_make_2d_spec_image_deprecated_and_matches(): with pytest.warns(AstropyDeprecationWarning): ccd = make_2d_spec_image(nx=300, ny=100, add_noise=False) + expected = ( + SynthImage(nx=300, ny=100, extent=(6500, 9500)) + .add_background(5) + .add_arcs(("OH_GMOS",)) + .add_source(profile=models.Moffat1D(amplitude=10, alpha=0.1)) + .to_array() + ) assert isinstance(ccd, CCDData) assert ccd.data.shape == (100, 300) + assert np.array_equal(ccd.data, expected) From 173cc002a0ca391d29860d0ca738911fc825c780 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 00:41:05 -0600 Subject: [PATCH 11/18] add dedicated synth_data guide; move SynthImage examples out of docstring --- docs/getting_started/science_spectrum.fits | Bin 20160 -> 20160 bytes docs/user_guide.rst | 3 +- specreduce/utils/synth_data.py | 100 ++------------------- 3 files changed, 9 insertions(+), 94 deletions(-) diff --git a/docs/getting_started/science_spectrum.fits b/docs/getting_started/science_spectrum.fits index 7eaa7c872d8a9a84c270cea4d22c4ad7fbe36204..d8bcf277448b6fa532a5d2fbd06c0675c57b88e3 100644 GIT binary patch literal 20160 zcmeIYXHb+|w=OD*C?CFr1wfwugEIZuSV-{5Q{&pJ)tcCd1|5$TK4H zJVHV|!WA}qgnB51h6j87y_Em>)y~9);Y9xZTYOqXUSL>&r&owVkoW(R>>J<}xZO7> zaQlA<`Im3zF3#q5X8%T>^S?a!U;O%~u>az-|2+1O2L92&KN|Q)1OI5?9}WDYfqyjc zf4>I)C?}m5_I4{2RA>x4XA_!@xyrwsN>jnpH_#&_Tp`rM)89+MYggz`i;A~@*seb${`c_xao!Pm|8wUZk+&%*EHG3d zFkrFve*@v4{ObN;&njVoo4i6oJ$wU0!~g#a{FC1idH-|g9g+7xa$eG(+nLy!|8ZV3 z6K9kEHoyP$=Oa${&)U(V-Q__HNGjB27dEkK>pfb(0(X4#5Uz4x;pH;|66Ma{MYw}>^=drsOH3cPZvbj ztM&9u%>n;g^wfR!%V4%7?U${X1fm;X2U|SQf^f&U*5muDVOk>lG|oXA(al+r{w|&n z(zDkp+|_{D^=*bhv*seY)m6$_(*eSUoC~EA8(?<#-0JryJ|MbnrcR$(D})gjUeopF z!?gI5e$fvPM7OuJ(b^gzbe)pzr7VP5onnqkIhKFd*0K#cf}{08PP2ziRR^P>7;+8w$>S$!K~Bs0pw!{s6TyvZa% zEEr}7S?XgvE+R%Iz%|QlI)snqMx}RK!AvIPw8GqC#K@Zp7HiIcaGki$sd3kc{N8V> z<8~rON!Cp%o(;j{^NVfW=fZ5-tfs}bl890NQMK1W1A=?&Q|C%ugel|Y-bE!@h|zgi zuJ>9G{Nm2kXQ`WD%6c+wlE{3-EV*EX_Vo~$@hhF~<-#=nUJg6=6=IMPb~sQS{Ee9O zZQ(i41}t3qUTY%8a_=(w8xinV&OI_HGY_V-$L`m;JcJlKpRDClX5hu7mguUFCidwf zT)l{em{pD$q4Q(FyJwZW`iC=2jZU3778QyZH|Tcf3czDN?l$|d2BuE z|6>s);C{Gx_jF8?NsF&8yvvXEhfwdnXF7a*l!PGa=|;u^edb z>-}63wjd^9e6R4aJ9s|TX6jD%puO6%y2@}oVz}eKrrd4?cMfgaB0pcyaw{zN`34{+ zbzbHQlQZ*W)yNd>of@YC_D(`FpZ1i~@I5#oQw)mqDWk<%JCOBc@pETK2(8P@apEf;0U< zo3rzB;CI4r8A_i|)&bRQo4MM?AGA18{SJ-Ih^caNWIod%@;|j*x~~P=rWFU(iX#wH zx8>k(!ya%?y-k|>q7}55j87lNO+-xNKEcPcG*CfF<4ea{gJ$e{SGsixVwyAcXBeLY z<+Eg3T=YTEL|coNO(S;FTAmcIa~RapV~X0>UxFt6^5OD2BB$*MwanZT+=c12MG}sn zEwoLWBQ=1SjscV9f+TPy`AQ~9Cy9J9mo36Mi0M_x=kniy3MqQ{X8tiEUvG@hvZsjo zLBFayuN>6u7u97szA$+%WvVOj4l%>qlXViBLDepLTU*`-lQv#o-GLK`9h*Hv=y?mA z>w^=oMel$~(fIu*o<}2g(mjV&6;nX%R<_kmv4Y9V)ReU^k05rMsA>G*3~)<|wlj)~ z2~LTthxM-@cFvNXg`c9p)k#<|;aD_i^2U})OT-bY7~*JNuLCY)<*J6slAwu|C4Z0b zL9FWO;I~sViG9ZBOx}7MG`n^4ONtJP^evFRfCJL{9}e}5qM-*3BmuoJOH>TPw7 zR^Zbobf3FM>^AcFx}`xh#I6i7nB~?B!MvEn)(h)EyEx;~#P>H5Yf;LOlJbOrT63#& z)@hh3^jRxS??tTb@SM``F%W2HEu5Ee7^dQ#=8Bt_A$AqxjMNVv1P5OkN@~o7sk6lz zgF=GS?SQD}W@89tLv}7$kqT4x)Sqj9iXe8=2icjuT@XU`sU_6}(>IBy+t-^TcI)yV zEf=aFJk>PE*OUvh=^e}5?o=Z-B;~*!-Md6i_2U;zBbbK84l(W#oZ(%+C9JI=EDhc| zu=hSpZ(VcCh4XL zu<~f~4hZ9Fw={1czpt5F$(S{S*yAgox6Y=>xmu+pE#qP8d%Wt#SPjIUKGopPTMS`v z#-=?H5->ecGw^6O1F@yDiL~J?2reWmm#o?j+S3IE=V``>t>8-w-o!x=xy^m8*fE$+ z9QZM-stB>yL<85(Q3n5u)sL>;kD#qN`(-fgAz~YLXBK@-0>69x&=sv`p#3~?(f;07 z#NPjrd1Q!!Ks~B(W!YBHK1*rdzUYP6*03OBpM~JB68GJcR0Y}!TFR8NG{m-jm&weI z0H4>KVJLG8G-}$UmeRpoweI!WZ=+PD20U&idn^_SRg{GT8u=6OKHK$4d}5uL5b%z(*}H)Um8t{_g{qh5+i0_Ta@q}g0Cn4Hb~EmA6qIHkc^;hWxrJNwnA zODn=)Vy_w#lhuT{#ThnU@`pfeNRRI}UkDQx-6uBA8gcqdHhwIY1eNBd^<~E`m<;GF zG|y%tZuzrS3ZbV!1-k8$JCFktmtJc7M zuId5(t%+j|)4n3Ema!n2u>gFzT|FG5>(GC^=+@y#PsFwS@pCfxw;eS?W*5O=#q7tVX=xFFg!6bPhQLL;u&9 zIi<^ABL4UYPU`c?@nK8G!C>Og;fq`*;wwjRQYmYBn*w<)Tu7KbmdeV02w;#V-sGW_( zh7p|H_9=;7Kd-@HRL7C6ZM8^zJ%W>0C(5(GWCsHs_sIg<8zlZ3!6{f4k!rGr4}(gp zg9xrblH>?Zp=^KNp{Ox1bf;Z+KY9^K$|E?FU4L|$6wQFqT1}6(^!G?wF@jT=yNJKc zT@^-t{hJ(hP9n*91gFsayVA*Ss94B;Ne&a&HhU2Wgpi8?kuRSnc< z#f(SSufTp?`nuTD(O|jeJkZSOA@yOclLbOVb|x<#FKP{o zm*=}~xEp}AzEfwxVVlGv*ssakTC?pq zSf09qhT?uuil6S+4tIbqen+HAdKFlk*HnKYbqMt^f1}&rJ6J?6NPMYP3YK@|r&$&I zz`NR;-fLb2Tc)aMHg7CgTaL=3)(rg7;TiMg=YuIqo!5w+50>9GTb~Y!r19N-ojB z+NRa|tc3}#wrYjIMhx~D`?9YLya8+bstIFWN`n{M9&K?Z3RaO7ukBJ~zzPdrqW1O? z_+xiS2Q0b`i{!J;nq9NN+IeV;!^Z>QE1z?^Qgj8@lh0SrIGhGnc+LJ@y>8$&tn3@d zPXSY0H+G#}7g&++>W_R^0AFcg-B0aK*oc_VlPC2cD@y!Z{lIGorcTT{W=?pSKIi;t zheKfP(^$|WavXx~lGYQvpTde;dy{E31FRTF=jr3NK(IPigEBAxbD6U;uiF!>xE-;} zCoh5^H#o`G>E?Dlfp&8xv1}mwmjXm=$1STcD zGDqiuak5Ki`yAmwrDF_BYi% zs+@@(9$t7+(vb;4fX$7}AHTuO)D4$9qz%?ln+`*lX%H9-oj>lK<)vA~Z#oF!t%eCd9=?XXk}6Xye=}IeFNU4+bwl#{ zpXX;R(uRZinCY79<-y8-p4S|t0O9NEbMBX4!FsEFL|-T2rBI}~KlC+(J6By%%UuAv zw76f0iY{2EmB-ABI1Zu3ut{QrJ8ZW;<0Kgp|2<=+L*KUr$vT=7a%~^O!QjBSGsfG& zDhc!mN>~KpM(K4&MH*mZXP`4rV>Vdlk_G9UUInB9%+ z$y<9w9nS>8R$Mgks6;D+OTns`Z>vLi5`Iyt&00(Fn@k*cr#BR=D$Cuj`Bg|(O6^iS9SLg_!@S0d zyI@_}TDY$$2Fa(SQhRoC;hM8@2q_6X(syJj7-9<}){lU1etHNmqk)?Ywx? zrST9p)tly9d4pcEVT`g|JXrNljJ{RhLh=d^v1c4ZSOqC&1PqY(-Wc2>TgOAPXWF}x zE*czYnzW5W#9tZ|_M0{^kzDlr)PSNVnC32x@}hPcW)Ym?iphr){B7N@xa%) zlEk+=TfPY&J%KRA=A=4&zz{d(%oE7qZa?;JTRj$qX>v1HT)KoT;{(?D~v-nt@Em+TJl99jsAZ%{H`qQ#-VDbjv z6fXknh4;#w--P#!yI+*IfVDfE3KR5IbbC3`plkA^x3&`EqC4(2&1l7-O`|77u!^@ zf+Y@C_a@eP71HP4cdTAf`5BIO?gu)0MuYV!vHa;`3kcP}W%&>1!2VNy!QSOzV0|v^ z7_CcqU-Ie4z0Q6x&bFjAoFnqT-jz``YJ_m!$KsyfX|OqAktb7t61Ge-z&G`r?Xcf!8M`sH)CBVhF}&pT{I{QYOI-^nMEpx;fd z?PQz;YtX&9&YtkDrSn-)dorw43>vqHh=Da6Ke}rb@pm`(`B+Kh&!;goln;VE>ZH!p zb%gg_^){VR)4{x_ysRv;0_@SZJ(haRgHXBERBF?6(1-hjSMS*ew&-Vphxc2ehXsp; z(u+V(Jyuz=-UaM&)2oyG^C5^Rd{%Im#KW=A6>lc(0b6`Y&zTTE@T2@D%E;~qGyO?b z{_$9_CvK2?v`Yg5%~R2~-#@?#@*bC!_kuk+#%5^GAM=NJ>+QmIpr3Fs)c&yzY>DH$ z7sQ?iKYPz0Z*ef}??+5KM*2p!)U85m)(-G5`5cXl@Q35%TIICPcCe@Rw1)7Pg1_Re zWr{i-wr|QJBd&LWJ!7hPM(QB=nxV8@WpOYD>IAZiJz&o?G`e!|3V6d+yFc}>g*`3S zV_6-MKYN`|`_UM1!}rk(c%E>G-ZVJiw+n39z5Ayew*Xg6ym{SdN7(e7{vs(Q4Yqu4 zz5Xe2a3kinsH@e&GH24&>lwemR=EDnt)vmW$eL?Z*?U;^PS@{KC3;ruS`b?%fIxT8 zqT&iUFv~`pF8GiEwvwb%L8T-3`<2I+KbC^sn6pK?Yc#-C(Tllvbq4qn)ge18%V2q^ z*|%bxE!e6qB|mRG0q@SeLA??h7+yv{GEU@!t-iZ$Zc`R`--UCQ_r3>{X}fZ%TMXEX zvn4GaYyj8LU1g)E8|(&x-zL$Gz}Bo=8PKW-?yPxD#cW+zK26`1dGQ+9+8+W_U%mx5 z_TVB7D<9C+KV=;z@t3VP`OxKe`QWW}HD`S$^_%=8X4{j&RzV7q`SD|U&*WofS zJ$2Jdlv}_yS~V#CSp&R1gFAXG7Q$-gkNC+iB*0!8p|17)GkBboV+V{~LC?OP`ZzfP z?B$2ot{FTJZqBygd`(~2TLkrLccz1FT+7-!dKWlfnKhPqEI5eY@OEvy2sYlAA0NMz z=wGwcWdAo<_xpT^%~k^2bVB>>$;04gh-P|yVuBt#snzw;eXz|nXY@_ILgr~=CgOry zBu}+SStJoXTRK_HnH5Lg`{nmu$va?n%GNl|9R)UHXBcg+1*nt0Z}f-4U_<#wiK?9k z+d4DPS7`#cBFoC>DG}aJ>GVX*CH7!j-7HXRBy$Y4U$I)Lu-4l@7#YO_+y3q7GA#kP zHZ}Gad6&V=$od%jO95=h@!Bs896?>q*qYt60Cw9m-+Xo<_V2vNW9;%7prSPC^(zWM z-+nu4LB=MqS2G2ROj!UVPB;~LDXWf z*B`khLGdBB!_qyQD!_!4z z3QCmKe%%Ax%eHleADIJxES>iJk`(AeqIHXoO#s^`M0|757f`966ld=-ho$m`+x5<# zVEd*Ual#1ia|)=ga8KC!)@d+~e*$~!C7-jAyFjH+k_$>Absf`f-l14+umfJCw?q@( zjdK>$gNk777j2RuGY;(F(e-AyWcEhYw= z^cR+L*Mh+gwO(MCOnBe3Cdb};EgVD(`<5Lr1A9l1Q`iAIsQW(m_FsAmTh{ndN}sL3 z-jx!Qb(q{+(ea4!?H1_zO)9I^gTRg`E4iB61ir<3!>oWbST71!uPdwud-u~?&Kmz!<`lhw59neA_|_l}a3I^zuXx82k1pW4E9<&Jje&BI_vFIZ_*Is=>sjW2Uc z>R~zB!{+w3Utq`512yOX(n$dryfp0JOMoF_I!aI(%1JZL1mZ@TCzw%8N&Ig=Z%r<;S#X;B-v zM|d~V8|_-_0Lu;iwK8t2!R8LFReW3kULn0N=}-&ob{2#e2>rko%wsX0`;)m`{kEH| zIxv^`9X*~#>?GN|Jot?kI3I)AebpLZyxp(6&*&Z4slM&|JHCLEx!&1UAq*B@3qF21 z5)O7cZ$?$m1#lu|g8Dbe!TLR~&fzWx>;q@b-}MoH&wHV|`)3F2+peb#iKT#j@KM;L z{^j6Q-Z!gRwiMP=Pi&mLtq$z0{yg2^!=U=?_Pp-+4vWP3M^&|-4FyD}3%F}Py~%N1a;t7XlU6C|#lIK%%T+XU)V@w5SDZCJjMjQW~L{Oe?M zwY*|7r~_P!jJb-iA06D8*P;b>;g25kh2);-k;Avg9n)afQ#h_fMi=ZNSvh}=nc%uz zm#`=!{mylI6+87I*u}6((Ix(#=VJDa`xKT&eo>v;qrpDAIr7qy!{A6dE-ej~fi1OR zJ*$QAbv~(}ZKXT0lXTb3>nDLJu~U6;feF|bPPL9VUjSZk0aKA*4MuUP%lwy;V3*$$ z*R&$MXD{77*}@E#lG0<|*AYCIz8S4%65baNpOXt43!8b544vEiz^J{mEZd#z87gaOL^I?^I>$8`B7}zx)bqyY0z%}R+eRp*?ELN+? ze?3Ct($)CBFFqGY-pM;)e`5+5_f%B_15LrMD_AfqfaHOW(2XinEMa#-@F=&CNIM8zr8x3GX+>C;D%g3p?|k zrhUWZVBe2jIYuxEoQH8ft*QaAT-fn~_e=!r2PXp6(we~GK8Q5z7=* zaL6Gu85}8-Z+G)|!gj*5T9v2CU_bfPxGUQkoZ}mKL3cO9uG6`~>Qma4J3LefjdgSE5;*6bk%C6gBzjSl$K6@D4a{k#ib01hI zJpGaW#UJdpeXMEao*)%|rJUdp>_fFH{eB$<`%PZ?vZ@8dKHn@=(a(idFw^4!dkNO(W%miy>%DQt3|*(N?EajnNSceu3=T;=Cg_s5YTx0iZ-vxd2hWXc1DX_ofjJA0n0q%>X9cTPKV1K7b>*En` zu>0z@x0CrFf6Vk(S|>>#B&HHQPlm+F?;qU{ei;Qp(yN{#av4c)UslsLR0Q^~Dg2rr zwGcS08d#&N2+L~<0|p*fz#h=4>KIA@kN2s6*UEUe!&UEm%h9M{3S{R^sYf4#8}DPj(&2FOg1(qG1K#%*}YUZ({@mTXXbNuDv9A ze|OW-?Qg)DS0!=6%NT;3p}Be|?t{56$UkbD7dZ1@8#iqwyzlvP$>fp&tdAL8%bfTP z9A&Y9Z^73f*r;!HlWhY_-z{q-`cHzhQ0CvAA%51b`7Go15| zLbzP^eHo_GYiS0;%Pt^wezNDa*?Rz|XQ-zAf+ zi{Rj}GUheo5I8H#b8nV*K$w=dsMqQjtUVt+zvFiU9Fu4FdoQ0NbDODrz33$*f1R4K zb~qCpGm+6UHRPU%ZW+>A?+=@(n-?oc-OsU5(l)-X1K~T(mga_oa8UgEVrj%maOezo zuf{$I^|pw2$xMZ9?bylBV;aG+3gA=s%OOm(lx~=L6?C`Kv5lv1gJUDCI`i`ARVtg%l>z3q@bMyDo56AQYdJJ91wx^UL_z-sIJg{kT~){eXB~fh-Ke`{ z?!^1xGiN?*3|{_>SWfJ7{W-%fu~f1@^Mq@t@dFIHymaVtqW_JLeWp%Y4dM8^@$Hur zVV6<7*7xoha6AUmmr4`fr z!SOcjOOhkJ-}>Hdy(SM<<-CTs{!78xvSt36`JPCQymNx_`Xub=i7)y2OtN=HIX+S) zygz;>BCV6{9=8;*)As$0v#mJ*mx&$;ygHtrhU}& zB>D_9G0re2_v3-a2R9b}gu~#Db5=*EgR|2s;EHuQnJa1^UAZ+MHl^(^^H%qR6VA$N zci4^O`Ru0S4|Kse%wPB^dk?AqPB%_jLwKKDX{lBD3wCRYVgp(2;Ox0Kpuc`J*{6y) zUYaiohcGJzy{WC>?ES9h=1F+BJ@`ve%pbN(C1o%DehW^_Y}eQ=g!k3Ks%Lx8fWAhQ z{nKd;IB~{F1%YHwV8ZWdUGXep|1T${9CQRH!L#gM=oIkpo!s<3kw*F#d%xP!?cgLO zy#E<~AA+_gYxJDbVc-7nR*(YmJ9goWxqDN=uRq28mOlaXYwZiA8iK%~?wDJ|yFf6j z(ah!mAB@O3YIc6lz~Oxj4d4*oPo1@A?_Hs7sx7y7FC(T{tb=GDGCo-0&E+BKH_?lNoht`9$ zKTcabM;W}{Q(lv;7r-Iv;5^HHhw=A1Qhzt{*PK2D zzBl{HftGP#etDX8^?(gHhkC1`&IJ=W-8bTQUWK);thZAlk#lrf_pysQ;CYnJx*|&CspE>uFOfH%8v{?O8UAY&W{9^?scLBa=Rl;_=Xwa`N@0xPS7o38dEn6Oq z2G9D{WW|T=u&;P?rt=5sV^8&r=Rd6n-_U2>qERMfe$qA9$0H4#Gt!0^+L9npE)(%g zSOTkx9a~eUN|V0V!25Z-4fw~_=IRvPfaMRz_0(sQf6uK=7ww(`p37GCU*g`dX{h;7 zctsza(x^K1FZaQlefapOXD7kDAvjQ@HV95xPM_n?RB%0%!`@JjV5^+4uO_mT^a4Q2#9~+kV_w#w7l7xzjOQbT+8lVgaQuOkw5m*16r7f&d?>{)Z9h|Etp zH?9S6l#9S^ocZ`_cL4Fzxvj~2M}c#zBkSy-Tmz zD6yAVW_|&uIo#FGbT_!Qaq`cSx?nfMSk{2Mgw)fAlEN4O_o24_K|6QQqvq-wrg6Y| zTvL{1KL*^VCq75*{|>8o)R$@4Am-d}aD2X~iIXZqObaM;|VUizsKoM#iJe_YE3 zXZX_mR~&aZY@54(?mJR%ywEV0awm4}5t~^$Wi6~npL2Z6(E;a`W2lk0B=P%{rxn{1=s&_;jQDvpewomQmb46PUpCd z12Kem!+ql%M^{)*qJb#%eMX6v$lcfb9N^PEh{OSL}bY2K34u zNA^#60nWDy>#8F~;Pm>wY6s+Qpv>^Kc<5x@X2x2F@fDI zkq)n+xuB+Ac=@V#Ex0||f+5{|WbRS>e0@9V&!<0?824Qk)K<>{{>jI%lIc6mw@3gb z^V?Ws;1%f;eK*VtBK<<|wt`1rSCIX}`2nj%NE~Gy*c#hv3~R&fJt|j>LCIQX?HT6{ zZeY5U&gD)ppDzo4!ZrjYzqKK6@Gl2<|Kqdj(2Wq0C41C$SU1X6@Nu!B9?ZyFM5VO5N-|PxCUU z`~x!ACW*pgAob|GCS6dAeW#b|M}fM2MtnEv`>kJ;tSwR&0j0?`d$x@1dA4;q%O2Pc z`z=q(SND*;F71-gF(#tKFTyihhIL>yoJ!5?egjJPVXm6R4RGVdVl9tp!ot^1Q*P23 zQ2M{_JJ_(neI7aCa@Q!(n?p1BE&8C0yYYA+_%~rm@YDD72{qAK!vfxkBnjN@V9?Srp9RXwx zPMQ6X^$aDvhn)2Jv{4jR^QSHSBzO(Va*j1Cg3P6!`n_9Xc^DSEAMYwTLFQZpg5|2H-)gV!f^xWPXtO^ATt2^_u3lkH~C4hKk~l*q3Lr_F3Ydg zuST|k=^EEFajEL=H{%xmK3OPnXX`~Eys zy`%wX39+1{^B<|eRzdM7;RQYom=sj3xU^S}*Iw0xN{3f9MiOqT5; zd(}G(UEcii25+ol`JkUHEIQ6RStvJy+O;le!f!S395=2UJH8zZdF7WyffGSR>@Cv~ z?E<&r+6J2?-LP-HDZVr+jqDlczIUBKc>fl-$zmLt&v*}>xvv~f^mt?XKFLtBAEKDQ z%gYdSvy|syX{Mm|b(x)*VE~?dUX)d$2rT8&_pZp=4k}i1d()hLaLsEZYh){76WwC{ zTWA6*UN7g{yb4k;$%U2%l064W{t?PjteX3B%mV<5wMWEs zF%8_u)2B@=-vJAu^0XZT#D6*2T7J4>;7J`h!!of3W9_L>n^Y1fxpf@XC!T@Mk|-aH*~ucI}7y%Xk7%$&)Kzn%RM`9bFl}N=WALP8yC1 zjpSZIyIahP?8B6q9{KJ@?#s^Q_c_#lC-(=dWf$8&08eD9z02yyu+4ILAJC9N=Ccvj z&a22i%tf=UUS+4iT={fytI-os2Md&^_g5bJ WbNT9@KmTap9}WDYf&Ys&@c#gSdeAul literal 20160 zcmeIYc{J7S+c%sN(LfO;L?uKSG8Xl*joY|Mp%fA+L&l;(neBOR%8)WPN~x3#DMN%p zgftIIN=Y=(_@3A0de;5B*Sdc1y5IM{*R!5yt?T=TwflCiS$G!(pd1o`d=4GdEV-=SbZbzY+o7VZ_a#Vd5n->fIEC{$+})%D-VGa&N3 zLPNcx6t;MUdntrR1^fKHl>c~ZZ{kW_L;m@<_%w;UpoqZDKA{RbeE*l^ZGk>PVcT{D zh5d(+fB9wZ>TGUr_HX1l|I5Pv;?F;Y{TJ8%v+N%Y{G)+?H1LlG{?Wic8u&*8|7hU< zehvJQPOhOk*c&S-Q>gaNCKOw9<$u|ghQj7;L0+Lz3gKRx1AG*GBE!d%(|=j=pMU;u z&7ZaZzvs_?=Hps_pFeBs5P1PUL0iN9{>aC-{I5)0fvo>`&!6(N|5l!^j^V$NN7nzl zhi2-q)}G|?|Al&CYv_)MV1;0>P_ICraG%h>uNN%;-{fiiC;j~?@9(Aj2U=$fs=2M% zn!m@l{6EPXHyPh;J^@?)y;NssH=6mM=b^mIYp2h+_W#ua|KrpDrC-h#4o-IC@<^5R z=kwo?rl8=w)}HD-E-xtXe?s0L`~A~ilzjstBL9&1-^2IEe#hng&+T_yp7)N3pm2qt zz@=LM4TOK5@A5zFSvexe+b1;KYg{=|JbjY ziL=Rn`+WbY&&My_KY#zxz&{%JM+5(-X#g?2_P*xN)4+cI#%H6VAt)DT+e>>UASOLx ziGq?P*v8U%Ul|E7HGdVk<(v~@vSmJ`K9vRgYi@(}15=nzJF(iGw-+(T>I~-GJOpmd z+L#gUZcx5YO#i&J7cpmo1oNH-gPVLtv$DMernj|bcNY>l7w4^85VsoK^~L7*w)?>J zh`|D3P99>)ZY|~bCxU11LcMot1g0~pr)j;WBIf!&=KZ1<;8uNlWTMg!Gev3f%Nu$S z)1VzPTDBa#0b_0R5iyt@v*qGo5n^t?dG>65F?e~JcV`DagsI1&oPxzFh`C#EKjE4f zQuHiC+M+(deBP0f_1C-*^T5;niC-D`&0n{?Huiv7C+DWa_j`zWKHvYDaSKvJzDlw_ zb;ER^Bv(LZ6=GgL7mnGt3H)?7rk`2`%r4Y8&wsQLF@wi7%WklcBIVBv)Rux-BX_$R zQyDRz{rmS#F9bio|CQ;2oiJ^=lK*VHUyHa#bxs>?8qU zhSx>#gGv^zvQ~uIqE|xm4BZhc(iU8&6Ab&-)+-g;yj1T&NR zs*HAf#40Vgmn}^L@8luXI>i~Pe?WRikX;6wEsYFv(Aa=XTkAr&M;Ow-Xx@$xml&@)fj~nD8Hbma@$&A}z z6_xJn;njdL*)DF=^?Qhol=5}(m3xX}U>}vS8Ov}b zI8VhL+Zu`3#2K#n=6k_9v1`oq0}Yg|dFmI6zay4C!!o9#4Q%?-i`^#eprr4(KdOHX zu?NK;#`|`Hv+GW*p2l-f)_wGMj`l@tmelORC7NK#CGib?=Ypc*ml2C) zV8s>1Q0@F+@_c#z)?zuto>n>F>HY=m#jn=nr@aDY%fX>3Mw=0PUiac~$T_fu+AG}J z%t6t$rYI;pM(h=HZec$UEZ0=ZdJz5C|0?ClTtMtKS7~OFJXmLs?aPoe10|p?!0_=& z#MTGeq=ggzz->}#7El6Z$Ct78$(s;+J3+nGAReq4eTCL4R-jmmUKo}2N9^64tGZRc zz_Ow$eVP;llTVKfBx@cZ_QBO}cW5PGwIo!njFkgr!RDm)EiV!K{L!1d>|Q=wK=p;-!zB{QZ9!Bi1U0GggM9$e8YvwC$gGpEU(H4Q3h?{&| zG`Y@|y#EoKC{zTK=9-66nUfGVtxI$A%ClgNw8|v(UxP_sRPB^0zY!-nW!j@pPr!~H z&2r2l`kd<~G<}T};^rCWOI`Q|&X(Hc8qS`exW_7*4@M$xLBvpbrU%&bn4W|10YxC$ zPCkZwSM~f@i_QyhY|nSxKe7xIObC8-$QN;%gV)?b`oP`4)lb}x14^9N;m={45oe&W z+Bc65?nVi}fp8&EGKv;9l!_s4b?CmIk#~syR zUtRfDtqLh>+;jK)T41Ieas84f1#zKi7H-rA@INWk?O)&s)3+xJ^+tve7xf|Sc;`2y z7#tQ!2_k+~@Rah-n~M<_V{`V=)Is8(zdihRl6;rdkXF2E3F49q*;7&uAw|htIq(&U z7oGQqhZg7|jyFx$AT$X4cS_P*hmB!s)&3^GOapQ0!6&0+GLYgd``Dv~=%Jy&rF# z-%Z4Ij4dc!_!(Tel+tJ3KS4SD#G-BB8{!`9e5#8n5Kv{2_V%2GexaVJk zP1luzBfOqo+ClEMrxQ_yoJ;hz%=sf!R z?PW745vfA8Rd$FQWh~jz*b0vQExL%-7%2PeXV|@;i}(prbAHyV6F*jyTK@epOj1u> z(sxWkyzsH-Y4cP_{;*_WDD^8Svbj6jTy+sIrk7xQdJb5(_njLqI|Y+-!K)6vK7)9v zTa$G!vxxsusg$SaL^Vo zOST5a6ds4ke0gL4*|LaVnprkDAktToG>;>lkJc3f7VR#|ItKVDgY1n(N_- zc-tkPrMsNLYR;P#I~oAW#I!xHl)oW9a>5@ESU0<+j|68ynNlTiFd!W9X97zf23`S2 z@VK@>RVZ|0z5SjTdn3M}F6@K;I>jx0Ez6N$Gmewj`B{8a-3IzIE^6pJ zzlVg)<2ZRUE4KaEPz}ApN{;TS3M5309YkT5Wgll3#^RzRN;bk7uoe0DvA#EIiLS=*kMWens)caM9M@cdm! zoI8$_Jz-Nn{iHT@6E}o6`#wgZ=3h7|tKJRuFM;k6)q(WHZX{ZaoaxlU3aQBaD6&x_d=#G+i}AQt&uV7G1%PPiMQl%EO<~ ztI)ev8gT~6;(zH=`+WIO;!o(u(}P8vtB|~S94Eg}Y2wmL6)+H+-6y-~B9e{Aai&Pe ze-_+n34;?QZ=g8|$@v*iVwLa%y>Eh|O{$(iFgIpW-7>vpe%-f6)GQLz#vpa07~#msf&-upvuS$qVg zA_2(_f9dnVf$aMoFQF%P()N6q6OteQf%73)$(p@NlfuF@i>;2x4c{dMjx8n+({ zzlNS!f>rMUdC*nE&}=* zz4FsJE?~V{Q0?w64$Gr`A7y8EfxdQwZmD=F((Kz{wu_B;qw04&^^OKM1*j zzFALJDv{)YiFHRMPA9=iUuk3Y{v^=1tWSURQ5l?i*;IAe&9JB$uBuSG1G?{iiGmGY z;Dx<#O}d&3yOArcMrAdi`xS7v)n|btwUuE#I0T2QUwV>P=z+eyajK>UiH}PbH&S}j zKn*=J@<`kY^niDXUk;0b>nf9eqD2E1_jDs3Ci{ZEW9GyvK`MAA8`u5zSPR=$zvlhc zexQfwM#lxkf%7rl+$13h4kPn7Y+1zxJH zaDPAFQs)~Fi{|f@VM}5_-Zy2?D$O|`L|({l|Ovo*$U93ZU!Ee4gj~jX7kNf zEzqdHS}t_gg1%?qjdu*m_m+OkcsIiUwyT~VODebqdbGHY>cBH_o2@oD^Ebg_^ug;% ztO3vu=sX)*Cj%~}`SQxVJ+NM^`Z;K-1L(1C?xoev;BHk|>Spf^nz`7vbHSTIkKfm` zZwUw7`nmmWBIjV2k~>+&^9ksQ`Oc;X8^KG}k^Z!!8kVQFNYtJ#1wFapj>zwc;59Fg zJZI7d>WxFo=7$^soiSi@XPXhpKQ%n^Y?pwh=eFO;A^>!@*v+)NJHfj>U48Q`AvnXK3i})qF7<%UcfG#wx^h0^2=l7Djer>;(UJMts?&A-=IV(X=Z_wV^mjRxs zOK_8@DXh8-PS}<3LC<`9&Y0E#p81y~)*bDUXr_qVZ#8wz+B8|WC1Ln_ZhtM4+j16jjy!J zH^ASQdS|MJ4{W+#CTJ@bfPUponDqQ#;9qw=Zd)e|N70x8Nuhn9S4eFOBX+S?00g@J#z^ysSzGhyYYmGk}BZqTo->)XEREck|)bZ^phSohQ& zv}~&fy?W1P&BoW@4=AjekaiD_nqOxe>m~l7Hs|rz<%&q1w25&_&K@);v%JZ<(V*YB z;Z}8Y1Ngt2FTIfofYsGL4G%^Y=nef{@xq6||Juo^o6rR6laQM_=ESa=raM}Pv>|ox z`}eV;i{W@sWPYS}BB4?pc^3WzX^Q>bA9@ z_ALkullB6=ZTF3%=l6j>df;4<_b~A@9XXPd+Cgv6!4{b-;Qa{o`lN6eG@JG-69O7R zzgu^8N&H8o96IYjQGEwTJGtcigFc{l_3ICPP$TCx`R%^pH?W?a`r^g2OwfC#m0WT2 z0ss3q2-Y&zvqhVT9tPGHE=q%>_I(IJovbgU!8NE9v!v}8I-erwgW-lW7i35?02 z4}OWOkUDY6M4LccSbl3tkG35FLs( zri?jnj8}poR%8BTp*6T=4+_VGufV=%;Ma@}Nib%yuULNQk=t_O$F3GC?yaSOp_aw@m}~^DRI|fA zH#N{AA2xsZ@&$~g)uLC1cM{%ybW^h@9ki)E%QqR`2Sej!Vzm2tu#SCk=6bcj%J``5 zoV_Q&&=#3!*6^0_P74m&tJJ}=f9Zza!zN(pERCM3Me5Mr_S753T5!0ooLByq=vmKk ztUYTJ*h!sPTC)d1^AqpVA590tU}r?eQ7hSBxFMGr`S0oL6%OE zMTNYrz%ZGz<-)lGVBe>@O{&L8vkDx4|+y6zagB9$s_HWg9r_ z2BJ+>f??-sMX}ZJ2E*oB8^5LitiBIpPpkIAGDWk47uODkU9XL++GDUZj3(5d-Un+V zbayUo1cN4gQ$9@^Y^NYY-h!Ouq$&lRD+q1V=C*NNMRbU5CRaYlAhKufTAo zT|2X#3ARpLr%>+)Sbb@)l=?~Re_h0ikj{Fr@9i&={Kbbgn$0H}_k!V?S!%HAH`uFx z`CXp&42~hC%Nsl*z*v7x>(^;Ru)Jn7Hfr+VxMXur!ii`w++Uonmk0xQjSxF)AOqH$ z>dqPp#ev}=tjgSX7Mzf(wSHOaKwVI0$m;$DhL_rL$JgW@bV)A>6w`zKjU=7MH~L_B z(-h}8DT1|kk;1;&#IEda-L5_71BOrdktbI-fVszOMP<1*sM{x=NOJ!G#@5U^#}^+0 zvm$qkN53>I9>xiO`28G=ZB+;Td2L`BZ&NC}mH}IroSCP53c>JyA+GgZ1dNWQ5-YZE z22K9C)|0#BT!Ms{qdul!$5gx7_N<4^NMsJ{CeeSen((!@ec(vUHD<^R!?IH_$w6f? z7@-bv35Hj|&S*0%R-6u7k^W<4Dpg>FhyAuG_y}fA!OBBR6+vqao?bq73XGi@duGj0 z1ItW6jel+t>~%~iz8^ioh$M-9qz_nKGUbkIykOaQ`ymg_(oMs(O%H-mC^)-l(h;!yIzR2Xn*{s8kDr9ny1*z>;fDs2{!!=UlkX|>L0kBy zGF4s~jB|F=4SK|hU71~}s%FCS`~H^aTHRor4@v$_u?9E?D==ur9!x;((BR33S{HULb6*1U4{L0Dh@ zw4v^jF&I_0JC0vJ0p`Q+hK}XL|1XlKX*v*pcs+O^K=m0|-;FhFe$RzX(NnWWlEk0a z9NMaVPzKCnw{p2-NihFDwpa59vBSFZ7o)$N!FcricKWh%KV1I8U34??JLTp+uik={HjH?=kFb20f530;JTN*R*QnE! zz@Ey{+<(6o7WBIPzf04=xHr*gxPtU+r0xVUc(Y&|JkXlcp$NwP#bxD-GQkS2F-}ul z4hQFRYb6E9`iC|;`#E?a|hSgSTA1hzf}d?{>_(_>9EC#v0p2=A?dq+8(?* zlHA9^R8ExjW76m9el_zV8O~Iis{~Zaj&D9iq*QCsCluc-N?d1E%;x z>!*YJh@R&^>h#+L2aA^PJMDIZIqP@B3D+y&zFOHf`|~MSv$oxm;FZl9wUWw*?JGIQpdBS(D*jfLfIYbJ%4usa7lTH7O8+h?3Fe}O$L<#< zgV(NNx;rZhw!M=>KQ^ucbFq~|u4FBE=a==WTFS$cB|bX2Lmy0)fFr)qKfvo>9<+T> z8P;n>Zk>NA1m+U{9Cc#%DG@t%3N5?>NA+hDvuUJnq<%SdWKA$q#EN>q&s{|9VfW06 zHlnv>55&r=$UH{<;X>DA7N8189{CnP1yl1E^MLY8GS~5?itRHI)RaD1`@Y#=E?+2Y zLHdGxhWzK<8b?9P`I_9>xkPF^UOM(kfEYUQ48FewYR zl|)?vZ!}oTO^xJBpUyuLSkDF1%;Hs4FX^9}EmqtfJ^?nhCZ8fg$i6K6yGITyV4LtF{K1R@}!zWt+iY_+(YF!!6inv#TCYkp|QHl55L?BjA6S(IjF14OR`) z6_Wa!z_jh@Ol5T=W%pHX{_Z2N4m>}+rZWUg`=2zoui{AAX09xt8w|&wh%+aB(!g|7 zyshA64SveEl`D-&y*&B0gX4$?m}@L7`)L6--y_pTht5Y+U<2`b+3UR z=i8$k0a)EQWgqT1jWCI;u~&w}H8 z^XI1)I)k~PM}I`8>~6lySup)dfSGDoxm&xj%5X_9vKO+me$6 z2f*K%qRG2;7moBH3$uAgz>MGs+U-yRuVZUtdt@z~MB9&1`_;kRb@712eGhQE$7UZ> zUI$w^!ArlKIkPtk_Ck(_9XA+(8SD3U zi#=qrgC{1(=L)52hd5MtG`QB@T`6K{b|8 zw%twk%bvgX?(YI{kM-|=R!Vr$-uqjN66(R^nzkR__L%U!a#HkHrJ!xpKYMTv@soT% zJI_vO@Vw;PQ-*@zbhUVaR>>SNQ&~;Qs~o|je^INFsDXo@w(*dm4w#27n7u#6B=Zwb z(;HW~?%fO8W9t_@hej|@bS+f)M&jUVpY>C6R6w zFpIa&T2WvI?j(nm!P66ABPlum%^oK(FEBX*)1$$i-n972gG5+1WV$EVk^PpO7rD8s z3|uwM9PwIzaNYq!XU<<&eG$ys;s6~F z5(l+7dy;-lg^h1tu!o*Jm^bde{??QXuHDs}Vj0AqU23#U;~#?AFzQ{iycQgeV0hRi z70{BF?3uGc5X>gor%6Xi9Mq`Vwe#`_oVK-!FZyf%<}H&A_EYr0zP++ox`G4Sk$|TB z#l%0i_}r5W4JP@+=*?9Esj#~6`Nt#M1TfneP7fXxg8d_X1=oed5w0gkD4y6!`?=OU ziu#*Lv#(d2o*8$QJUW;egk>42wvf)AhAs**sS;Kd!at zsJaO&s|jrbw@4n+BWv<0AqC9v@U64W`#{~aR*biJG2us0Rq=t;gRI#@UAn~1k+IDy ziUsB)AH$gSLVw`ZRrOo}$HPrS9RUl$d_uo$c>`eSOgfZeDGQt4N~5)M6<|I)r!`|4 ziRZTRyBu1DLH%kcCi!hTm@hhtIx{YTq2d+2gv_fsyxpBiX2`&N`Bf!r;@}_iW4b4D zYv9=ZY}eq~ATVFgJLVNg=2fVHo*^geVdvPEwme=0%r{U_zPA}{ZIh3Kd(~jIXydKo z8iIddOV*$XsRvu<^>sX(4eKXYHJg5uIQyP1ed)|SunY_L*9%_<&DZQjw+gZ64`);N zNX-N*RW58~-+9=xk(Z zRbGJADh=V|CyT-SGEb=0`8ugr^08W!^bME2S8^;Ter9xa?7?eeWZq_EUPDV6Xm1U( zj@SasZ{EK)s4gdUt7qtg=`OIbooKU`=!5wqdADL}Ac-fJ{L-g3z$!&&lA)hFn7@iX zzy3-7e_~?(VXJ5o@4_Pk&g=neLVL)$t!05IA9Gx5`wo zrWd(7zcL1EN9=Zn?_$zdx>R`bB$?+EyVE((VHa3$Bvu~yR1J$?(?0Mv3xYNCv%~Z9 z3t+P)H!padNzV81g339mU`fusb!y=NSnn5AR;n9-Cb@max8PE+q*hr5aY#Kl+iiN? z;%>rTW=IxfP{EqBxlU))1FQj;237A6Sk6A$yz>4KSTaegzHcV=;IrKVPoycZu3n|N z&x7EQJyTJ0`!3iar56jPRDxDf^@6Y`D|HAzrSR#Z zeIDVBzFBnE2j7I{&!X4DPNYs=ICt6XsWD)#a*uFwAp9poCAsSF4MlR!e~T z(mY&ou_#z-XXKV$BK6?9 z53&t^&aS4y!n@aK`ZP9JdS1yhZZg5z6t80bQVN#87?ehPQ?LvYCv|BNUi507vxpM| z4o{8TCVb8XYenJytRoXhU6%DUTt)y6{Xz$)B@mv%sP+41p%vh~i9Ivm{}@&$^$JcI z6W(Ff$S#%OaB!6ymKtp%a|K6e6Af3*0Lw&rc<{kFaK-C7{Vs!juodiuqbxHPc0sBZFDTtNExw}OR7i5#2Ojjitq&vB*hq;41C4-fkL zKhqHa%WnAoAy-jwLKmFRGba3Gs9R7euM8}j^x6%zW?=T+|Hyye081&8-oxjWz;asI zzCfMS4FPMic�<6td42wIzPX+0(W^t&H%BC5C5nWMTKB!?5!{v6FT2jl~m&!Imyw zIn9Ug72G3VtcwYc?0VWP!k^6VdG&T56qf_l={9qlkpsDJE!8W#NF2O6TK*t(4=kd# z9e@$(E4U9EOecEoOEowy-D{J?tp*-h%hU~OwzaO&z05-3lJ}<#`ZFRZmVGY{k`R_bJNIxJvj-7l* z9&E!o6Q3S6h0XPCIuj>06F+!L#9;-QA0D1?`>yODEXq#4x|=%!R^;u3*~bVkDjKC~ zNiQVvSN~Qej|SH6j}sovXdyfXeo|9XK%4ZY-z>)#tbMch=0^#Gz5L#FZC7JBxX173 zC~AYX-*D8g*95GuPcuVYtBHQnX1NUW!HU@wu3@(a?CI16wP{me@$%{>);x0Vaj_ph zU%5>BIFz}qS5--TPYu{7NA{m^(*LUB2XGfX6)ZQe15IV8m{waCvF}^2V_B-;EwYqS zl}Le2tVP#@31mO?A#bZM-r%%`eJb<%42!PlJSRco&smaB#XOG^{AUe%9%X~-Iu_q^ z+YBs@;rg!Td*J%YN}QU%1~h~0<`u`tdfvwFBf2xezC1Rh)#CvB!1IT{H<3BZlvt=MAKGa`+JVtnrX)gzUgoBp9WyIT^_?Ik6lf72U!A|_iA7%X_2a}zXq(FjfNt!C&8PM+P*4c4Xh8|DP6v;6Rf Date: Mon, 1 Jun 2026 01:09:47 -0600 Subject: [PATCH 12/18] fix type annotation with astropy modeling; add synth_data.rst --- docs/synth_data.rst | 182 +++++++++++++++++++++++++++++++++ specreduce/utils/synth_data.py | 12 +-- 2 files changed, 188 insertions(+), 6 deletions(-) create mode 100644 docs/synth_data.rst diff --git a/docs/synth_data.rst b/docs/synth_data.rst new file mode 100644 index 00000000..2a0e5c2d --- /dev/null +++ b/docs/synth_data.rst @@ -0,0 +1,182 @@ +.. _synth_data: + +Synthetic Data +============== + +The `specreduce.utils.synth_data` module provides `~specreduce.utils.synth_data.SynthImage`, +a composable builder for generating synthetic 2D spectroscopic images. It is primarily used to +create test data and documentation examples, but is also a convenient way to explore how the +reduction tools respond to known inputs. + +The builder +----------- + +`~specreduce.utils.synth_data.SynthImage` is **immutable** and **chainable**. You start from an +empty canvas of a given size, add *signal layers* and *noise*, and finally render the result with +one of the ``to_*`` methods. Each ``add_*`` call returns a *new* ``SynthImage``, so the original is +never modified and a base configuration can be branched safely:: + + base = SynthImage(nx=1024, ny=512).add_background(5) + arc = base.add_arcs(["HeI"]) # base is unchanged + spec = base.add_source() # an independent branch + +Signal layers are additive and stackable: + +- :meth:`~specreduce.utils.synth_data.SynthImage.add_background` adds a constant level. +- :meth:`~specreduce.utils.synth_data.SynthImage.add_source` adds a continuum source whose spatial + profile follows a Chebyshev trace (call it more than once for multiple sources). +- :meth:`~specreduce.utils.synth_data.SynthImage.add_arcs` adds emission lines from one or more + ``pypeit`` calibration line lists, with an optional cross-dispersion tilt. +- :meth:`~specreduce.utils.synth_data.SynthImage.add_skylines` is a convenience wrapper around + ``add_arcs`` for night-sky airglow (OH) line lists. + +Noise is applied last, regardless of the order in which it is added, in physical order +(Poisson, then read noise): + +- :meth:`~specreduce.utils.synth_data.SynthImage.add_poisson_noise` applies photon (Poisson) noise. +- :meth:`~specreduce.utils.synth_data.SynthImage.add_read_noise` (alias + :meth:`~specreduce.utils.synth_data.SynthImage.add_rdnoise`) adds Gaussian read noise. + +Both noise stages draw from a single generator seeded by the ``seed`` argument, so a seeded +``SynthImage`` renders reproducibly; with ``seed=None`` the noise is non-deterministic. + +Finally, render the image with +:meth:`~specreduce.utils.synth_data.SynthImage.to_array` (a `~numpy.ndarray`), +:meth:`~specreduce.utils.synth_data.SynthImage.to_ccddata` (a `~astropy.nddata.CCDData`), or +:meth:`~specreduce.utils.synth_data.SynthImage.to_spectrum` (a `~specutils.Spectrum`). + +The following example builds a traced continuum source with a constant background, Poisson noise, +and read noise, then renders it to a `~astropy.nddata.CCDData`: + +.. plot:: + :include-source: + + import matplotlib.pyplot as plt + from astropy.modeling import models + from specreduce.utils.synth_data import SynthImage + + image = ( + SynthImage(nx=1024, ny=400, seed=42) + .add_background(5) + .add_source(profile=models.Moffat1D(amplitude=20, alpha=0.1)) + .add_poisson_noise() + .add_read_noise(3) + .to_ccddata() + ) + plt.figure(figsize=(10, 4)) + plt.imshow(image, origin="lower", aspect="auto") + +Tilted and curved arc lines +--------------------------- + +This is an example of modeling a spectrograph whose output is curved in the +cross-dispersion direction: + +.. plot:: + :include-source: + + import matplotlib.pyplot as plt + import numpy as np + from astropy.modeling import models + import astropy.units as u + from specreduce.utils.synth_data import SynthImage + + model_deg2 = models.Legendre1D(degree=2, c0=50, c1=0, c2=100) + im = ( + SynthImage(nx=3000, ny=1000) + .add_background(5) + .add_arcs(['HeI', 'ArI', 'ArII'], line_fwhm=3, tilt_func=model_deg2) + .add_poisson_noise() + .to_ccddata() + ) + fig = plt.figure(figsize=(10, 6)) + plt.imshow(im) + +Modeling a non-linear dispersion relation +----------------------------------------- + +The FITS WCS standard implements ideal world coordinate functions based on the physics +of simple dispersers. This is described in detail by Paper III, +https://www.aanda.org/articles/aa/pdf/2006/05/aa3818-05.pdf. This can be used to model a +non-linear dispersion relation based on the properties of a spectrograph. This example +recreates Figure 5 in that paper using a spectrograph with a 450 lines/mm volume phase +holographic grism. Standard gratings only use the first three ``PV`` terms: + +.. plot:: + :include-source: + + import numpy as np + import matplotlib.pyplot as plt + from astropy.wcs import WCS + import astropy.units as u + from specreduce.utils.synth_data import SynthImage + + non_linear_header = { + 'CTYPE1': 'AWAV-GRA', # Grating dispersion function with air wavelengths + 'CUNIT1': 'Angstrom', # Dispersion units + 'CRPIX1': 719.8, # Reference pixel [pix] + 'CRVAL1': 7245.2, # Reference value [Angstrom] + 'CDELT1': 2.956, # Linear dispersion [Angstrom/pix] + 'PV1_0': 4.5e5, # Grating density [1/m] + 'PV1_1': 1, # Diffraction order + 'PV1_2': 27.0, # Incident angle [deg] + 'PV1_3': 1.765, # Reference refraction + 'PV1_4': -1.077e6, # Refraction derivative [1/m] + 'CTYPE2': 'PIXEL', # Spatial detector coordinates + 'CUNIT2': 'pix', # Spatial units + 'CRPIX2': 1, # Reference pixel + 'CRVAL2': 0, # Reference value + 'CDELT2': 1 # Spatial units per pixel + } + + linear_header = { + 'CTYPE1': 'AWAV', # Grating dispersion function with air wavelengths + 'CUNIT1': 'Angstrom', # Dispersion units + 'CRPIX1': 719.8, # Reference pixel [pix] + 'CRVAL1': 7245.2, # Reference value [Angstrom] + 'CDELT1': 2.956, # Linear dispersion [Angstrom/pix] + 'CTYPE2': 'PIXEL', # Spatial detector coordinates + 'CUNIT2': 'pix', # Spatial units + 'CRPIX2': 1, # Reference pixel + 'CRVAL2': 0, # Reference value + 'CDELT2': 1 # Spatial units per pixel + } + + non_linear_wcs = WCS(non_linear_header) + linear_wcs = WCS(linear_header) + + # this re-creates Paper III, Figure 5 + pix_array = 200 + np.arange(1400) + nlin = non_linear_wcs.spectral.pixel_to_world(pix_array) + lin = linear_wcs.spectral.pixel_to_world(pix_array) + resid = (nlin - lin).to(u.Angstrom) + plt.plot(pix_array, resid) + plt.xlabel("Pixel") + plt.ylabel("Correction (Angstrom)") + plt.show() + + nlin_im = ( + SynthImage(nx=600, ny=512, wcs=non_linear_wcs) + .add_background(5) + .add_arcs(['HeI', 'NeI'], line_fwhm=3, wave_air=True) + .add_poisson_noise() + .to_ccddata() + ) + lin_im = ( + SynthImage(nx=600, ny=512, wcs=linear_wcs) + .add_background(5) + .add_arcs(['HeI', 'NeI'], line_fwhm=3, wave_air=True) + .add_poisson_noise() + .to_ccddata() + ) + + # subtracting the linear simulation from the non-linear one shows how the + # positions of lines diverge between the two cases + plt.imshow(nlin_im.data - lin_im.data) + plt.show() + +.. note:: + + The older ``make_2d_trace_image``, ``make_2d_arc_image``, and ``make_2d_spec_image`` + functions are deprecated in favor of `~specreduce.utils.synth_data.SynthImage` and will be + removed in a future release. diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 315cb254..fea7daed 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -136,7 +136,7 @@ def add_background(self, level: float) -> "SynthImage": def add_source( self, - profile: Model | None = None, + profile: Model = None, trace_center: float | None = None, trace_order: int = 3, trace_coeffs: dict | None = None, @@ -153,7 +153,7 @@ def add_arcs( line_fwhm: float = 5.0, amplitude_scale: float = 1.0, wave_air: bool = False, - tilt_func: Model | None = None, + tilt_func: Model = None, ) -> "SynthImage": """Add emission lines from one or more pypeit calibration line lists.""" if tilt_func is None: @@ -329,7 +329,7 @@ def make_2d_trace_image( trace_center: int | float | None = None, trace_order: int = 3, trace_coeffs: dict | None = None, - profile: Model | None = None, + profile: Model = None, add_noise: bool = True, ) -> CCDData: """Deprecated. Use :class:`SynthImage` instead. @@ -367,7 +367,7 @@ def make_2d_arc_image( line_fwhm: float = 5.0, linelists=("HeI",), amplitude_scale: float = 1.0, - tilt_func: Model | None = None, + tilt_func: Model = None, add_noise: bool = True, ) -> CCDData: """Deprecated. Use :class:`SynthImage` with ``.add_arcs(...)`` instead.""" @@ -401,11 +401,11 @@ def make_2d_spec_image( line_fwhm: float = 5.0, linelists=("OH_GMOS",), amplitude_scale: float = 1.0, - tilt_func: Model | None = None, + tilt_func: Model = None, trace_center: int | float | None = None, trace_order: int = 3, trace_coeffs: dict | None = None, - source_profile: Model | None = None, + source_profile: Model = None, add_noise: bool = True, ) -> CCDData: """Deprecated. Use :class:`SynthImage` with ``.add_arcs(...)`` and From 01d5366be1e0333ceb9365462a11142b7a2fcdb2 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 16:18:26 -0600 Subject: [PATCH 13/18] minor text and formatting edits --- specreduce/utils/synth_data.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index fea7daed..f7494810 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -92,7 +92,7 @@ class SynthImage: .to_ccddata() ) - See the :ref:`synth_data` guide for worked examples, including tilted arc + See the :ref:`synth_data` guide for more examples, including tilted arc lines and modeling a non-linear dispersion relation. """ @@ -252,7 +252,8 @@ def to_spectrum(self) -> Spectrum: @dataclass(frozen=True) class SourceLayer: - """A continuum source whose spatial profile follows a Chebyshev trace. + """ + A continuum source whose spatial profile follows a Chebyshev trace. The dispersion axis is the X (column) axis, matching the historical ``make_2d_trace_image`` behaviour (the trace ignores any WCS). @@ -276,7 +277,8 @@ def render(self, ctx: _RenderContext) -> np.ndarray: @dataclass(frozen=True) class ArcLayer: - """Emission lines from one or more pypeit calibration line lists. + """ + Emission lines from one or more pypeit calibration line lists. Requires a resolvable WCS (supplied on ``SynthImage`` or built from its ``extent``). ``tilt_func`` applies a cross-dispersion tilt to simulate @@ -332,7 +334,8 @@ def make_2d_trace_image( profile: Model = None, add_noise: bool = True, ) -> CCDData: - """Deprecated. Use :class:`SynthImage` instead. + """ + Deprecated. Use :class:`SynthImage` instead. Equivalent to ``SynthImage(nx, ny).add_background(background) .add_source(...)`` followed by ``.add_poisson_noise()`` when ``add_noise``, From fb6cd3ac0f9a58a9b64bb61ecbeb4a0363c45aef Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 17:12:45 -0600 Subject: [PATCH 14/18] add spectrum argument to SynthImage.add_source for wavelength-dependent flux --- docs/synth_data.rst | 7 ++- specreduce/tests/test_synth_data.py | 63 ++++++++++++++++++++++++ specreduce/utils/synth_data.py | 76 ++++++++++++++++++++++++++--- 3 files changed, 137 insertions(+), 9 deletions(-) diff --git a/docs/synth_data.rst b/docs/synth_data.rst index 2a0e5c2d..450a5317 100644 --- a/docs/synth_data.rst +++ b/docs/synth_data.rst @@ -23,8 +23,11 @@ never modified and a base configuration can be branched safely:: Signal layers are additive and stackable: - :meth:`~specreduce.utils.synth_data.SynthImage.add_background` adds a constant level. -- :meth:`~specreduce.utils.synth_data.SynthImage.add_source` adds a continuum source whose spatial - profile follows a Chebyshev trace (call it more than once for multiple sources). +- :meth:`~specreduce.utils.synth_data.SynthImage.add_source` adds a source whose spatial + profile follows a Chebyshev trace (call it more than once for multiple sources). By default the + source is a flat continuum, but passing a 1D `~specutils.Spectrum` modulates the flux along the + dispersion axis: the spectrum is resampled onto the image wavelength grid, normalized to a peak of + one, and zeroed outside its wavelength range. - :meth:`~specreduce.utils.synth_data.SynthImage.add_arcs` adds emission lines from one or more ``pypeit`` calibration line lists, with an optional cross-dispersion tilt. - :meth:`~specreduce.utils.synth_data.SynthImage.add_skylines` is a convenience wrapper around diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index adcb3be1..14107a07 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -73,6 +73,69 @@ def test_add_source_stackable(): assert two.to_array().sum() > one.to_array().sum() +def test_add_source_spectrum_modulates_and_normalizes(): + nx, ny = 100, 50 + extent = (3000, 6000) + wave = np.linspace(extent[0], extent[1], 400) * u.AA + flux = np.exp(-0.5 * ((wave.value - 4500) / 300) ** 2) * u.count + spectrum = Spectrum(flux=flux, spectral_axis=wave) + img = ( + SynthImage(nx=nx, ny=ny, extent=extent) + .add_source( + profile=models.Gaussian1D(amplitude=10, stddev=5), + trace_coeffs={"c0": 0, "c1": 0, "c2": 0}, + spectrum=spectrum, + ) + .to_array() + ) + # peak value equals the profile amplitude (normalized flux peak == 1) + assert np.isclose(img.max(), 10.0) + # flux varies along the dispersion (X) axis -- not a flat continuum + col_sums = img.sum(axis=0) + assert not np.allclose(col_sums, col_sums[0]) + # most flux concentrated near the 4500 A peak + waves = _linear_wcs(nx, extent=extent).spectral.pixel_to_world(np.arange(nx)).to(u.AA).value + assert col_sums.argmax() == np.abs(waves - 4500).argmin() + + +def test_add_source_spectrum_cropped_to_extent(): + nx, ny = 120, 40 + extent = (3000, 6000) + wcs = _linear_wcs(nx, extent=extent) + waves = wcs.spectral.pixel_to_world(np.arange(nx)).to(u.AA).value + wave = np.linspace(4000, 5000, 200) * u.AA + flux = np.ones(200) * u.count + spectrum = Spectrum(flux=flux, spectral_axis=wave) + img = SynthImage(nx=nx, ny=ny, wcs=wcs).add_source(spectrum=spectrum).to_array() + # flux outside the spectrum's range is zero + assert np.all(img[:, waves < 4000] == 0) + assert np.all(img[:, waves > 5000] == 0) + # flux present within the spectrum's range + assert img[:, (waves > 4400) & (waves < 4600)].sum() > 0 + + +def test_add_source_spectrum_builds_wcs_from_extent(): + extent = (3000, 6000) + wave = np.linspace(extent[0], extent[1], 100) * u.AA + flux = np.ones(100) * u.count + spectrum = Spectrum(flux=flux, spectral_axis=wave) + img = SynthImage(nx=80, ny=30, extent=extent).add_source(spectrum=spectrum) + arr = img.to_array() + assert arr.shape == (30, 80) + assert arr.max() > 0 + # the built WCS is carried through to the rendered CCDData + assert img.to_ccddata().wcs is not None + + +def test_add_source_spectrum_requires_wcs_or_extent(): + wave = np.linspace(3000, 6000, 50) * u.AA + flux = np.ones(50) * u.count + spectrum = Spectrum(flux=flux, spectral_axis=wave) + img = SynthImage(nx=60, ny=20, extent=None).add_source(spectrum=spectrum) + with pytest.raises(ValueError, match="Must specify either a wavelength extent or a WCS"): + img.to_array() + + def _linear_wcs(nx, extent=(3000, 6000), wave_unit=u.Angstrom): wcs = WCS(naxis=2) wcs.wcs.ctype[0] = "WAVE" diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index f7494810..70810ecd 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -140,11 +140,35 @@ def add_source( trace_center: float | None = None, trace_order: int = 3, trace_coeffs: dict | None = None, + spectrum: Spectrum = None, ) -> "SynthImage": - """Add a continuum source with a Chebyshev-traced spatial profile.""" + """ + Add a source with a Chebyshev-traced spatial profile. + + Parameters + ---------- + profile + Astropy model describing the cross-dispersion spatial profile of the + source. Defaults to ``Moffat1D(amplitude=10, alpha=0.1)``. + trace_center + Central cross-dispersion position of the trace. Defaults to the + middle of the image. + trace_order + Polynomial order of the Chebyshev trace. + trace_coeffs + Coefficients of the Chebyshev trace, e.g. ``{"c0": 0, "c1": 50}``. + spectrum + Optional 1D `~specutils.Spectrum` describing the dispersion-axis flux + of the source. Its flux is resampled onto the image wavelength grid + (requiring a resolvable WCS) and normalized so its peak within the + image extent is one. Wavelengths outside the spectrum's range are set + to zero. The normalized flux multiplies the spatial profile column by + column, so the source amplitude varies with wavelength. When ``None`` + (default) the source has a flat continuum. + """ if profile is None: profile = models.Moffat1D(amplitude=10, alpha=0.1) - layer = SourceLayer(profile, trace_center, trace_order, trace_coeffs) + layer = SourceLayer(profile, trace_center, trace_order, trace_coeffs, spectrum) return self._clone(_layers=self._layers + (layer,)) def add_arcs( @@ -180,14 +204,18 @@ def add_read_noise(self, sigma: float) -> "SynthImage": add_rdnoise = add_read_noise def _resolve_wcs(self): - has_arc = any(isinstance(layer, ArcLayer) for layer in self._layers) + needs_wcs = any( + isinstance(layer, ArcLayer) + or (isinstance(layer, SourceLayer) and layer.spectrum is not None) + for layer in self._layers + ) if self._wcs is not None: wcs = self._wcs if wcs.spectral.naxis != 1: raise ValueError("Provided WCS must have a spectral axis.") if wcs.naxis != 2: raise ValueError("WCS must have NAXIS=2 for a 2D image.") - elif has_arc: + elif needs_wcs: if self._extent is None: raise ValueError("Must specify either a wavelength extent or a WCS.") if len(self._extent) != 2: @@ -253,15 +281,18 @@ def to_spectrum(self) -> Spectrum: @dataclass(frozen=True) class SourceLayer: """ - A continuum source whose spatial profile follows a Chebyshev trace. + A source whose spatial profile follows a Chebyshev trace. The dispersion axis is the X (column) axis, matching the historical - ``make_2d_trace_image`` behaviour (the trace ignores any WCS). + ``make_2d_trace_image`` behaviour (the trace ignores any WCS). When a + ``spectrum`` is supplied, its normalized flux modulates the profile along + the dispersion axis; otherwise the source has a flat continuum. """ profile: Model trace_center: float | None = None trace_order: int = 3 trace_coeffs: dict | None = None + spectrum: Spectrum | None = None def render(self, ctx: _RenderContext) -> np.ndarray: trace_center = ctx.ny / 2 if self.trace_center is None else self.trace_center @@ -272,7 +303,38 @@ def render(self, ctx: _RenderContext) -> np.ndarray: ) trace_mod = models.Chebyshev1D(degree=self.trace_order, **trace_coeffs) trace = ctx.yy - trace_center + trace_mod(ctx.xx / ctx.nx) - return self.profile(trace) + image = self.profile(trace) + if self.spectrum is not None: + image = image * self._normalized_flux(ctx)[np.newaxis, :] + return image + + def _normalized_flux(self, ctx: _RenderContext) -> np.ndarray: + """ + Resample the source spectrum onto the image's dispersion-axis wavelengths. + + Returns a 1D array of length ``ctx.nx``: the spectrum's flux interpolated + onto the image wavelength grid, zero outside the spectrum's range, and + normalized so its peak within the image extent is one. + """ + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", message="No observer defined on WCS.*") + image_waves = ctx.wcs.spectral.pixel_to_world(np.arange(ctx.nx)) + spec_waves = self.spectrum.spectral_axis.to( + image_waves.unit, equivalencies=u.spectral() + ) + flux = np.asarray(self.spectrum.flux.value, dtype=float) + order = np.argsort(spec_waves.value) + resampled = np.interp( + image_waves.value, + spec_waves.value[order], + flux[order], + left=0.0, + right=0.0, + ) + peak = resampled.max() + if peak > 0: + resampled = resampled / peak + return resampled @dataclass(frozen=True) From 7712786c43722f8b6b877f4db8b424aba53227de Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 17:14:35 -0600 Subject: [PATCH 15/18] add spectrum-modulated source example to synth_data guide --- docs/synth_data.rst | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/docs/synth_data.rst b/docs/synth_data.rst index 450a5317..3c74b73d 100644 --- a/docs/synth_data.rst +++ b/docs/synth_data.rst @@ -69,6 +69,51 @@ and read noise, then renders it to a `~astropy.nddata.CCDData`: plt.figure(figsize=(10, 4)) plt.imshow(image, origin="lower", aspect="auto") +A source with a wavelength-dependent spectrum +--------------------------------------------- + +By default a source is a flat continuum, but passing a 1D `~specutils.Spectrum` to +:meth:`~specreduce.utils.synth_data.SynthImage.add_source` modulates its flux along the dispersion +axis. The spectrum is resampled onto the image wavelength grid, normalized to a peak of one, and +zeroed outside its wavelength range. The example below builds a continuum with two emission lines and +renders a traced source whose brightness follows it (no network access required): + +.. plot:: + :include-source: + + import matplotlib.pyplot as plt + import numpy as np + import astropy.units as u + from astropy.modeling import models + from specutils import Spectrum + from specreduce.utils.synth_data import SynthImage + + # a rising continuum with two emission lines + wave = np.linspace(4000, 7000, 1000) * u.Angstrom + continuum = 1.0 + 0.3 * (wave.value - 4000) / 3000 + emission = ( + 2.0 * np.exp(-0.5 * ((wave.value - 5000) / 20) ** 2) + + 1.5 * np.exp(-0.5 * ((wave.value - 6300) / 30) ** 2) + ) + spectrum = Spectrum(flux=(continuum + emission) * u.count, spectral_axis=wave) + + image = ( + SynthImage(nx=1024, ny=300, extent=(4000, 7000), seed=42) + .add_background(5) + .add_source(profile=models.Moffat1D(amplitude=50, alpha=0.1), spectrum=spectrum) + .add_poisson_noise() + .to_ccddata() + ) + + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5)) + ax1.plot(wave, spectrum.flux) + ax1.set_xlabel("Wavelength (Angstrom)") + ax1.set_ylabel("Input flux") + ax2.imshow(image, origin="lower", aspect="auto") + ax2.set_xlabel("Dispersion axis (pix)") + ax2.set_ylabel("Cross-disp. axis (pix)") + fig.tight_layout() + Tilted and curved arc lines --------------------------- From 90eae444a46a71e45e8c692e3d2c108e373a7642 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 17:35:01 -0600 Subject: [PATCH 16/18] add test to help complete code coverage --- specreduce/tests/test_synth_data.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index 14107a07..8383388b 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -114,6 +114,22 @@ def test_add_source_spectrum_cropped_to_extent(): assert img[:, (waves > 4400) & (waves < 4600)].sum() > 0 +def test_add_source_spectrum_outside_extent_is_zero(): + nx, ny = 80, 30 + extent = (3000, 6000) + # spectrum lies entirely outside the image's wavelength range -> no overlap + wave = np.linspace(8000, 9000, 100) * u.AA + flux = np.ones(100) * u.count + spectrum = Spectrum(flux=flux, spectral_axis=wave) + img = ( + SynthImage(nx=nx, ny=ny, extent=extent) + .add_source(profile=models.Gaussian1D(amplitude=10, stddev=5), spectrum=spectrum) + .to_array() + ) + # nothing overlaps, so the normalized flux (and the source) is all zero + assert np.all(img == 0) + + def test_add_source_spectrum_builds_wcs_from_extent(): extent = (3000, 6000) wave = np.linspace(extent[0], extent[1], 100) * u.AA From a8d432f86dd1aea99ef40e885236edbf0a1d3cd1 Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Mon, 1 Jun 2026 17:52:48 -0600 Subject: [PATCH 17/18] cover disp_axis==1 arc rendering and deprecated-shim noise branches --- specreduce/tests/test_synth_data.py | 52 +++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index 8383388b..a517f9ba 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -165,6 +165,20 @@ def _linear_wcs(nx, extent=(3000, 6000), wave_unit=u.Angstrom): return wcs +def _transposed_wcs(ny, extent=(3000, 6000), wave_unit=u.Angstrom): + """A WCS whose spectral axis is the *second* axis (dispersion along Y).""" + wcs = WCS(naxis=2) + wcs.wcs.ctype[0] = "PIXEL" + wcs.wcs.ctype[1] = "WAVE" + wcs.wcs.cunit[0] = u.pixel + wcs.wcs.cunit[1] = wave_unit + wcs.wcs.crval[0] = 0 + wcs.wcs.cdelt[0] = 1 + wcs.wcs.crval[1] = extent[0] + wcs.wcs.cdelt[1] = (extent[1] - extent[0]) / ny + return wcs + + # --- local (no network) validation tests --- def test_arc_requires_wcs_or_extent(): @@ -233,6 +247,18 @@ def test_arc_image_supplied_wcs_and_tilt(): assert ccd.data.shape == (100, 300) +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_arcs_dispersion_along_y_axis(): + # spectral axis is the second WCS axis -> disp_axis == 1 branches in ArcLayer + ny = 300 + wcs = _transposed_wcs(ny) + tilt = models.Chebyshev1D(degree=2, c0=20, c1=0, c2=40) + img = SynthImage(nx=100, ny=ny, wcs=wcs).add_arcs(["HeI"], tilt_func=tilt).to_array() + assert img.shape == (ny, 100) + assert img.max() > 0 + + @pytest.mark.remote_data @pytest.mark.filterwarnings("ignore:No observer defined on WCS") def test_arcs_stackable(): @@ -349,3 +375,29 @@ def test_make_2d_spec_image_deprecated_and_matches(): assert isinstance(ccd, CCDData) assert ccd.data.shape == (100, 300) assert np.array_equal(ccd.data, expected) + + +def test_make_2d_trace_image_add_noise(): + # exercises the add_noise=True (Poisson) branch of the deprecated shim + with pytest.warns(AstropyDeprecationWarning): + ccd = make_2d_trace_image(nx=100, ny=50, add_noise=True) + assert isinstance(ccd, CCDData) + assert ccd.data.shape == (50, 100) + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_make_2d_arc_image_add_noise(): + with pytest.warns(AstropyDeprecationWarning): + ccd = make_2d_arc_image(nx=300, ny=100, add_noise=True) + assert isinstance(ccd, CCDData) + assert ccd.data.shape == (100, 300) + + +@pytest.mark.remote_data +@pytest.mark.filterwarnings("ignore:No observer defined on WCS") +def test_make_2d_spec_image_add_noise(): + with pytest.warns(AstropyDeprecationWarning): + ccd = make_2d_spec_image(nx=300, ny=100, add_noise=True) + assert isinstance(ccd, CCDData) + assert ccd.data.shape == (100, 300) From 079c63d6a3bd5ddb5b39c6e901cab56cc89d334c Mon Sep 17 00:00:00 2001 From: "T. E. Pickering" Date: Sat, 6 Jun 2026 15:28:05 -0600 Subject: [PATCH 18/18] fix: adapt default trace_coeffs to trace_order in SynthImage.add_source --- specreduce/tests/test_synth_data.py | 9 +++++++++ specreduce/utils/synth_data.py | 15 ++++++++++----- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/specreduce/tests/test_synth_data.py b/specreduce/tests/test_synth_data.py index a517f9ba..22e87e5a 100644 --- a/specreduce/tests/test_synth_data.py +++ b/specreduce/tests/test_synth_data.py @@ -65,6 +65,15 @@ def test_add_source_default_profile(): assert arr.max() > 0 +@pytest.mark.parametrize("trace_order", [0, 1, 2, 3]) +def test_add_source_default_coeffs_low_trace_order(trace_order): + # The default trace_coeffs must adapt to trace_order; orders < 2 previously + # raised because the hard-coded defaults included c1/c2. See PR #311. + arr = SynthImage(nx=100, ny=50).add_source(trace_order=trace_order).to_array() + assert arr.shape == (50, 100) + assert arr.max() > 0 + + def test_add_source_stackable(): one = SynthImage(nx=200, ny=100).add_source( profile=models.Gaussian1D(amplitude=100, stddev=10) diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 70810ecd..1d5f28de 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -296,11 +296,16 @@ class SourceLayer: def render(self, ctx: _RenderContext) -> np.ndarray: trace_center = ctx.ny / 2 if self.trace_center is None else self.trace_center - trace_coeffs = ( - {"c0": 0, "c1": 50, "c2": 100} - if self.trace_coeffs is None - else self.trace_coeffs - ) + if self.trace_coeffs is None: + # Default to a curved trace, but only keep coefficients that fit + # within the requested order so that trace_order < 2 still works. + trace_coeffs = { + f"c{i}": c + for i, c in enumerate((0, 50, 100)) + if i <= self.trace_order + } + else: + trace_coeffs = self.trace_coeffs trace_mod = models.Chebyshev1D(degree=self.trace_order, **trace_coeffs) trace = ctx.yy - trace_center + trace_mod(ctx.xx / ctx.nx) image = self.profile(trace)