Skyfield: HomeTable of ContentsAPI Reference

Example Plots

This section of the documentation will gradually accumulate example scripts for producing images from Skyfield computations.

Note that these example scripts are written for fairly recent versions of matplotlib. If you try running them on a system with an older version of the library, you might see errors — in particular with how they specify colors, in which case you can try omitting those parameters to get the script running. In any case, these are only intended to be a starting point for building your own scripts, either with matplotlib or whatever other plotting library you prefer.

Drawing a finder chart for comet NEOWISE

Here is a stand-alone script that brings together four different data sources — a planetary ephemeris, a comet orbit database, a large star catalog, and constellation diagrams — to plot the course of Comet NEOWISE across Ursa Major over one week of July 2020:

_images/neowise-finder-chart.png

Its code includes many design decisions and presentation tweaks that you will probably want to adjust for your own project. Use the script as a starting point:

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection

from skyfield import api
from skyfield.api import Star, load, tau
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.data import hipparcos, mpc, stellarium
from skyfield.projections import build_stereographic_projection

# The comet is plotted on several dates `t_comet`.  But the stars only
# need to be drawn once, so we take the middle comet date as the single
# time `t` we use for everything else.

ts = load.timescale()
t_comet = ts.utc(2020, 7, range(17, 27))
t = t_comet[len(t_comet) // 2]  # middle date

# An ephemeris from the JPL provides Sun and Earth positions.

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

# The Minor Planet Center data file provides the comet orbit.

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

comets = comets.set_index('designation', drop=False)
row = comets.loc['C/2020 F3 (NEOWISE)']
comet = sun + mpc.comet_orbit(row, ts, GM_SUN)

# The Hipparcos mission provides our star catalog.

with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)

# And the constellation outlines come from Stellarium.  We make a list
# of the stars at which each edge stars, and the star at which each edge
# ends.

url = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
       '/skycultures/western_SnT/constellationship.fab')

with load.open(url) as f:
    constellations = stellarium.parse_constellations(f)

edges = [edge for name, edges in constellations for edge in edges]
edges_star1 = [star1 for star1, star2 in edges]
edges_star2 = [star2 for star1, star2 in edges]

# We will center the chart on the comet's middle position.

center = earth.at(t).observe(comet)
projection = build_stereographic_projection(center)
field_of_view_degrees = 45.0
limiting_magnitude = 7.0

# Now that we have constructed our projection, compute the x and y
# coordinates that each star and the comet will have on the plot.

star_positions = earth.at(t).observe(Star.from_dataframe(stars))
stars['x'], stars['y'] = projection(star_positions)

comet_x, comet_y = projection(earth.at(t_comet).observe(comet))

# Create a True/False mask marking the stars bright enough to be
# included in our plot.  And go ahead and compute how large their
# markers will be on the plot.

bright_stars = (stars.magnitude <= limiting_magnitude)
magnitude = stars['magnitude'][bright_stars]
marker_size = (0.5 + limiting_magnitude - magnitude) ** 2.0

# The constellation lines will each begin at the x,y of one star and end
# at the x,y of another.  We have to "rollaxis" the resulting coordinate
# array into the shape that matplotlib expects.

xy1 = stars[['x', 'y']].loc[edges_star1].values
xy2 = stars[['x', 'y']].loc[edges_star2].values
lines_xy = np.rollaxis(np.array([xy1, xy2]), 1)

# Time to build the figure!

fig, ax = plt.subplots(figsize=[9, 9])

# Draw the constellation lines.

ax.add_collection(LineCollection(lines_xy, colors='#00f2'))

# Draw the stars.

ax.scatter(stars['x'][bright_stars], stars['y'][bright_stars],
           s=marker_size, color='k')

# Draw the comet positions, and label them with dates.

comet_color = '#f00'
offset = 0.002

ax.plot(comet_x, comet_y, '+', c=comet_color, zorder=3)

for xi, yi, tstr in zip(comet_x, comet_y, t_comet.utc_strftime('%-m/%d')):
    text = ax.text(xi + offset, yi - offset, tstr, color=comet_color,
                   ha='left', va='top', fontsize=9, weight='bold', zorder=-1)
    text.set_alpha(0.5)

# Finally, title the plot and set some final parameters.

angle = np.pi - field_of_view_degrees / 360.0 * np.pi
limit = np.sin(angle) / (1.0 - np.cos(angle))

ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.set_aspect(1.0)
ax.set_title('Comet NEOWISE {} through {}'.format(
    t_comet[0].utc_strftime('%Y %B %d'),
    t_comet[-1].utc_strftime('%Y %B %d'),
))

# Save.

fig.savefig('neowise-finder-chart.png', bbox_inches='tight')

If you choose a different rendering engine instead of the venerable but rather ornery and complicated matplotlib, then of course the plotting calls you make will be completely different. But the basic data loading and filtering will be the same, so hopefully the script will still help get you started in targeting a more modern plotting library.

Plotting satellite altitude during re-entry

Here is the decreasing altitude of a satellite as its orbit decayed and it re-entered the atmosphere above the Pacific Ocean:

_images/goce-reentry.png

The code to produce the diagram using matplotlib, including custom tick marks that are based on the date, is:

from matplotlib import pyplot as plt
from matplotlib.dates import HourLocator, DateFormatter

from numpy import arange

def label_dates_and_hours(axes):
    axes.xaxis.set_major_locator(HourLocator([0]))
    axes.xaxis.set_minor_locator(HourLocator([0, 12]))
    axes.xaxis.set_major_formatter(DateFormatter('\n%a %d'))
    axes.xaxis.set_minor_formatter(DateFormatter('%Hh'))

from skyfield.api import load, EarthSatellite

# Load the satellite's final TLE entry.

sat = EarthSatellite(
    '1 34602U 09013A   13314.96046236  .14220718  20669-5  50412-4 0   930',
    '2 34602 096.5717 344.5256 0009826 296.2811 064.0942 16.58673376272979',
    'GOCE',
)

# Build the time range `t` over which to plot, plus other values.

ts = load.timescale()
t = ts.tt_jd(arange(sat.epoch.tt - 1.0, sat.epoch.tt + 3.0, 0.01))
reentry = ts.utc(2013, 11, 11, 0, 16)
earth_radius_km = 6371.

# Start a new figure.

fig, ax = plt.subplots()

# Draw the blue curve.

x = t.toordinal()
y = sat.at(t).distance().km - earth_radius_km
ax.plot(x, y)

# Label the official moment of reentry.

x = reentry.toordinal()
y = sat.at(reentry).distance().km - earth_radius_km
ax.plot(x, y, 'ro')
ax.text(x, y + 10, 'Moment of re-entry')

# Grid lines and labels.

label_dates_and_hours(ax)
ax.grid()
ax.set(title='GOCE satellite altitude', ylabel='km above sea level')

# Render the plot to a PNG file.

fig.savefig('goce-reentry.png', bbox_inches='tight')