Refactor synthetic data utilities

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

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
-}
+}

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(
@@ -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),
@@ -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
 

docs/synth_data.rst

@@ -0,0 +1,230 @@
+.. _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 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
+  ``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")
+
+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
+---------------------------
+
+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.

docs/user_guide.rst

@@ -18,4 +18,5 @@ User Guide
 
     extinction
     specphot_standards
-    mask_treatment/mask_treatment
+    mask_treatment/mask_treatment
+    synth_data