Repository navigation
Planetary WCS implementation in Astropy #11170
Description
Description
In the framework of the EuroplanetH2020 some work has been done to define planetary standards for WCS.
Several research communities are involved in planetary coordinate standardization.
In particular Geologists and Remote Sensing specialists move from Earth to Planets using Geographical Information Systems (GIS) and proj based descriptions (if not ESRI!).
Astronomers move from Sky to Planets, using World Coordinate System and the Unit Sphere.
To improve interoperability between those two worlds we proposed a basic implementation of planetary WCS in GDAL (the FITS driver was already there and has been improved thereafter thanks to CNES funding). Planetary WCS are not yet a FITS standard but asking around (fitsbits list, IVOA workshops,...) everybody agrees on a pragmatic approach: once implemented in GDAL,
astropy and WCSlib they will just become a standard... :)
I like this approach... and after reading #5898, I'm reaching out to see if this is doable.
The novelties introduced with respect to the current WCS standard (which already support planetary coordinates in some way) are:
- the mapping of the
xysuffix forxyLN - xyLTCTYPE to the OBJECT keyword - the definition of a linear CTYPE
xyPX - xyPY(Projected X - Projected Y) which units are meters - the definition of the three keywords A_RADIUS, B_RADIUS, C_RADIUS to describe the shape of the reference system (DATUM)
- Ideally a switch between planetographic and planetocentric coordinates using the WCSNAME keyword value.
I'm an astropy user not a contributor (yet), but I will be happy to take care of this with some guidance if the team would like to consider it for integration.
Additional context
ping @thareUSGS and @Erard as co-authors of the paper and directly involved in the definition of planetary standards on GIS and IVOA side
Activity
cmarmo commentedon Mar 17, 2023
Dear maintainers, I may be able to be funded to work on this issue in the next months. I would be happy to offer some pull requests upstream if you are interested in.
I need a step-by-step approach to the issue (as suggested in #5898 (comment)). I have figured out some actions that would be useful to move forward, and I would be grateful for any comment and suggestion.
Is this an acceptable starting point?
Thank you so much for your attention.
mhvk commentedon Mar 17, 2023
This sounds great! While rereading some of the discussion, it also occurred to me that
sunpysurely has gone through something similar, so let me also ping @Cadair and @ayshih to see if they can provide some guidance of integrating coordinates on a body with WCS, etc.cmarmo commentedon Apr 21, 2023
Hello, thanks for the answer and sorry for the long silence: administrative details took more time than expected.
I have open a first pull request liberfa/erfa#94 to address item 1.
mhvk commentedon Apr 27, 2023
Further to my comment over at
erfa- liberfa/erfa#94 (comment) - my sense is that the first step from the astropy side is to enable a generalPlanetLocationclass that can take user-supplied data to provide -centric to -detic conversions.Let me also try pinging @Cadair and @ayshih again to see if they can point us to equivalent things for the Sun - it would be nice to have something consistent!
p.s. I think I would start not with
EarthLocation, but rather theBaseGeodeticRepresentationinastropy.coordinates.earth- I think that itself could start depending instead on a new base class that implements conversion toCartesianwith the general erfa routines. I think in that case, one likely wants a new subclass for every planet one wants to model.cmarmo commentedon Apr 27, 2023
@mhvk thanks for your answer.
I've already started working on a
BodyLocationclass as a generalization of theEarthLocation, but I was leveraging my pull request to erfa.I'm going to put on hold liberfa/erfa#94, and offer a draft PR on astropy following your suggestion.
But at some point the question about where the constants should be hosted will come up again.
Rocky bodies are really close to the Earth implementation in astropy, and less challenging than the Sun: while ad hoc spheroid radii and flattening are nice to have, most of the time rocky planets surface projections will make use of standard spheroids.
mhvk commentedon Apr 28, 2023
@cmarmo - sounds good! I must say I also think that once we have
BodyLocation, it makes little sense not to have the appropriate constants as well - for where to put those, it does seem mostly a choice betweenERFAandastropy.cmarmo commentedon May 11, 2023
As #14763 seems on its way in, next step would be the implementation of
BodyLocationitself as suggested here.While Earth standards about latitudes and longitudes are quite clear and homogeneous, in planetary applications multiple conventions are used.
The Gazetteer of Planetary Nomenclature provides all the possible conventions for each planetary feature (see eg for Abalos Mensa on Mars).
My plan was to add new parameters (east or west positive, odetic-ocenter,...) to the
BodyLocationmethods, so that the user may chose the convention for the output.EarthLocationwill inherit the methods with default parameters, so its behavior will not change.Another possibility could be to keep
BodyLocationsimple and then create aSolarSystemBodyLocationwhich add the new parameters.I have a slight preference for the first solution: please let me know if this is ok the maintainers.
Thanks for your attention.
mhvk commentedon May 11, 2023
@cmarmo - I think those changes are really differences in representation, and are most simply implemented as such. Certainly, using planetocentric latitudes is, as it needs a different conversion formula to cartesian - since
erfadoes not provide that, I guess we'd have to provide this ourselves.In principle, a different orientation (east or west positive) effectively implies different transformation formulae and thus would seem most logically implemented as a different representation (just like we have
PhysicsSphericalRepresentation, which differs fromSphericalRepresentationonly in swapping in latitude for co-latitude. For the differentwrap_angle, this should really not be necessary, but it still would seem most logically be done at the representation level. Indeed, this is a long-standing request - see #11767.Now for the bigger picture: I don't think
EarthLocationis necessarily the best model to use here. There is a long-standing issue #4261 about making it more like aBaseFrame. Whether that is best, I'm not sure, but it would seem more logical to not convert coordinates toCartesianupon creation, but rather as needed, i.e., be much more like aBaseFramewith itsrepresent_asmethod. Also, I think there really is little reason for it to be aQuantitysubclass - most of theQuantitymethods do not make much sense for it!So, I would suggest taking this as an opportunity to think through what a
PlanetLocationclass really needs to provide, and how it should allow transformations to other frames, etc. Specifically, it may be good to start with thinking about implementing the WCS transformations you mentioned at the head of this issue! What would they need?cmarmo commentedon May 11, 2023
The current funding should lead to the definition of a "body-fixed" coordinate system frame which I would like to use when reading or writing planetary FITS images. Ideally, I would like to instantiate body-fixed frames reading the RADIUS_? in the FITS header, the choice between planetodetic or planetocentric representation would be also determined by a specific keyword.
Also If I understand correctly a
BodyLocationwill prepare the transformation between the body-fixed frame and something like a planetary Alt-Az: this was my understanding for what is currently done for the Earth, that's why I am trying so hard in that direction... :)Happy to follow any other path as long as it'll lead me at the right place.
Have I answered to your question? Or at least facilitated a solution?
mhvk commentedon May 11, 2023
The first item you mention, of having a nice coordinate system which can be used to read/write FITS images of planetary bodies using various coordinate systems, seems the more pressing and perhaps also the easier, since surely this is a solved problem for the Sun! Indeed, certainly the sunpy documentation suggests so: https://docs.sunpy.org/en/stable/guide/coordinates.html
But obviously what you are hoping to do is to implement the basics needed to generalize such coordinates to other bodies. I think this would need to be some type of mash-up of how astropy converts WCS to sky images and how sunpy creates its
Mapinstances. I'm not very familiar with either, so perhaps best is to ping @Cadair and @ayshih to ask if they can give pointers - indeed, maybe forsunpyit would be helpful if some of there machinery (and overrides) could move upstream to astropy...mhvk commentedon May 11, 2023
Your second aim we have thought a bit about, but right now there is no clear solution. Apart from #4261, there is a lot more discussion in #10996 about how we can generalize the
locationattribute in coordinates away fromEarthLocation(which is what would be needed for something like a Mars AltAz). From that issue, I think the best approach would again be to follow sunpy and introduce a more generalobserverattribute; see #10996 (comment) - in many ways, this is like a regular coordinate frame, but I think it somehow keeps track of what it is fixed to. The discussion in #10996 was really quite good, and I haven't really caught back up to it...cmarmo commentedon May 12, 2023
Thanks @mhvk for your comments and the pointers: they are really helpful!
I have open #14820 as draft pull request, to move on on the implementation of the body-fixed reference system.
This is still work-in-progress and tests do not pass.
It is more to show that I have been able to progress thanks to your follow-up... :)
Obviously, any comment, suggestion, improvement is welcome, from anyone interested in the topic.
cmarmo commentedon May 18, 2023
Dear maintainers,
it is my understanding that sunpy makes use of the auxiliary information available in the auxprm structure from WCSLIB.
Its documentation states that "It is anticipated that this struct will expand in future to accomodate additional parameters."
Right now it contains parameters related to the Sun and its observation, eg, the reference radius.
I have been in contact with Mark Calabretta in the past and he was open to suggestions about planetary WCS implementation. Wouldn't be reasonable to have there the radii of the planetary body and the position of the observer?
mhvk commentedon May 18, 2023
Certainly seems reasonable, but let me ping the WCS maintainers, @mcara and @nden, as well as again @Cadair and @ayshih.
cmarmo commentedon May 19, 2023
While experimenting with planetary fits images working on the WCS extraction, I realized that I really needed some planetary options in the
BaseGeodeticRepresentation, in order to efficiently move forward.Following @mhvk suggestion, I have open #14851. Thanks for considering it.
cmarmo commentedon May 24, 2023
I've taken the liberty of interpreting the 👍 by @Cadair as a green light to contact Mark Calabretta: I'm preparing the list of keywords and corresponding additional parameters to be added to the auxiliary WCS structure.
cmarmo commentedon May 25, 2023
For those interested, the list of additional parameters, together with the related Planetary Data System (PDS) version 3 and 4 is available here.
nden commentedon Jun 3, 2023
This all looks reasonable to me. Sorry for the late response.
cmarmo commentedon Jun 24, 2023
I have marked #14820 ready for review: the shape of the spheroid is still missing.
Mark Calabretta already added the keywords described here in a beta release of wcslib. I'm planning to implement the reading of the shape of the spheroid in a subsequent pull request (as I am still having issues in the compilation of astropy with an external wcslib...). Thanks for your patience.
cmarmo commentedon Jul 6, 2023
Linking here #14990 as part of the same path towards reference frames for planetary bodies.
mhvk commentedon Sep 28, 2023
@cmarmo - I think we managed to do everything for this issue, right? If so, let's close it! (If not, list what still needs to be done here or in a new issue...)
cmarmo commentedon Sep 28, 2023
There are two "nice-to-have"s related to this issue:
The first should be quite easy, just needs some time allocation.
The second is more challenging.
I feel better when I close issues... :)
Perhaps I may open two issues about the previous.
Related to #5898 , transformations between reference frames would be wonderful to have ... but they need some funding... :)
mhvk commentedon Oct 4, 2023
@cmarmo - belatedly: I agree with closing issues, and think the two leftovers are better as new issues, as they are quite concrete and will not be easily found here. So, if you can raise those issues, that would be great! And then close this one, obviously!
cmarmo commentedon Oct 5, 2023
I have opened #15437 and #15438 to keep track of the remaining issues.
Hence, I'm closing this one.