Kepler Orbits¶

Skyfield now offers basic support for computing the position of a comet or minor planet whose elliptical, parabolic, or hyperbolic orbit is provided as Kepler orbital elements.

Beware that the internal routines supporting Kepler orbits are rudimentary and subject to change — only the interface documented here is guaranteed to work in future Skyfield versions.

Skyfield loads orbital elements from text files using the Pandas library. Install it before trying any of the the examples below:

pip install pandas

Comets¶

The IAU Minor Planet Center distributes a CometEls.txt file of orbital elements for predicting comet positions. The file is plain text, so feel free to open it with a text editor to see the comets for which it offers orbital elements. To build a dataframe of comets:

from skyfield.api import load
from skyfield.data import mpc

with load.open(mpc.COMET_URL) as f:
    comets = mpc.load_comets_dataframe(f)

print(len(comets), 'comets loaded')

864 comets loaded

Since the comets file has no explicit expiration date, load.open() will only download the file once. Subsequent calls re-open the copy of the file already on your filesystem. To force a fresh download and receive updated orbits and new comets, pass reload=True.

To generate a comet’s position, first select its row from dataframe. There are several Pandas techniques for selecting rows, but most Skyfield users will simply index their dataframe by comet designation.

# Keep only the most recent orbit for each comet,
# and index by designation for fast lookup.
comets = (comets.sort_values('reference')
          .groupby('designation', as_index=False).last()
          .set_index('designation', drop=False))

# Sample lookups.
row = comets.loc['1P/Halley']
row = comets.loc['C/1995 O1 (Hale-Bopp)']

When computing the position of a comet from Earth, there is a complication: cometary orbits are not measured from the Solar System barycenter but are instead centered on the Sun. You will therefore need to add the barycenter→Sun vector to the Sun→comet vector to produce a position that you can pass to the observe() method, which always measures positions from the Solar System barycenter.

# Generating a position.

from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN

ts = load.timescale()
eph = load('de421.bsp')
sun, earth = eph['sun'], eph['earth']

comet = sun + mpc.comet_orbit(row, ts, GM_SUN)

t = ts.utc(2020, 5, 31)
ra, dec, distance = earth.at(t).observe(comet).radec()
print(ra)
print(dec)

23h 59m 16.85s
-84deg 46' 57.8"

Hopefully Skyfield will in the future support generating positions for whole arrays of comets in a single efficient operation, but for now your code should expect to operate on one comet at a time.

Minor Planets¶

There are nearly a million minor planets in the IAU Minor Planet Center’s database of orbital elements, thanks to the prodigious output of automated sky surveys over the past few decades.

The database can be downloaded as a single MPCORB — “Minor Planet Center orbits” — file that offers each minor planet’s orbital elements as plain text. But the raw file requires a bit of preprocessing before Skyfield is ready to load it:

The first 43 lines of the file are paragraphs that explain its contents, state the terms under which software programs may include the data, and provide links to documentation. Skyfield will need these non-data lines ignored or removed.
While an uncompressed version of the file is available for download, most users opt to download the 55 MB compressed version from the Minor Planet Center to save bandwidth and storage. Decompressing the full 190 MB of data stored inside can require more than 1 second of computing time depending on your platform and processing speed.
The complete catalog lists nearly 1 million objects, which can take several seconds to load and index.

For all of these reasons, it usually makes the most sense to download, uncompress, and filter the file before starting your application.

If your operating system provides tools for pattern matching, they might be the fastest tool for selecting specific orbits. Here’s how to extract the orbits for the first four asteroids to be discovered — (1) Ceres, (2) Pallas, (3) Juno, and (4) Vesta — on a Linux system:

zgrep -P '^(00001|00002|00003|00004) ' MPCORB.DAT.gz > MPCORB.excerpt.DAT

If your operating system lacks such tools, you can build them yourself using Python. Note that mass operations that Python implements in C, like reading an entire file’s contents with read() and scanning the full contents with a regular expression findall(), will be much faster than using a Python loop to read every line. Here’s an example script for performing the same search as the zgrep command shown above:

# mpc_make_excerpt.py

"""Search the MPCORB file for minor planets, given their packed designations."""

import argparse
import re
import sys
import zlib

from skyfield.api import load
from skyfield.data import mpc

def main(argv):
    parser = argparse.ArgumentParser(description='Grep MPCORB.DAT.gz')
    parser.add_argument('designations', nargs='+', help='packed designations'
                        ' of the minor planets whose orbits you need')
    args = parser.parse_args(argv)

    designations = [re.escape(d.encode('ascii')) for d in args.designations]
    pattern = rb'^((?:%s) .*\n)' % rb'|'.join(designations)
    r = re.compile(pattern, re.M)

    data = load.open(mpc.MPCORB_URL).read()
    data = zlib.decompress(data, wbits = zlib.MAX_WBITS | 16)
    lines = r.findall(data)

    sys.stdout.buffer.write(b''.join(lines))

if __name__ == '__main__':
    main(sys.argv[1:])

The same four asteroid orbits could then be extracted with:

python mpc_make_excerpt.py 00001 00002 00003 00004 > MPCORB.excerpt.DAT

Note that the minor planets file has no explicit expiration date, so load.open() in the above script will only download the file once. Subsequent calls re-open the copy of the file already on your filesystem. To force a fresh download, pass reload=True.

In either case, the resulting file — shorn of its text header, and containing only minor planet orbits — is ready for Skyfield to load.

with load.open('MPCORB.excerpt.DAT') as f:
    minor_planets = mpc.load_mpcorb_dataframe(f)

print(minor_planets.shape[0], 'minor planets loaded')

4 minor planets loaded

Some Skyfield users have encountered Minor Planet Center files with bodies whose orbital elements are incomplete, presumably because their orbits are still being determined. To avoid receiving an EphemerisRangeError exception when Skyfield tries to compute a position for these bodies, you can ask Pandas to filter them out of your dataframe:

# Filtering the orbits dataframe to avoid triggering
# an `EphemerisRangeError` on ill-defined orbits.

bad_orbits = minor_planets.semimajor_axis_au.isnull()
minor_planets = minor_planets[~bad_orbits]

As was demonstrated in the previous section on comets, you can ask Pandas to index the dataframe by minor planet designation for quick lookup.

# Index by designation for fast lookup.
minor_planets = minor_planets.set_index('designation', drop=False)

# Sample lookups.
row = minor_planets.loc['(1) Ceres']
row = minor_planets.loc['(4) Vesta']

Finally, generating a position involves the same maneuver necessary for comets: since minor planet orbits are centered on the Sun, the Sun’s position vector must be added to theirs to build a complete position.

ceres = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)
ra, dec, distance = earth.at(t).observe(ceres).radec()
print(ra)
print(dec)

05h 51m 45.85s
+22deg 38' 50.2"