Adding functionality to represent non linear spectral dispersion - #713
Adding functionality to represent non linear spectral dispersion#713A-Derks wants to merge 35 commits into
Conversation
| ] | ||
|
|
||
|
|
||
| def _grating_equation_constants( |
| ) | ||
| alpha_in = m.Const1D(amplitude=adjusted_incident_angle_sine) | ||
| wavelength_from_grating = WavelengthFromGratingEquation( | ||
| groove_density=wavelength_parameter / spectral_order, |
There was a problem hiding this comment.
If I'm reading the math right, I think groove_density = wavelength_parameter / spectral_order only in the case where refractive_index_derivative = 0 and out_of_plan_angle = 0. Let's just set groove_density to the input grating_density, which is the same thing.
There was a problem hiding this comment.
But here's the part I don't get. If groove density is just G, then where is the denominator in Greisen eq 70? I'm looking at step 3 on page 756: "Compute the wavelength \lambda from \Gamma via Eq. 70" That equation has the two numerator parts, an adjusted sin \alpha and a sin \gamma, which we pass in the return of this function. The denominator in eq 70 is what is called wavelength_parameter above. Are we supposed to scale the sin(angles) before passing them to the model?
| adjusted_incident_angle_sine = ( | ||
| refractive_index - (refractive_index_derivative * reference_wavelength).decompose() | ||
| ) * np.sin(alpha) | ||
| wavelength_parameter = grism_parameter - refractive_index_derivative * np.sin(alpha) |
There was a problem hiding this comment.
Why is this variable called wavelength_parameter? When reading through these functions this name caused me a lot of confusion because I thought it had something to do with reference_wavelength or just wavelength in general. The units are 1 / m, which doesn't feel like a wavelength to me.
Based on my comment down in generate_grating_spectral_transform I'd argue we don't actually need this variable outside of this function. Given that, we could sidestep any naming issue by just removing this variable altogether and writing out angle_slope with the other variables directly.
There was a problem hiding this comment.
I agree with Arthur, but if you do want to name these parts of the grism equation (Greisen eq 70), then I would stick with names that imply units. You already do that with adjusted_incident_angle_sine, so you could name this part with something that includes ...inverse_wavelength.
| adjusted_incident_angle_sine, | ||
| wavelength_parameter, | ||
| reference_refracted_angle, | ||
| theta, |
There was a problem hiding this comment.
I don't think we need to return theta from this function. It's just a variable rename of the input parameter camera_angle, which (I think) just muddies the path of variables in generate_grating_spectral_transform. I.e., it makes it seem like theta is meaningfully different from camera_angle, which it isn't.
True, a .to(u.rad) has happened, but the return value is still a Quantity and so, e.g., np.sin will work correctly in any case.
| ( | ||
| adjusted_incident_angle_sine, | ||
| wavelength_parameter, | ||
| reference_refracted_angle, | ||
| theta, | ||
| angle_slope, | ||
| ) = _grating_equation_constants( | ||
| reference_wavelength=reference_wavelength, | ||
| grating_density=grating_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, | ||
| ) |
There was a problem hiding this comment.
It appears this function (_grating_equation_constants) exists to save a little bit of code repetition when computing some of the inputs to _refracted_angle_sine_model, alpha_in, and WavelengthFromGratingEquation, but, in doing so, it separates the inputs to generate_grating_spectral_transform from their actual usage and adds another layer of variable names and calculations. I think this added layer causes confusion (see my comments above) and obfuscates the flow of information through generate_grating_spectral_transform.
My suggestion would be to remove _grating_equation_constants entirely and rewrite _refracted_angle_sine_model to accept grating_density, spectral_order, refractive_index, and reference_wavelength. adjusted_incident_angle_sine could be computed directly in the function body.
| ) * np.sin(alpha) | ||
| wavelength_parameter = grism_parameter - refractive_index_derivative * np.sin(alpha) | ||
| reference_refracted_angle = np.arcsin( | ||
| (grism_parameter * reference_wavelength).decompose() - refractive_index * np.sin(alpha) |
There was a problem hiding this comment.
Is the .decompose() call necessary here? It looks like the results of this calculation stay in Quantity land for the duration of their life.
| epsilon = out_of_plane_angle.to(u.rad) | ||
| theta = camera_angle.to(u.rad) | ||
|
|
||
| grism_parameter = (grating_density * spectral_order) / np.cos(epsilon) |
There was a problem hiding this comment.
Not to be That Guy, but Greisen already defines something called the "grism parameter" in Eq. 75, which is different than the definition here.
There was a problem hiding this comment.
"grism constant" would be an OK alternative
| ) | ||
|
|
||
| pixels = np.array([0, 100, 217, 300, 511], dtype=float) | ||
| expected = WCS(header).all_pix2world(pixels, 0)[0] * u.m |
There was a problem hiding this comment.
As mentioned on the related dkist-inventory PR, a disadvantage of all_pix2world is that the results are unit-less and require us to apply the assumed-correct unit to its output. (It also has this weird signature and output that require non-obvious 0 and [0]).
Instead, I’d suggest using WCS(header).spectral.pixel_to_world(pixels), which will have proper units based on the CUNIT key. Then expected.to_value(u.nm) will work directly without needing to assume anything in the meantime.
There was a problem hiding this comment.
Sorry to keep bouncing parts of this code between here and dkist-inventory!
The thing this PR is missing is the asdf converter, schema and corresponding entries in the manifest, tests etc.
This is a thing, so this should be the last thing to do, but it will be needed to successfully test a roundtrip to ASDF file for this (which there should be a test of in dkist-inventory). We do have a helper script for this process which is here.
You should (for this PR) be able to run it like this:
$ pipx run tools/update_schema.py --manifest dkist-wcs refracted-angle-sine-modelHowever, it seems like this is bugged, so we'll have to fix it first. It's now fixed enough it runs, it should get you started.
| def _refracted_angle_sine_model( | ||
| reference_pixel, dispersion, reference_refracted_angle, theta, angle_slope | ||
| ): | ||
| @custom_model |
There was a problem hiding this comment.
Can you refactor this into class form, I think that with it like this it will be hard to add the required asdf schemas.
| assert u.allclose(shift1.offset, 0) | ||
|
|
||
|
|
||
| def test_generate_grating_spectral_transform(): |
There was a problem hiding this comment.
Can you also add tests for just the new model in isolation? I also suspect this test should probably be in dkist-inventory and not here.
| epsilon = out_of_plane_angle.to(u.rad) | ||
| theta = camera_angle.to(u.rad) | ||
|
|
||
| grism_parameter = (grating_density * spectral_order) / np.cos(epsilon) |
There was a problem hiding this comment.
"grism constant" would be an OK alternative
| adjusted_incident_angle_sine = ( | ||
| refractive_index - (refractive_index_derivative * reference_wavelength).decompose() | ||
| ) * np.sin(alpha) | ||
| wavelength_parameter = grism_parameter - refractive_index_derivative * np.sin(alpha) |
There was a problem hiding this comment.
I agree with Arthur, but if you do want to name these parts of the grism equation (Greisen eq 70), then I would stick with names that imply units. You already do that with adjusted_incident_angle_sine, so you could name this part with something that includes ...inverse_wavelength.
| reference_refracted_angle = np.arcsin( | ||
| (grism_parameter * reference_wavelength).decompose() - refractive_index * np.sin(alpha) | ||
| ) | ||
| angle_slope = wavelength_parameter / ( |
There was a problem hiding this comment.
I would like to rename this. If we use the name "grism parameter" for Γ (i.e. rename Gm/cosϵ) then this is the slope of the grism parameter wrt wavelength and wrt the distance along the detector (ξ in the paper) scaled. Γ itself is dimensionless. I think grism_parameter_slope would be better. Or maybe grism_parameter_per_wavelength, which would make its use in _refracted_angle_sine_model more clear.
| ): | ||
| alpha = incident_angle.to(u.rad) | ||
| epsilon = out_of_plane_angle.to(u.rad) | ||
| theta = camera_angle.to(u.rad) |
There was a problem hiding this comment.
If these are already Quantities, then I don't think we need to convert them to radians. In this function, we only ever do trig on them, not use them directly in algebra, so the conversion isn't necessary.
| adjusted_incident_angle_sine, | ||
| wavelength_parameter, | ||
| reference_refracted_angle, | ||
| theta, |
| output_angle = ( | ||
| np.arctan(-np.tan(theta) + wavelength_offset * angle_slope) | ||
| + reference_refracted_angle | ||
| + theta |
There was a problem hiding this comment.
Here is where theta aka camera_angle does need to be in radians. If it isn't passed through, the conversion should happen here.
| def refracted_angle_sine(pixel): | ||
| wavelength_offset = ((pixel - reference_pixel) * u.pix) * dispersion | ||
| output_angle = ( | ||
| np.arctan(-np.tan(theta) + wavelength_offset * angle_slope) |
There was a problem hiding this comment.
Renaming angle_slope to something like grism_parameter_per_wavelength (see comment here) then makes this line clearer about what it is doing.
| name="Spectral", | ||
| ) | ||
|
|
||
| return m.Mapping((0, 0)) | (alpha_in & refracted_angle_sine) | wavelength_from_grating |
There was a problem hiding this comment.
+1 to same name convention for both
| ] | ||
|
|
||
|
|
||
| class _refracted_angle_sine_model: |
There was a problem hiding this comment.
This should be refactored to be a subclass of Model like some of the other classes in this file, and not use @custom_model, like this example: https://docs.astropy.org/en/stable/modeling/new-model.html#a-full-example-of-a-linemodel
There was a problem hiding this comment.
Pull request overview
This PR adds support for representing FITS grating/grism spectral WCS (“-GRA/-GRI”) as an Astropy modeling transform, enabling non-linear spectral dispersion handling in dkist.wcs.
Changes:
- Added a grating/grism spectral transform builder (
build_grating_spectral_transform) plus a public compatibility wrapper (generate_grating_spectral_transform). - Introduced an ASDF-stable
refracted_angle_sine_modelcustom model to preserve a stable import path for serialization. - Added a regression test that compares the generated transform against
astropy.wcs.WCS(...).spectral.pixel_to_world().
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
dkist/wcs/models.py |
Adds the grating/grism spectral transform builder + ASDF-stable refracted-angle model; wires in gwcs spectroscopy model. |
dkist/wcs/tests/test_models.py |
Adds a test validating the new spectral transform matches Astropy WCS for a representative header. |
changelog/713.feature.rst |
Documents the addition of generate_grating_spectral_transform() as a compatibility wrapper. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Merging this PR will degrade performance by 6.94%
Warning Please fix the performance issues or acknowledge them on CodSpeed. Performance Changes
Tip Investigate this regression by commenting Comparing Footnotes |
…dated gwcs grating equation
… a thin refracted_angle_sine_model subclass in here
… a thin refracted_angle_sine_model subclass in here
…r_wavelength. Theseare now computed in gwcs
…which inputs have default values, removed generate_grating_spectral_transform(), adn updated changelog
…KISTDC/dkist into non_linear_spectral_dispersion
…which inputs have default values, removed generate_grating_spectral_transform(), adn updated changelog
eigenbrot
left a comment
There was a problem hiding this comment.
Approving, but the current PR description is out of date
…_pure_python.yml (#731) Bumps [OpenAstronomy/github-actions-workflows/.github/workflows/publish_pure_python.yml](https://github.com/openastronomy/github-actions-workflows) from 2.5.0 to 3.0.1. - [Release notes](https://github.com/openastronomy/github-actions-workflows/releases) - [Commits](OpenAstronomy/github-actions-workflows@a138926...e5af21e) --- updated-dependencies: - dependency-name: OpenAstronomy/github-actions-workflows/.github/workflows/publish_pure_python.yml dependency-version: 3.0.1 dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] <support@github.com> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Drew Leonard <andy.j.leonard@gmail.com>
…ts to load_dataset, which silently dropped it and only honored it for Path input (#737) Co-authored-by: Drew Leonard <andy.j.leonard@gmail.com>
Bumps [actions/checkout](https://github.com/actions/checkout) from 6 to 7. - [Release notes](https://github.com/actions/checkout/releases) - [Commits](actions/checkout@v6...v7) --- updated-dependencies: - dependency-name: actions/checkout dependency-version: '7' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] <support@github.com> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Drew Leonard <andy.j.leonard@gmail.com>
Bumps [CodspeedHQ/action](https://github.com/codspeedhq/action) from 4.17.0 to 4.18.1. - [Release notes](https://github.com/codspeedhq/action/releases) - [Changelog](https://github.com/CodSpeedHQ/action/blob/main/CHANGELOG.md) - [Commits](CodSpeedHQ/action@9d332c4...a4a36bb) --- updated-dependencies: - dependency-name: CodspeedHQ/action dependency-version: 4.18.1 dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] <support@github.com> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Drew Leonard <andy.j.leonard@gmail.com>
updates: - [github.com/zizmorcore/zizmor-pre-commit: v1.25.2 → v1.26.1](zizmorcore/zizmor-pre-commit@v1.25.2...v1.26.1) - [github.com/astral-sh/ruff-pre-commit: v0.15.15 → v0.15.20](astral-sh/ruff-pre-commit@v0.15.15...v0.15.20) - [github.com/PyCQA/isort: 8.0.1 → 9.0.0a3](PyCQA/isort@8.0.1...9.0.0a3) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Drew Leonard <andy.j.leonard@gmail.com>
# Conflicts: # .github/workflows/build_workshop.yml # .github/workflows/ci.yml # .github/workflows/codspeed-walltime.yml # .github/workflows/codspeed.yml # .github/workflows/dkist_downstreams.yml # .github/workflows/prepare-release.yml # .github/workflows/sub_package_update.yml # .pre-commit-config.yaml
Background:
DKIST instruments (Cryo SP, DL-NIRSP, ViSP) produce FITS headers with CTYPE values of AWAV-GRA, indicating a non-linear pixel-to-wavelength transform based on the grating/grism equation described in Greisen et al. (2006, FITS WCS Paper 3, sections 5.1.2 and 5.2). The PV{n}_m keywords encode the physical grating parameters.
The gwcs library (via gwcs#745) provides two models for evaluating the FITS grism transform: RefractedAngleSineModel, which computes the pixel-dependent refracted-angle sine from the detector geometry and grating parameters, and WavelengthFromGrismEquation, which converts the incident-angle sine and refracted-angle sine into wavelength using the FITS grism formalism. WavelengthFromGratingEquation remains the classic two-input grating-equation model and is not used by this PR. This PR provides the builder function that wires the FITS grism models together using the parameters extracted from the FITS header.
What changed:
dkist/wcs/models.py
build_grating_spectral_transform — the core implementation. Given the full set of grating/grism parameters, it:
All four optional grism parameters (refractive_index, refractive_index_derivative, out_of_plane_angle, camera_angle) default to their identity values, so the transform degrades to a standard grating equation when they are absent from the header.
RefractedAngleSineModel and WavelengthFromGrismEquation are defined in gwcs.spectroscopy (added via gwcs#745), which also provides the ASDF converter that serializes them. This PR imports both directly from gwcs rather than defining or serializing a local model in dkist.
Notes:
The CI test failures are due to the gwcs dependency: RefractedAngleSineModel and WavelengthFromGratingEquation do not exist in any released gwcs version yet. Tests will pass once gwcs#745 merges and the dependency is pinned or updated.
The CodSpeed regression (-6.9% on test_plot_dataset) is unrelated to the spectral transform code — it affects a plotting benchmark and is likely noise or a pre-existing issue.