Planetary orientation models from text PCK file

[psbaltar]

This adds a PlanetaryOrientation class to planetarylib as mentioned in #796. The class is able to provide the planet's orientation (RA and DEC of the pole and the prime meridian), as well as the planet's shape (three radii for an ellipsoid).

The models are read from a text PCK file. There's notes about how binary PCK files should be preferred if available as they are more precise. In fact, the Earth orientation model in the text PCK file doesn't include nutation. From what I can tell, there's binary files only for Earth and Earth's moon.

from skyfield.api import load, PlanetaryConstants, PlanetaryOrientation

ts = load.timescale()
t = ts.utc(2020,1,1)

pc = PlanetaryConstants()
pc.read_text(load('pck00010.tpc'))

po = PlanetaryOrientation(pc, '399')
print(po.orientation_model, po.nutation_model, po.shape_model) # True False True

po.at(t)
print(po.orientation) # [-0.12819123924090417, 89.8886076127033, 2637009.9229112]
print(po.radii) # [6378.1366, 6378.1366, 6356.7519]

The models are described in detail in Archinal, B.A., A’Hearn, M.F., Bowell, E. et al. Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009. Celest Mech Dyn Astr 109, 101–135 (2011). https://doi.org/10.1007/s10569-010-9320-4 Also available here: https://aa.usno.navy.mil/downloads/reports/Archinaletal2011a.pdf

A SPICE PCK file that PlanetaryOrientation reads the models from is available here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc

[brandon-rhodes]

Thanks very much for putting together this code! Having Skyfield be able to compute and return these values will be very useful. I'll see if I can write up a quick code review this week; if not, then look for something from me by mid-October.

[psbaltar]

I made some updates to vectorize the calculations. I hadn't been using Skyfield's vectorization features because I wasn't very familiar with the topic. It's something I hadn't worked with before in numpy.

I also wanted to take the opportunity to thank you for this great library! I'm not a professional astronomer nor a professional programmer, so I've learned a lot while working with Skyfield.

[billyauhk]

I should point out that for Earth, Skyfield should already supports using Earth Orientation Parameters; for Moon, libration angles from an ephemeris should be used. The PCK is not necessary for Earth/Moon, and might not be the currently-accepted values.

Since WGCCRE writes a new report every 3 years, values from a PCK may not be updated to the current standard at https://astrogeology.usgs.gov/groups/IAU-WGCCRE The 2015 standard is here: https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf

(P.S. Not a Skyfield user on my own, but found this when I am trying to check if Skyfield could compute physical ephemerides of planets.)

Cross-ref #417 #709 as both will involve using WGCCRE values at some point.

[psbaltar]

Since WGCCRE writes a new report every 3 years, values from a PCK may not be updated to the current standard at https://astrogeology.usgs.gov/groups/IAU-WGCCRE The 2015 standard is here: https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf

I've been using pck00010.tcp for my own purposes (non-research) because it's the latest I can find from NAIF. It's dated 2011-10-21 and says it's based on the 2009 WGCCRE report. I had come across the 2015 report, but I haven't gone through the exercise of checking what changed. However, I did see that they removed the orientation model for Earth because "accuracy was poor, and the expressions failed near the fundamental epoch (J2000.0), yet they were sometimes used as a recommended model".

Since the PCK file is plain text, it wouldn't be hard for a user to make any necessary updates if they wanted to.

[billyauhk]

For Earth/Moon, the models are kind-of removed from the 2015 report in order to prevent people misusing them and treat it as the "correct" expressions. Precession-nutation models for Earth are already inside Skyfield but one may need some additional code to retrieve what you want (I don't know what you want at the very end).

Some variation of the code at https://github.com/skyfielders/python-skyfield/issues/46#issuecomment-95299079 will be useful in deriving a rotation matrix (and thus a milli-arcsecond precise orientation of the Earth), so adding new code with known and inferior precision and add confusion to users sounds wary to me.

@psbaltar I don't think "my own purposes (non-research)" is a responsible way to add code into a public repository, without highlighting the possible side-effect which may affect others. Known problems/deficiencies should at least be documented (this is why I added my opinion in the first place). I have no opinion on using a PCK for any other Solar System Bodies other than Earth/Moon, but please try to prevent user unknowingly using an inferior model when the correct one is already there.

P.S. I know there are some PCK containing post 2015 values floating around, but I don't have the copyright to release them. Try one's luck searching for "oriccre2015_20210225.tpc" or "pck00011.tpc" on the web.

[psbaltar]

@billyauhk Sorry for the misunderstanding. I thought that the warnings against using the PCK files were meant as a blanket statement to not use them because they may be out of date.

I wasn't trying to imply/recommend that the PCK data should be used for any particular cases. This is just provided as generic code that's capable of reading a model from a given file. While there are already better models for Earth/Moon, this is a way to get models for the other planets. The data file it is intended to read just happens to include Earth and Moon models as well, even though they aren't current/best. In fact, the PCK file itself cautions against using its Earth rotation model.

@brandon-rhodes would have to weigh in on how best to make the user aware of the limitations.

[brandon-rhodes]

@psbaltar — My apologies for taking so long to respond to this pull request! The last few months have been very busy, but I'm returning to my GitHub notifications at last.

First, I want to reassure you that I don't find this contribution at all irresponsible! It’s not as though you are adding new pages to the Skyfield documentation, telling users to ignore the normal wgs84.latlon() and ITRS techniques for Earth orientation and replacing them with a PCK file. When this feature is documented, we can simply highlight its use in things like the orientation of Mars and Jupiter and Saturn, and steer people away from old PCK files that might list the Earth.

Second, I wanted to ask if you're interested in my writing up a full code review—or whether I have delayed too long in answering, and by now you've lost interest and moved on to other projects? If you're still around, I would be happy to make comments on the code—just let me know either way, thanks!

[psbaltar]

@brandon-rhodes Thanks for getting back to me. A code review would be great. Planetary imaging season is over, so I've back burnered this. I do plan on circling back to it though.

Either way, I already have the code and would like to contribute to Skyfield since it's been so helpful!

[brandon-rhodes]

Thanks for the reply! The next couple of weeks will have me busy with preparations for PyCon, but I'll do the code review right after the conference if I don't get to it beforehand.