How to Plot Satellite Ground Tracks with Python Using Skyfield ?

Introduction

Plotting satellite ground tracks is a common task in Earth observation, orbit analysis, and remote sensing validation workflows. Ground tracks are widely used to visualize satellite motion, estimate overpass times, and understand orbital geometry.

In this article, we demonstrate how to use the Python package Skyfield to plot satellite ground tracks accurately. We also clarify an important and often misunderstood concept: the difference between a ground track and an instrument swath.

The examples focus on:

• Suomi-NPP (SNPP), a sun-synchronous Earth-observing satellite
• The International Space Station (ISS), a low-inclination platform

Ground Track vs Swath: The Critical Distinction

Ground Track (Orbital Concept)

A ground track is defined as:

The path traced on Earth’s surface by the satellite’s nadir sub-point as it orbits the Earth.

Key characteristics:

• A 1-dimensional curve
• Depends only on the satellite orbit
• Independent of the onboard instrument
• Computed from orbital elements such as Two-Line Elements (TLEs)

Ground tracks are commonly used for:

• Orbit visualization
• Repeat-cycle analysis
• Overpass timing estimation
• Educational and mission-planning tools (e.g., NOAA iSTRaK)

Swath (Instrument Concept)

A swath is defined as:

The 2-dimensional area on Earth observed by an instrument during an overpass.

For VIIRS onboard Suomi-NPP:

• Cross-track scanning radiometer
• Approximate swath width: ~3040 km
• Defined by instrument field-of-view and scan geometry

Swaths are obtained from:

• VIIRS geolocation products: VNP03MOD, VNP03IMG
• VIIRS fire products: VNP14MOD, VNP14IMG (a subset of the swath)

A swath contains the ground track, but the two are not the same thing.

This distinction is critical: many plots labeled “VIIRS ground track” are actually showing swath coverage or pixel locations, which is technically incorrect.

Installing Skyfield in Python

Skyfield is a Python library for high-precision satellite orbit propagation based on JPL ephemerides. It is well suited for computing satellite ground tracks, sub-satellite points, and observation geometry for Earth-observing missions such as Suomi-NPP (VIIRS).

Option 1: Install with Conda (recommended)

Using conda is the easiest and most reliable way to install Skyfield because all dependencies are handled automatically.

1

conda install -c conda-forge skyfield

This installs:

• skyfield
• numpy
• jplephem
• sgp4

All of which are required for satellite orbit propagation.

If you are working in a dedicated environment:

1
2
3

conda create -n skyfield_env python=3.11
conda activate skyfield_env
conda install -c conda-forge skyfield

Option 2: Install with pip

You can also install Skyfield using pip, either system-wide or inside a virtual environment:

1

pip install skyfield

To avoid dependency conflicts, it is strongly recommended to use a virtual environment:

1
2
3

python -m venv skyfield_env
source skyfield_env/bin/activate   # Linux / macOS
pip install skyfield

Optional: Verify the Installation

You can verify that Skyfield is correctly installed by running:

1
2

from skyfield.api import load, EarthSatellite
print("Skyfield installed successfully")

If no error is raised, the installation is complete.

Notes on TLE Data

Skyfield propagates satellite orbits using Two-Line Element (TLE) sets. For accurate ground tracks, always use recent TLEs, for example from:

• CelesTrak (NORAD catalog)
• Space-Track.org

Outdated TLEs can lead to noticeable positional errors, especially for low-Earth-orbit satellites like Suomi-NPP.

Full Working Example: Suomi-NPP Ground Track

Below is a complete, NOAA-iSTRaK-like example showing a few full Suomi-NPP orbits, plotted correctly using geodesic curves.

Python Code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

from skyfield.api import load
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# Load up-to-date SNPP TLE
ts = load.timescale()
tle_url = "https://celestrak.org/NORAD/elements/gp.php?CATNR=37849&FORMAT=TLE"
sat = load.tle_file(tle_url)[0]

print("TLE epoch:", sat.epoch.utc_strftime())

# Parameters
START_TIME = (2025, 12, 18, 0, 0, 0)  # UTC
ORBIT_PERIOD_SEC = 99 * 60           # ~99 minutes
DT = 20                              # sampling (seconds)
N_ORBITS = 4                         # number of full orbits to plot

# Time array
t = ts.utc(
    START_TIME[0], START_TIME[1], START_TIME[2],
    START_TIME[3], START_TIME[4],
    np.arange(0, ORBIT_PERIOD_SEC * N_ORBITS, DT)
)

# Satellite subpoint (ground track)
sp = sat.at(t).subpoint()
lats = sp.latitude.degrees
lons = (sp.longitude.degrees + 180) % 360 - 180

# Split into individual orbits
points_per_orbit = int(ORBIT_PERIOD_SEC / DT)
tracks = [
    (lons[i*points_per_orbit:(i+1)*points_per_orbit],
     lats[i*points_per_orbit:(i+1)*points_per_orbit])
    for i in range(N_ORBITS)
]

# Plot
fig = plt.figure(figsize=(13, 6))
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
ax.gridlines(linewidth=0.5, linestyle="--")

for i, (lon, lat) in enumerate(tracks):
    ax.plot(
        lon, lat,
        transform=ccrs.Geodetic(),
        linewidth=1.8,
        label=f"Orbit {i+1}"
    )

ax.legend()
ax.set_title("SNPP Orbital Ground Tracks (Nadir Subpoint)")
plt.savefig("SNPP_Orbital_Ground_Tracks.png", bbox_inches="tight", dpi=100)
plt.show()

Code Explanation (Step by Step)

1. Why Skyfield?

Skyfield performs high-precision orbital propagation using the SGP4 model and TLEs, making it ideal for satellite ground-track calculations.

2. Why the Date Matters (Very Important)

1

START_TIME = (2025, 12, 18, 0, 0, 0)

TLEs are time-dependent. They are typically valid only within ~1–2 weeks of their epoch.

If you:

• Use an old TLE
• Propagate it far into the future or past

➡️ The computed orbit will be physically incorrect, even though the code runs.

This is why we:

• Load current TLEs from CelesTrak
• Match the propagation date to the TLE epoch

You can verify TLE validity with:

1

print(sat.epoch.utc_strftime())

3. Choosing Orbit Parameters

1
2
3

ORBIT_PERIOD_SEC = 99 * 60  # ~99 minutes
DT = 20                    # sampling (seconds)
N_ORBITS = 4               # number of orbits

• ORBIT_PERIOD_SEC: Approximate orbital period of SNPP
• DT: Temporal resolution of the ground track
• N_ORBITS: Controls plot readability

➡️ Plotting a few complete orbits is far clearer than plotting continuous tracks over many hours.

4. Why ccrs.Geodetic() Matters

1

ax.plot(lon, lat, transform=ccrs.Geodetic())

This ensures:

• Curved paths are plotted as great-circle arcs
• The track resembles official NOAA visualizations
• No artificial straight-line artifacts appear

Comparison with NOAA iSTRaK

The resulting plot can be directly compared with the official NOAA tool:

https://viirs.umd.edu/WebGL_traj/orb_npp/disp_npp_20251218.php

Both show:

• Individual orbit passes
• Westward drift due to Earth rotation
• Smooth, curved ground tracks

ISS Ground Track Example

The same approach can be applied to other satellites, such as the International Space Station (ISS).

Key Differences

• Low inclination (~51.6°)
• Shorter orbital period (~92 minutes)
• No polar coverage

ISS Example Code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

# Load ISS TLE
tle_url_iss = "https://celestrak.org/NORAD/elements/stations.txt"
sats = load.tle_file(tle_url_iss)
iss = [s for s in sats if s.name == "ISS (ZARYA)"][0]

# Parameters
START_TIME = (2025, 12, 18, 0, 0, 0)
ORBIT_PERIOD_SEC = 92 * 60
DT = 20
N_ORBITS = 3

t = ts.utc(
    START_TIME[0], START_TIME[1], START_TIME[2],
    START_TIME[3], START_TIME[4],
    np.arange(0, ORBIT_PERIOD_SEC * N_ORBITS, DT)
)

sp = iss.at(t).subpoint()
lats = sp.latitude.degrees
lons = (sp.longitude.degrees + 180) % 360 - 180

points_per_orbit = int(ORBIT_PERIOD_SEC / DT)
tracks = [
    (lons[i*points_per_orbit:(i+1)*points_per_orbit],
     lats[i*points_per_orbit:(i+1)*points_per_orbit])
    for i in range(N_ORBITS)
]

fig = plt.figure(figsize=(13, 6))
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
ax.gridlines()

for i, (lon, lat) in enumerate(tracks):
    ax.plot(lon, lat, transform=ccrs.Geodetic(), label=f"ISS Orbit {i+1}")

ax.legend()
ax.set_title("ISS Orbital Ground Tracks (Nadir Subpoint)")
plt.show()

References