SGP4 propagation in Python with skyfield

A TLE describes where the orbit is; an SGP4 propagator computes where the satellite is at any given time. This week we connect the two using Python — the same code path that powers every satellite tracking app on the internet.

What SGP4 does

Simplified General Perturbations 4 is the orbital propagation model standardized by NORAD for use with TLEs. It accounts for the first-order perturbations of low-Earth orbit: Earth's J2 oblateness, atmospheric drag, solar/lunar gravity, and solar radiation pressure. It's accurate to within ~1 km for a fresh TLE and degrades over time as the underlying model assumptions diverge from reality.

SGP4 is not the most accurate propagator available — for centimeter-precision (GPS satellites, geodetic missions) you'd use a numerical propagator with high-fidelity models — but it's accurate enough for almost every satellite-tracking and ground-station application, and it's standardized across the industry.

skyfield: Python's reference implementation

skyfield is a high-precision astronomy + satellite-propagation library that wraps the canonical SGP4 implementation and adds proper time-system handling (UT1 vs UTC vs TT vs TDB), Earth orientation, and a clean Python API. It is the right default choice; the older sgp4 standalone package is lower-level and lacks the time-system polish.

from skyfield.api import EarthSatellite, load, wgs84

tle_lines = [
    "1 25544U 98067A   24130.50145833  .00018539  00000-0  33188-3 0  9994",
    "2 25544  51.6406 348.5395 0006703 117.9568 358.1729 15.50289267449420",
]
ts = load.timescale()
iss = EarthSatellite(tle_lines[0], tle_lines[1], "ISS", ts)
t = ts.now()
geocentric = iss.at(t)
subpoint = wgs84.subpoint(geocentric)
print(f"ISS sub-satellite point: ({subpoint.latitude.degrees:.4f}, {subpoint.longitude.degrees:.4f}), alt {subpoint.elevation.km:.1f} km")

Time matters

Satellite propagation is exquisitely time-sensitive. A 1-second error means a 7 km position error for an ISS-like LEO orbit. skyfield handles four time systems: UTC, UT1, TT (Terrestrial Time), and TDB (Barycentric Dynamical Time). TLE epochs are in UT1 (Universal Time 1), which is Earth-rotation-based and differs from UTC by up to ±0.9 seconds — skyfield handles the conversion automatically via the IERS (International Earth Rotation Service) tables it downloads on first use.

Always use skyfield's ts.utc(...), ts.tt(...), etc., to construct times. Never compute time offsets manually.

Sub-satellite points and ground tracks

The sub-satellite point is the lat/lon directly below the satellite (the point on Earth's surface intersected by the line from Earth's center through the satellite). For a 24-hour ground track, you propagate to a series of times, extract the sub-satellite point at each, and stitch them into a polyline.

from skyfield.api import load
import numpy as np

t0 = ts.now()
times = [ts.tt_jd(t0.tt + i/1440) for i in range(0, 1440, 1)]  # 1-minute steps for 24h
positions = [(wgs84.subpoint(iss.at(t)).longitude.degrees,
              wgs84.subpoint(iss.at(t)).latitude.degrees) for t in times]

One subtlety: ground tracks cross the antimeridian (180° / -180°) and must be split there into separate line segments for proper visualization, or the line will draw across the entire map. The lab notebook handles this with a simple jump-detection split.

The lab

You'll generate the ISS's next 24-hour ground track at 1-minute resolution (1,440 points) as a GeoJSON LineString with embedded timestamps. Then you'll visualize it in QGIS over a Natural Earth basemap, color-coded by altitude (use a small vertical exaggeration since ISS altitude varies only ~10 km over an orbit).

This is the same pipeline LaunchDetect uses for the live satellite tracker in production — the only differences are streaming updates (Week 19) and a 3D globe instead of a 2D map (Week 18).