Implement APE 14 low-level FITS-WCS API and high-level WCS class

[astrofrog]

About

APE 14 describes a common low-level API for inter-operability, and #7325 adds the abstract base class to astropy.wcs. The aim of this pull request is to add the implementation of the low-level API for FITS-WCS (aka astropy.wcs.WCS) as well as the high-level class.

Current status

Example of use with the current PR:

In [1]: from astropy.wcs import WCS

In [2]: wcs = WCS('L1448_13CO.fits')

In [3]: wcs.world_axis_physical_types
Out[3]: ['pos.eq.ra', 'pos.eq.dec', 'spect.dopplerVeloc.opt']

In [4]: wcs.world_axis_units
Out[4]: ['deg', 'deg', 'm s-1']

In [5]: wcs.pixel_to_world_values(0, 1, 2)
Out[5]: [array(51.73971735), array(30.30833358), array(-9826.59656305)]

In [6]: wcs.pixel_to_world(0, 1, 2)
Out[6]: 
[<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
     (51.73971735, 30.30833358)>, <Quantity -9826.59656305 m / s>]

In [7]: from astropy import units as u

In [8]: from astropy.coordinates import SkyCoord

In [9]: coord = SkyCoord.from_name('L1448')

In [10]: vel = 10 * u.km / u.s

In [11]: wcs.world_to_pixel(coord, vel)
Out[11]: [array(147.95422821), array(43.99084675), array(300.48718796)]

In [12]: wcs.world_to_pixel(vel, coord)
Out[12]: [array(147.95422821), array(43.99084675), array(300.48718796)]

Future work

• Rebase this once Added base class for low-level WCS API from APE 14 #7325 is merged

• I need to try and implement world_axis_object_components/classes for time and spectral axes. The latter requires a decision in terms of what classes to use (see Decisions required). This will be done in a future PR

• This needs tests. A lot of tests. In particular I think it would be good if people could contribute examples of headers and what the output should be for the different properties.

• A changelog entry of course :)

[astropy-bot]

Hi there @astrofrog 👋 - thanks for the pull request! I'm just a friendly 🤖 that checks for issues related to the changelog and making sure that this pull request is milestoned and labeled correctly. This is mainly intended for the maintainers, so if you are not a maintainer you can ignore this, and a maintainer will let you know if any action is required on your part 😃.

Everything looks good from my point of view! 👍

If there are any issues with this message, please report them here.

[Cadair]

would it be worth putting functools.LRUCache on this or something, so it's not always computed twice per use in the high level API?

[eteq]

OMG. @Cadair, you've just blown my mind by pointing out the existence of functools.lru_cache. Where has this been all my life??

[astrofrog]

Good idea! Though the issue is that WCS objects are mutable, and it looks like the hash returned for the object doesn't change when the content of the WCS does.

[astrofrog]

I think I'd like to delay this to another PR because of the subtleties of mutable WCS

[astrofrog]

I added caching here, but not using LRUCache

[Cadair]

This shouldn't be hardcoded to degree right? If there is no projection then WCS will return in native units and not in deg. Unfortunately the only way I have found to reliably determine the output units is post-transform (https://github.com/sunpy/sunpy/blob/master/sunpy/map/mapbase.py#L871)

[astrofrog]

Will has_celestial be True in that case though? My understanding was that if WCS interprets the projection as celestial (which is what this branch is) then the units are degrees.

[eteq]

A few items while doing some experimentation:

• rename _wcs to low_level_wcs astrofrog/astropy#73
• add check that input to HighLevelWCS is a low level WCS astrofrog/astropy#74

And a bigger question for @astrofrog: is it possible to make the WCS object itself low-level API-compliant? I thought that was the goal, but maybe I missed somewhere that this was being prevented somehow.

Right now it's a lot of layers of indirection, so I'm pretty sure most users who actually understand WCS's would say "why would I bother, it's too hard to take control when I need it". So changing it to just one layer I think reaps a lot of rewards.

Put another way: right now to access the "real" fits wcs object, the user has to do highlevel._wcs._wcs.<whatever>, which is two layers of private - very undesirable if the user wants to modify the wcs, look at some of its state, etc. So instead, I think they should be able to do highlevel.low_level_wcs.<whatever>

[astrofrog]

@nden - I've pushed some fixes - does it seem ok now? Note that I've added a new property has_distortion on the WCS object to centralize the code to check if distortions are present.

[astrofrog]

I've also added more tests so the coverage of the astropy.wcs.wcsapi module is now 96%, so I think we could consider merging this soon as I think it should be reasonably stable.

First though, one question - are we comfortable adding the new API to the WCS object, or do we want to make that opt-in for one version? (it's apparently possible/easy to add inheritance to a class once it's initialized). So that would take the form:

wcs = WCS(...)
wcs.enable_ape14_api()

or something like that. Or shall we just go all in and have it be available by default straight away? (but not necessarily document it in the narrative docs yet)

The main reason I'm suggesting this is that some aspects of the API are not stable - e.g. pixel_to_world will eventually return Time and Frequency etc. objects instead of Quantity.

[astrofrog]

Does anyone understand why docstrings are not inherited from parent classes in the API docs?

https://11590-2081289-gh.circle-artifacts.com/0/home/circleci/project/docs/_build/html/api/astropy.wcs.wcsapi.fitswcs.FITSWCSAPIMixin.html#astropy.wcs.wcsapi.fitswcs.FITSWCSAPIMixin

I thought this was supposed to work, but apparently not.

[nden]

@astrofrog I'd say let's go with wcs.enable_ape14_api() first and start using it at least in a couple packages to see if anything needs to be changed. Can we merge this soon (now)?

Does nddata need any changes for this to be useful from other packages?

[eteq]

are we comfortable adding the new API to the WCS object, or do we want to make that opt-in for one version?

I think we should just add it, no opt-in needed. If we are concerned we should put a red-box warning in the docs that the API will change, but this won't present an immediate backward compatibility problem since the methods are all new, right?

I don't see the problem of eventually returning Frequency since that will be a Quantity subclass. I think everything will be except for Time, right? If so then maybe we raise a warning for times specifically?

[astrofrog]

@eteq - I added a warning. It turns out that WCSLIB 6.1 has just been released and includes support for Time, so we may be able to fix this sooner rather than later.

[eteq]

@astrofrog - just to clarify: you don't mean sooner as in in this PR, right?

[astrofrog]

@eteq - no I don't want to delay this PR. First we'll need to update WCSLIB in astropy which may not be trivial and might not happen for 3.1

Note to self - I need to implement some of @Cadair's suggestions higher up

[astrofrog]

I've added caching of _get_components_and_classes for performance

[eteq]

Just a couple minor things so I'm approving on the basis that basically it's all here.

One thing that's missing but I suspect might be planned as a second PR: a docs update that shows this functionality. I think we do need to at least demonstrate that this now works at a minimum, but I also think the following two things are probably needed as "detail" bits as well:

1. a discussion of the pixel conventions - i.e. the "array index" vs "pixel" meaning.
2. a bit on the custom_ctype_to_ucd_mapping decorator since it's conceptually rather complex

Do you plan to make that PR, @astrofrog? I get that it might be better to do that as a follow-on PR (and am fine with), but I think making those docs are a release blocker (which again I'm fine with as long as we can get it done soon after feature freeze).

[astrofrog]

@eteq - yes I can do the documentation improvements as a separate PR after this one. I think I've addressed your other requests.

[eteq]

The failing tests are unrelated so... Merging! 🎆 🍾