Planetary Reference Frames¶

Skyfield has limited support for planetary reference frames as defined in Jet Propulsion Lab data files. It supports:

Loading .tf and .tpc text files defining a wide array of planetary constants.
Loading .bcp binary arrays that predict the orientation in space of Solar System bodies over a range of dates.
Given a latitude and longitude on a Solar System body like the Moon that is (a) spherical and (b) whose orientation is defined by a data series in a .bcp file, computing that surface position’s location in space — allowing observations either from the point of view of that surface position, or observations of that surface position for other observers in the Solar System.

This leaves several features of such files still unsupported, though. Skyfield:

Cannot yet support a planetary reference frame whose orientation is specified as a numeric series in a text file.
Cannot yet compute latitude and longitude positions on the surface of non-spherical objects, whose three axes in the text file are not the same length.

It is expected that support for these features will be added someday, making Skyfield’s support for planetary constants complete.

Observing a Moon location¶

Here is how you would load up enough data to predict where in the sky you would point a telescope to see a particular latitude and longitude on the Moon — in this example, the famous Aristarchus crater.

from skyfield.api import PlanetaryConstants, load

ts = load.timescale()
t = ts.utc(2019, 12, 20, 11, 5)

eph = load('de421.bsp')
earth, moon = eph['earth'], eph['moon']

pc = PlanetaryConstants()
pc.read_text(load('moon_080317.tf'))
pc.read_text(load('pck00008.tpc'))
pc.read_binary(load('moon_pa_de421_1900-2050.bpc'))

frame = pc.build_frame_named('MOON_ME_DE421')
aristarchus = moon + pc.build_latlon_degrees(frame, 26.3, -46.8)

apparent = earth.at(t).observe(aristarchus).apparent()
ra, dec, distance = apparent.radec(epoch='date')
print(ra)
print(dec)

13h 03m 22.96s
-00deg 55' 27.3"

If your Moon location has a nonzero elevation placing it above or below the Moon’s “sea level”, you can provide build_latlon_degrees() with an extra elevation_m argument.

Observing from a Moon location¶

You can also ask Skyfield: where in the sky would an astronaut standing on the Moon look to see a particular Solar System body? The answer can be provided in either right ascension and declination coordinates against the background of stars, or in altitude and azimuth measured against the astronaut’s horizon.

from skyfield.api import PlanetaryConstants, load

ts = load.timescale()
t = ts.utc(2019, 12, 20, 11, 5)

eph = load('de421.bsp')
earth, moon = eph['earth'], eph['moon']

pc = PlanetaryConstants()
pc.read_text(load('moon_080317.tf'))
pc.read_text(load('pck00008.tpc'))
pc.read_binary(load('moon_pa_de421_1900-2050.bpc'))

frame = pc.build_frame_named('MOON_ME_DE421')
aristarchus = moon + pc.build_latlon_degrees(frame, 26.3, -46.8)

apparent = aristarchus.at(t).observe(earth).apparent()
ra, dec, distance = apparent.radec(epoch='date')
print(ra)
print(dec)

alt, az, distance = apparent.altaz()
print(alt, 'degrees above the horizon')
print(az, 'degrees around the horizon from north')

01h 03m 22.96s
+00deg 55' 27.3"
32deg 27' 09.7" degrees above the horizon
118deg 12' 55.9" degrees around the horizon from north

Computing the sub-solar point on the Moon¶

This works not only for the Sun, but for any target body. You can learn the Moon latitude and longitude directly beneath the body by observing the target from the Moon’s center and then asking the lunar reference frame for the latitude and longitude.

sun = eph['Sun']

p = moon.at(t).observe(sun).apparent()
lat, lon, distance = p.frame_latlon(frame)
lon_degrees = (lon.degrees + 180.0) % 360.0 - 180.0
print('Sub-solar latitude: {:.1f} degrees'.format(lat.degrees))
print('Sub-solar longitude: {:.1f} degrees'.format(lon_degrees))

Sub-solar latitude: 0.3 degrees
Sub-solar longitude: -104.9 degrees

Computing lunar libration¶

The Moon’s libration is expressed as the latitude and longitude of the Moon location that is currently nearest the Earth. The convention seems to be that the simple geometric difference between the Earth’s and Moon’s positions are used, rather than the light-delayed position. Thus:

p = (earth - moon).at(t)
lat, lon, distance = p.frame_latlon(frame)
lon_degrees = (lon.degrees + 180.0) % 360.0 - 180.0
print('Libration in latitude: {:.3f}'.format(lat.degrees))
print('Libration in longitude: {:.3f}'.format(lon_degrees))

Libration in latitude: -6.749
Libration in longitude: 1.520

The only subtlety is that the libration longitude is not expressed as a number between 0° and 360°, as would be more usual for longitude, but instead as an offset positive or negative from zero, which the above code accomplishes with some quick subtraction and modulo.

Computing a raw rotation matrix¶

If you are directly manipulating vectors, you might simply want Skyfield to compute the NumPy rotation matrix for rotating vectors from the ICRF into the frame of reference of the Solar System body. The frame object returned above can return these matrices directly. If given a single time t, the result will be a simple 3×3 matrix.

from skyfield.api import PlanetaryConstants, load

ts = load.timescale()
t = ts.utc(2019, 12, 20, 11, 5)

pc = PlanetaryConstants()
pc.read_text(load('moon_080317.tf'))
pc.read_binary(load('moon_pa_de421_1900-2050.bpc'))

frame = pc.build_frame_named('MOON_ME_DE421')

R = frame.rotation_at(t)
print(R.shape)

(3, 3)

An array of times, by contrast, will return an array of matrices whose last dimension is as deep as the time vector is long.

t = ts.utc(2019, 12, 20, 11, range(5, 15))
R = frame.rotation_at(t)
print(t.shape)
print(R.shape)

(10,)
(3, 3, 10)

The transpose R.T of the rotation matrix can be used to rotate vectors that are already in the reference frame of the body back into a standard ICRF vector.