diff --git a/.github/workflows/build_workshop.yml b/.github/workflows/build_workshop.yml index b429f10b..6c0e9fd8 100644 --- a/.github/workflows/build_workshop.yml +++ b/.github/workflows/build_workshop.yml @@ -31,13 +31,13 @@ jobs: if: ${{ !startsWith(github.ref_name, 'v') }} run: echo 'BRANCHNAME=latest' >> "$GITHUB_ENV" - - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + - uses: actions/checkout@9c091bb21b7c1c1d1991bb908d89e4e9dddfe3e0 # v7.0.0 with: fetch-depth: 0 persist-credentials: false - name: Setup Python 3.13 - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0 + uses: actions/setup-python@a309ff8b426b58ec0e2a45f0f869d46889d02405 # v6.2.0 with: python-version: "3.13" @@ -50,7 +50,7 @@ jobs: run: pip install -U pipx - name: Clone workshop repository - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + uses: actions/checkout@9c091bb21b7c1c1d1991bb908d89e4e9dddfe3e0 # v7.0.0 with: persist-credentials: true repository: 'DKISTDC/DKIST-Workshop' diff --git a/changelog/713.feature.rst b/changelog/713.feature.rst new file mode 100644 index 00000000..4d462e24 --- /dev/null +++ b/changelog/713.feature.rst @@ -0,0 +1 @@ +Added `dkist.wcs.models.build_grating_spectral_transform()` to build FITS ``-GRA``/``-GRI`` spectral transforms by composing a constant incident-angle term, `~gwcs.spectroscopy.RefractedAngleSineModel`, and `~gwcs.spectroscopy.WavelengthFromGrismEquation` using the FITS-facing grating parameters. diff --git a/dkist/wcs/models.py b/dkist/wcs/models.py index f19fb75b..85cdda42 100755 --- a/dkist/wcs/models.py +++ b/dkist/wcs/models.py @@ -14,6 +14,7 @@ import astropy.units as u from astropy.modeling import CompoundModel, Model, Parameter, separable from astropy.utils.decorators import deprecated_renamed_argument +from gwcs.spectroscopy import RefractedAngleSineModel, WavelengthFromGrismEquation from dkist.utils.decorators import deprecated from dkist.utils.exceptions import DKISTDeprecationWarning @@ -30,11 +31,63 @@ "VaryingCelestialTransform", "VaryingCelestialTransform2D", "VaryingCelestialTransform3D", + "build_grating_spectral_transform", "generate_celestial_transform", "varying_celestial_transform_from_tables", ] +def build_grating_spectral_transform( + reference_pixel: float, + reference_wavelength: u.Quantity, + dispersion: u.Quantity, + groove_density: u.Quantity, + spectral_order: u.Quantity, + incident_angle: u.Quantity, + refractive_index: u.Quantity = 1 * u.one, + refractive_index_derivative: u.Quantity = 0 / u.m, + out_of_plane_angle: u.Quantity = 0 * u.deg, + camera_angle: u.Quantity = 0 * u.deg, +) -> CompoundModel: + """ + Build a FITS grating spectral transform from header-derived parameters. + + Composes a constant incident-angle sine term (`~astropy.modeling.models.Const1D`), + a pixel-dependent refracted-angle sine model + (`~gwcs.spectroscopy.RefractedAngleSineModel`), and + `~gwcs.spectroscopy.WavelengthFromGrismEquation`, following the FITS + grating/grism spectral-coordinate formalism described by Greisen et al. + (2006). The input and output angles are computed from the Greisen + relations within the component models: + https://scixplorer.org/abs/2006A%26A...446..747G/abstract + """ + model = WavelengthFromGrismEquation( + groove_density=groove_density, + spectral_order=spectral_order, + reference_wavelength=reference_wavelength, + refractive_index=refractive_index, + refractive_index_derivative=refractive_index_derivative, + out_of_plane_angle=out_of_plane_angle, + ) + + alpha_in = m.Const1D(amplitude=np.sin(incident_angle)) + + alpha_out = RefractedAngleSineModel( + reference_pixel=reference_pixel, + reference_wavelength=reference_wavelength, + dispersion=dispersion, + groove_density=groove_density, + spectral_order=spectral_order, + incident_angle=incident_angle, + refractive_index=refractive_index, + refractive_index_derivative=refractive_index_derivative, + out_of_plane_angle=out_of_plane_angle, + camera_angle=camera_angle, + ) + + return m.Mapping((0, 0)) | (alpha_in & alpha_out) | model + + def generate_celestial_transform( crpix: Iterable[float] | u.Quantity, cdelt: Iterable[float] | u.Quantity, diff --git a/dkist/wcs/tests/test_models.py b/dkist/wcs/tests/test_models.py index b96e164e..8effebc2 100755 --- a/dkist/wcs/tests/test_models.py +++ b/dkist/wcs/tests/test_models.py @@ -9,10 +9,13 @@ from astropy.coordinates.matrix_utilities import rotation_matrix from astropy.modeling import CompoundModel from astropy.modeling.models import Tabular1D +from astropy.wcs import WCS +from gwcs.spectroscopy import RefractedAngleSineModel from dkist.wcs.models import (AsymmetricMapping, Ravel, Unravel, VaryingCelestialTransform, VaryingCelestialTransform2D, VaryingCelestialTransform3D, - generate_celestial_transform, update_celestial_transform_parameters, + build_grating_spectral_transform, generate_celestial_transform, + update_celestial_transform_parameters, varying_celestial_transform_from_tables) @@ -52,6 +55,49 @@ def test_generate_celestial_unitless(): assert u.allclose(shift1.offset, 0) +def test_build_grating_spectral_transform() -> None: + header = { + "CTYPE1": "AWAV-GRA", + "CUNIT1": "nm", + "CRPIX1": 218, + "CRVAL1": 854.1738582455826, + "CDELT1": 0.0022975580183395555, + "PV1_0": 23000.0, + "PV1_1": 90, + "PV1_2": 65.696, + "PV1_3": 1.25, + "PV1_4": 1000.0, + "PV1_5": 1.5, + "PV1_6": 0.8, + } + transform = build_grating_spectral_transform( + reference_pixel=header["CRPIX1"] - 1, + reference_wavelength=header["CRVAL1"] * u.nm, + dispersion=header["CDELT1"] * u.nm / u.pix, + groove_density=header["PV1_0"] / u.m, + spectral_order=header["PV1_1"] * u.one, + incident_angle=header["PV1_2"] * u.deg, + refractive_index=header["PV1_3"] * u.one, + refractive_index_derivative=header["PV1_4"] / u.m, + out_of_plane_angle=header["PV1_5"] * u.deg, + camera_angle=header["PV1_6"] * u.deg, + ) + + pixels = np.array([0, 100, 217, 300, 511], dtype=float) + expected = WCS(header).spectral.pixel_to_world(pixels) + result = transform(pixels) + + assert isinstance(transform, CompoundModel) + # The transform should contain the gwcs refracted-angle model. + assert any( + isinstance(sm, RefractedAngleSineModel) + for sm in transform.traverse_postorder() + ) + np.testing.assert_allclose( + result.to_value(u.nm), expected.to_value(u.nm), rtol=1e-10, atol=1e-10 + ) + + def test_update_celestial(): trsfm = generate_celestial_transform( crpix=[0, 0] * u.pix,