Implement option in WCS class to preserve the original units - #18338

Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for rebase and squash.
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the bot.
• Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) before merge.

👋 Thank you for your draft pull request! Do you know that you can use [ci skip] or [skip ci] in your commit messages to skip running continuous integration tests until you are ready?

I'm seeing a test failure in test_wtbarr_print which is a test unrelated to the changes here but what is happening is that the print_contents() call there includes some output from the new test I added in this PR, so somehow something isn't getting cleared in memory.

To illustrate why this PR is lovely, below is a WCSAxes plot with basically everything default:

Before:

After:

@mcara - are you happy with this in general? If so, I will carry out extensive testing in different packages and can then write more tests and some documentation.

Since it was in the draft mode I thought it was not yet ready for review. I'll take a look.

@mcara - yes just to be clear it works, it's just in draft mode because there are no docs yet.

@astrofrog While I review this, may I ask whether you had contacted Mark to see if this could be fixed in WCSLIB? It seems it would be a much more natural place for the fix as it would avoid unnecessary copying of Wcsprm structure. It seems unnatural/unexpected for WCSLIB to (kind of) arbitrarily convert units in the Wcsprm structure.

Section 2.3 of Paper I says:

2.3. Keyword value units
The original FITS paper (Wells et al. 1981) assumed that the units along each axis could be implied simply by the contents of the CTYPEi keyword. This has not turned out to be true in general. Therefore, we propose adding a new indexed, keyword CUNITi (character-valued) with which the units of CRVALi and CDELTi may be specified. Restrictions on the nature and range of units, if any, will be determined by the agreements applying to the specific axes. If they are not so limited, units should conform with the recommendations of the IAU Style Manual (McNally 1988).

Although in Paper II (section 1) it says:

Wells et al. (1981) established that all angles in FITS were to be measured in degrees and this has been entrenched by the AIPS convention and confirmed in the IAU-endorsed FITS standard (Hanisch et al. 2001). Paper I introduced the CUNITia keyword to define the units of CRVALia and CDELTia. Accordingly, we require CUNITia = 'deg' for the celestial CRVALia and CDELTia, whether given explicitly or not. Consequently, the (x; y) coordinates in the plane of projection are measured in degrees. For consistency, we use degree measure for native and celestial spherical coordinates and for all other angular measures in this paper.

This is confusing (to me). So I looked at the Hanisch et al. 2001 paper, section 5.3, and it says:

The units of all FITS header keyword values, with the exception of measurements of angles, should conform with the recommendations in the IAU Style Manual (McNally 1988). For angular measurements given as floating point values and specified with reserved keywords, degrees are the recommended units** (with the units, if specified, given as ’deg’).

recommended != required. I think it is worth asking why WCSLIB chose to behave like that and whether it is a bug and regardless of bug or not - can't this be handled there?

While I think it is worth contacting Mark Calabretta about this, I think we can go ahead and implement this here.

I also wonder how this works with the -TAB non-linear algorithm (see section 6 in Paper III - spectral). Maybe add an example using -TAB?

recommended != required. I think it is worth asking why WCSLIB chose to behave like that and whether it is a bug and regardless of bug or not - can't this be handled there?

I was in touch with Mark last year about this and this is what he said:

Correct - internally, WCSLIB only deals with SI units plus degrees. It would overly complicate matters to try to do otherwise.

Applications programmers see it when composing world coordinates to be transformed to pixel coordinates, or in the result of transforming pixel to world coordinates - I believe the documentation is clear enough on this. However, I feel applications programmers could protect higher-level users (i.e. the users of the applications) by doing unit conversions for them using the tools provided, namely wcsunitse().

This PR does pretty much that - we actually use wcsunitse in this PR to determine the correct scaling of units etc. I have followed up with Mark to see if he agrees that the approach in this PR is the right one, but at least from what he said before I don't think he would envisage changing WCSLIB.

In any case, thanks for the review! I will try and see if I can address these and wrap it up today or tomorrow.

What is the performance cost for this for the case of a single coordinate transformation (a pair of coordinates) and for a large array of coordinates?

Using the 5-D WCS from the test added in this PR (wcs_default is without preserve_units=True, and wcs_preserve is with), and doing a lot of repeated measurements, I get:

Before this PR

wcs_default.all_pix2world(2, 3, 8, 1, 2, 0)        0.014 +/- 0.002 ms
wcs_default.all_pix2world(pixel_coords, 0)               1.012 +/- 0.076 ms
wcs_default.all_world2pix(2, 3, 8, 1, 2, 0)        0.028 +/- 0.008 ms
wcs_default.all_world2pix(world_coords_si, 0)               1.178 +/- 0.146 ms

With this PR

wcs_default.all_pix2world(2, 3, 8, 1, 2, 0)        0.016 +/- 0.001 ms
wcs_default.all_pix2world(pixel_coords, 0)         1.128 +/- 0.060 ms
wcs_default.all_world2pix(2, 3, 8, 1, 2, 0)        0.019 +/- 0.002 ms
wcs_default.all_world2pix(world_coords_si, 0)      1.208 +/- 0.028 ms
wcs_preserve.all_pix2world(2, 3, 8, 1, 2, 0)       0.029 +/- 0.003 ms
wcs_preserve.all_pix2world(pixel_coords, 0)        1.234 +/- 0.001 ms
wcs_preserve.all_world2pix(2, 3, 8, 1, 2, 0)       0.037 +/- 0.004 ms
wcs_preserve.all_world2pix(world_coords_native, 0) 0.271 +/- 0.028 ms

Time times for wcs_default are consistent within errors, so no performance impact for the default case.

With preserve units, it seems pix2world is a little slower for pix2world, which is to be expected, but I am a bit surprised by world2pix actually being faster. The results are the same between wcs_default and wcs_preserve, so this is very puzzling!

Script to reproduce: https://gist.github.com/astrofrog/a0c88718992029198b79245114469d81

@mcara - I spent some time today working on making sure we can assign units to WCS.wcs.cunit, and I have something that seems to work. I have added a few tests for different scenarios. I've also made it so that if preserve_units=True, and WCS.wcs.set() has been called, users can't modify individual units eg WCS.wcs.cunit[1] in-place since that is difficult to deal with otherwise, and better to just have an explicit error. This won't break any existing usage since it only applies to when preserve_units=True.

In any case, I was wondering if you could take this for another spin and see if you are happy with the behavior under different circumstances?

@mhvk thank you for taking a look at this! Regarding the implementation i would prefer for now to keep the code paths reasonably separate for preserve_units False or True - mainly to avoid introducing bugs for regular users of the default WCS class - but we could always try and tighten up the implementation once it's had a good amount of real world testing?

I would like to try and add a ScienceState to control the default behavior, but would prefer to keep this for a separate PR.

I would like to try and add a ScienceState to control the default behavior, but would prefer to keep this for a separate PR.

Certainly fine, though I was thinking more about a configuration item. But given that, perhaps all the more reason to discuss separately.

Thanks for the review @mcara! I'll go ahead and merge, and will try and address @mhvk's suggestion regarding configuration in a follow-up PR.

Thanks for doing this @astrofrog this is something which has been a pain point for a long time.

I dare say this might be a killer feature for the next release. 😆 😭

I'm sure if we turned it on by default nothing would break.

nothing would break

Is that sarcasm? 😉

I think it would be important to have a method/property on the WCS object like is_preserving_units for downstream code to test which behavior is the object using. I forgot to mention this when reviewing this. Would it be possible to add this in a follow-up PR, @astrofrog ?

Yes of course!

@mcara - just to check, would it be better to simply have a WCS.preserve_units property? The reason I ask is because often a pattern is that things passed as kwargs to classes can then be accessed as properties, so there is some symmetry:

>>> wcs = WCS(..., preserve_units=True)
>>> wcs.preserve_units
True

Or would you prefer a method for this?

@mhvk @Cadair if you have opinions feel free to chime in!

+1 property which can't be set.

nothing would break

Is that sarcasm? 😉

From me?! Never.

I like property

For the record, drizzlepac regression tests all passed: https://github.com/spacetelescope/RegressionTests/actions/runs/16503552056

With the default behaviour or opting into preserve units? 😉

With the default behaviour or opting into preserve units? 😉

Default as drizzlepac is not yet aware of preserve_units and likely it won't ever use this since, I believe, HST WCS for all instruments use SI (but I may be wrong).

If we get a global switch, it'd certainly be good to check to see that there's no unexpected behaviour, even with SI units.

I think a global switch is bound to have consequences on other packages. For example, a package that is expecting RA & DEC to be in degrees will not be aware that, due to a global switch, it now gets angles in arcsec (for example).

Yeah, for sure, but it would be good to experiment with finding out where those assumptions are. Including inside astropy itself.

@mcara @Cadair lets discuss in #18463 ?