API Document for WCSAxes

[astrofrog]

A number of packages, such as pywcsgrid2, APLpy, and kapteyn, have developed ways of making plots in world coordinate systems (WCS) in Matplotlib, and since this is a common issue in Astronomy, we are interested in implementing the general problem of plotting WCS axes in the Astropy core, so that we can avoid duplicating this in other packages. The present document describes the proposed API.

The aim is not to provide the simplest possible user-friendly API to the user, but to design a flexible API consistent with the matplotlib API upon which user-friendly tools can be developed.

This API was written by @astrofrog and @leejjoon, and benefited from feedback from @ebressert and @ChrisBeaumont. It is written in markdown, since we wanted to include images. The latest rendered version of the API document (with images) can be seen here:

https://github.com/astrofrog/astropy-api/blob/7fba6015426e1604440b85f2cf0964ff7b6c9d35/wcs/wcs_api.md

At this time, we would be very interested in getting feedback from other developers (or users) regarding this API. In some places, the API may not be optimal, so please let us know if you have any suggestions! Also, if you see anything missing, please let us know (and offer suggestions if possible!).

(edited to update link)

[keflavich]

What about ax[other_wcs].pcolormesh and its ilk?

[astrofrog]

I think that would work because it works by plotting polygons. It's just that imshow doesn't support transformations.

[leejjoon]

In the case pywcsgrid2, ax[other_wcs] returns a new axes but with different transform. So, everything vector-based is supposed to work. Showing an image is quite different though.

pcolormesh, which is vector based, should work. But pcolormesh is quite slow unless you use agg backends.

Matplotlib now supports affine-transformed images (not for every backend but most of them), so we can support imshow by locally approximating the transform as an affine-transform. And my experience is that this is enough for most cases, but will not applicable for things like allsky plot, for example. Or we can regrid the image on the fly of course.

[mdboom]

Arbitrary image resampling (a la drizzle) has been on the TODO list for matplotlib for a long time. I'd love to find the time to do it (or find someone else's time to do it). In the meantime, I think, as @leejjoon points out, our best bets are to approximate with affines or use pcolormesh if the transform is too non-linear.

Oh -- and bonus points for providing a method that will use one or the other depending on whether the error rate of using imshow is above some threshold.

[pelson]

On this front, cartopy can handle inter-coordinate reference system transformations of image via a spherical nearest neighbour interpolation scheme. On top of this, we have Frankenstein-ed the pcolormesh method to capture any non-quadilateral (or wrapping) data to allow plotting of data which might very well cross the dateline without a hitch.

I can put together some examples (sadly documentation is still lacking for cartopy) if interested.

[keflavich]

Overall looks great. I left some inline comments. One example request would be something like the figures shown here: http://stackoverflow.com/questions/13633254/offset-polar-axes-grid-display-coordinate-display, i.e. a face-on view of the Galaxy. This isn't a native WCS format, so maybe it doesn't make sense to include as a WCSaxes, but perhaps it does make sense to include as a more general transformation? @astrofrog, you make plots in these coordinates with hyperion, I think, so I suspect you already have this capability in mind?

[astrofrog]

I think the biggest remaining thing to decide relates to @mdboom's comment regarding the ax['fk5'] notation. We have several alternatives:

ax.scale('fk5').method(...)
ax.transform('fk5').method(...)
ax.frame('fk5').method(...)
ax.method(..., transform='fk5')
ax.fk5.method(...)

how do people feel about this? I kind of like the idea of using transform because then we aren't adding anything to matplotlib (except to automagically interpret 'fk5' as a transformation). It might be the hardest to implement, but maybe we should try that first before resorting to anything else? There is something to be said about the other syntaxes though, because they allow

ax_gal = ax.frame('gal')

which can then be used, and everything done on it will work in Galactic coordinates. However, maybe that's a little confusing, because by default, all ticks and labels would be hidden (since maybe they just requested to overplot some polygons, but they don't care about getting a coordinate overlay). Maybe better to use the transform keyword for that, and then if users want a coordinate overlay, they do

ax.get_coords_overlay('galactic')

Any thoughts?

Finally, the only other tricky bit is to do with labels, in that there are tick and axis labels. However, maybe we can simplify by saying that any reference to labels for the coordinate objects is to do with tick labels, and we don't actually touch anything to do with axis labels, since the standard API is fine for that?

[mdboom]

I think I prefer the transform as a keyword option, as that avoids adding anything new to the API that isn't in matplotlib, and having a new method like get_coords_overlay as you suggest to do the grid drawing. I think the ax.frame('fk5').method(...) syntax is potentially interesting, but it's a bit confusing and I haven't seen a killer use case for it that makes it compelling.

[astrofrog]

I decided to switch to transform= as the API is a lot cleaner that way, and more consistent with matplotlib.

For now, I decided to leave the slightly different methods for axis and tick labels, because I think it makes it much more readable, even if slightly inconsistent with matplotlib.

Finally, I changed the section for multi-dimensional WCS so that this is now something the WCS class has to implement.

I think that concludes the comments I had to address - I will send one last email to the mailing list before merging this, if there are no objections.

[astrofrog]

The latest rendered version can be seen here: https://github.com/astrofrog/astropy-api/blob/7fba6015426e1604440b85f2cf0964ff7b6c9d35/wcs/wcs_api.md

[leejjoon]

If we want to do something like transform="fk5", we need to change the internals of matplotlib. Or we need to override all the Axes methods.

Most easy way is to have an Axes method that returns a new transform, but again this introduces a new method, and not very elegant I guess.

Personally, I am not a big fan of using this api. I think this, in some sense, abuses the transform keyword ("abuse" maybe too strong). But if others agree to go this way, that's completely fine with me. But I think we should discuss how to implement this!

To me, if we have different transform (transData, more exactly), it is natural to introduce a new axes. This is somewhat similar to twinx and twiny in the vanilla matplotlib , although there are a bit of difference.

And, I think the api like below

ax_gal = ax.frame('gal')

have good advantage over the current api due to similar reasons that already mentioned by @astrofrog.

[leejjoon]

Are we going to put WCSAxes inside the wcs module? They are quite different in nature, and I am not sure if this is a best idea, although this not critical.

[astrofrog]

I think they belong together because one is really just a graphical representation of the other - and WCS has to be able to return a WCSAxes instance from the _as_mpl_axes method. In fact, users shouldn't even have to invoke WCSAxes directly, since they can do projection=WCS(...). But I'm interested to hear your opinion about why you feel they should belong in a different module?

[leejjoon]

Well, I think my view is that wcs.WCS is one implementation of the world coordinate system, and WCSAxes is one implementation of the graphical representation of the world coordinate system. Anyhow, this is nothing important.

[mdboom]

I think I'm with @leejjoon on this one, but even further ;) astropy.wcs.WCS is one implementation (wcslib) of one representation (FITS WCS) of WCS. In fact, once we can make such a change, I'd almost want to rename it to fits_wcs to help distinguish it from the more general WCS that are sure to come... I think I'd prefer they were independent models. I can envision this plotting work being extended in the future to have different transformation backends, the first being astropy.wcs, but many to come in the future.

[astrofrog]

Right, I should have said that when I was talking about the WCS class, it doesn't have to be FITS WCS only. In my mind, astropy.wcs was intended to contain any WCS system in the long term, so that's why I see WCSAxes living there. However, if you still think it would be better to have it elsewhere, we could have an Astropy sub-package dedicated to plotting, e.g. astropy.plotting or astropy.graphics where other Astronomy-specific visualizations could live? I don't think this is an important point for now though, as I think we'll start implementing it as an affiliated package anyway until it's ready for inclusion, at which point we can decide where it should live (and by then, maybe we'll have made a decision concerning the generalized WCS).

[eteq]

I agree with @mdboom and @leejjoon that that location implies a tie to the WCS object rather than to a future "generalized" WCS. The question is where the gWCS ends up - if we end up creating astropy.gwcs or something, it's a problem that WCSAxes is in wcs. But @astrofrog is right that we may want to end up putting the generalized WCS in astropy.wcs anyway.

So if we are sure we don't want the future gWCS stuff in wcs, than it makes sense for WCSAxes to be elsewhere. But if there's a chance of it being in wcs, I agree with @astrofrog that it's more expedient and understadable to put it here for now.

[astrofrog]

The thing is that the recommended usage for this functionality is actually going to be

 ax = fig.add_subplot(1, 1, 1, projection=wcs.WCS('image.fits'))

so it's actually the WCS class that's passed to matplotlib, and this calls WCS._as_mpl_axes which has to return an Axes class, so doesn't it make sense for the axes class to live alongside the WCS class? The point is that I think the WCSAxes class may even be 'private' and never accessible to the user, which is why I'm not so keen on creating a new sub-package for it.

Since it'll be a private class anyway, then it doesn't matter where it lives so we can just put it alongside WCS.

[eteq]

Well, I don't see why we would make it private. It will end up user-facing because if you do pyplot.gca() or similar that's what object you will get. So I think it should still be treated as a "public" class in terms of API decisions. (Also, I could certainly imagine someone wanting to subclass it in some circumstances)

Don't get me wrong, I agree with you that wcs is a sensible place right now... I just don't think we want to treat it as "subordinate" to the current WCS given that we'd rather it eventually be tied more closely to the generalized WCS. (I wish there was an emoji 🤷)

[astrofrog]

You are right that the class needs to be accessible directly too.

For now, I'll develop it alongside the current WCS, but of course we can rethink this when we come up with the generalized WCS.

[eteq]

👍

[leejjoon]

I think the tick and label related API is a bit ambiguous.

In one example, I see

glon.set_ticklabel_position('t')
glon.set_ticklabel_color('green')

Are you suggesting that we implement methods like these for all (or most) properties, e.g., fontsize, alpha, etc?

Maybe one option is to just adopt set_tick_params method as in matplotlib.

[astrofrog]

@leejjoon - ah, I did not know about set_tick_params - that makes sense!

[astrofrog]

Thinking about this more, I agree with @leejjoon that using transform is going to be hard, and especially in cases where we need to specify additional parameters, e.g. epoch and equinox. Maybe the following makes more sense after all:

ax.frame('fk5', equinox='J2010')

Furthermore, since we can show and hide the coordinate grid with:

ax.coords.set_visible(False/True)

then it makes sense that one could simply do

ax2 = ax.frame('fk5', equinox='J2010')
ax2.coords.set_visible(True)

to show an overlay, so we don't need a new method for that.

@mdboom @leejjoon - what do you think?

[leejjoon]

Both approaches have pros and cons. So, here is my suggestion. I think the first step toward this is to have an API such as WCSAxes.get_sky_transform which returns a matplotlib transform instance of given coordinate. It should take a string argument such as "fk5", "gal", etc. It also will be able to take a wcs.WCS instance.

tr = ax.get_sky_transform("fk5")

The returned transform will be essentially = "transform from fk5 to pixel coordinate of ax" + ax.transData.

When the argument is a WCS instance

tr = ax.get_sky_transform(wcs.WCS("another.fits"))

then the returned transform is = "transform from pixel coordinate of 'another.fits' to the pixel coordinate of ax" + ax.transData.

(1) Given that we have this method, we will be able to do something like below.

tr_fk5 = ax.get_sky_transform("fk5")
ax.plot(pos_ra, pos_dec, "o-", transform=tr_fk5)

(2) Once Matplotlib is modified to support such a usage like transform="fk5", it will be quite straightforward to adopt it.

(3) With a method like get_sky_transform implemented, it is also straightforward to create an axes with that transform.

So, we first implement (1). We then implement (3) for advanced users. And then implement (2) when matplotlib is ready. How others think?

Anyhow, the name of get_sky_transform may not be a good choice, and please propose a new name.

[astrofrog]

@leejjoon - I like this approach - very clean! I wonder whether get_sky_transform (or whatever it's called) could actually be a method of WCS? (since it's just basically extending the WCS transformation to a different system). What do you think?

[leejjoon]

With the current API I proposed, get_sky_transform is supposed to return an instance of matplotlib's transform class. So, I think it is better to be a method of WCSAxes. But, WCS itself may implement a similar method, and then WCSAxes.get_sky_transform just utilizes it. But I am not sure what the WCS's method should return in this case. Is there a general transform class within astropy?

[astrofrog]

@leejjoon - ah I see, you don't want to tie the WCS class to Matplotlib. In that case, what you suggested makes sense.

[astrofrog]

@leejjoon - I've implemented your suggestion, which I really like, although I renamed it get_transform instead of get_sky_transform because e.g. in the case of passing it a WCS object, there's no reason it has to be a sky projection.

[astrofrog]

@mdboom @leejjoon - do you have any final comments before I propose that this be adopted? Thanks for all your feedback!

[mdboom]

No additional comments from me.

[eteq]

"pyplot-style" to mean means the functional paradigm of using the pyplot functions, like matplotlib.pyplot.axes and matplotlib.pyplot.subplot . But below it looks like you're talking about the OO-paradigm of doing fig.add_suplot(...). Do you intend to add functions to support the functional interface?

[astrofrog]

My understanding is that once we implement WCS. _as_mpl_axes, we can call projection=WCS(...) in both the functional (matlab-style) interface, as well as the OO-like and full-OO API.

@mdboom - can you confirm whether this is the case?

[astrofrog]

By the way, if this is the case, I will clarify this in the document.

[astrofrog]

I clarified this in the document.

[eteq]

Aside from my comment above (and my response to one of the old comments further up, which may not require any action), this seems good to me!

[astrofrog]

I will now merge this - thank you to everyone for your feedback! Now comes the hard part, actually implementing it :) I will start co-ordinating this at some point over the next few weeks.

If you have any suggestions for improvements, please open a new issue, or even better, a pull request!