I wanted to make a quick ground-track plot, so I used the Python package Skyfield to propagate a TLE and return Earth-Centered Inertial coordinates. Then I created a Topos() object fixed to the Earth's surface at lat/lon = 0, 0 and used it to generate the rotation of the earth in order to "unwind" the ISS position on the surface.

side note: This is not a good way to do this, but it gives results good enough to make a small map for an illustration. Problems include assuming the earth's axis is still in the z direction, and of course treating the earth's surface as spherical.

Is there a better way to do this within Skyfield without using a method that starts with an underscore - in other words using methods intended to be completely public for the user? Also, is there a way to get the Earth centered, Earth fixed (i.e. rotating Earth) coordinates directly in Skyfield without unwinding like this?

EDIT: I've adjusted the script since it's been over a year and Skyfield v 1.0 has been released.

ISS_TLE = """1 25544U 98067A   16341.96974289  .00003303  00000-0  57769-4 0  9996
2 25544  51.6456 276.4739 0005937 300.1004 104.8148 15.53811586 31866"""
L1, L2 = ISS_TLE.splitlines()

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, EarthSatellite, Topos

degs     = 180./np.pi

r_earth  = 6371.  # for ROUGH approx. ground track, just use a spherical Earth

load    = Loader('~/Documents/YourNameHere/SkyData')
data    = load('de421.bsp')
earth   = data['earth']
topoZZ  = Topos(latitude_degrees=0.0, longitude_degrees=0.0)

location = earth + topoZZ

ISS      = earth + EarthSatellite(L1, L2)

ts       = load.timescale()

minutes  = np.arange(0, 140, 0.5)
time     = ts.utc(2016, 12, 7, 12, minutes)

Epos     = earth.at(time).position.km
ZZpos    = topoZZ.at(time).position.km    ## Position of (0.0N, 0.0E) to get rotation
ISSpos   = ISS.at(time).position.km - Epos

theta_ZZ = np.arctan2(ZZpos[1], ZZpos[0])   # calculate Earth's rotaion

sth, cth         = np.sin(-theta_ZZ), np.cos(-theta_ZZ) # unwind
xISS, yISS, zISS = ISSpos
xISSnew, yISSnew = xISS*cth - yISS*sth, xISS*sth + yISS*cth # rotate ISS data to match Earth
ISSnew           = np.vstack((xISSnew, yISSnew, zISS))

x, y, z = ISSnew
r       = np.sqrt((ISSpos**2).sum(axis=0))
rxy     = np.sqrt(x**2 + y**2)
ISSlat, ISSlon   = np.arctan2(z, rxy), np.arctan2(y, x)

plt.plot(degs*ISSlon, degs*ISSlat, 'ok')

enter image description here


1 Answer 1


I know this is an old question but for funsies, here's a quick script.

This question was asked when the .subpoint() feature wasn't supported in Skyfield to grab the long/lat for a satellite orbit projected onto the Earth.

Here's a quick script that shows how to use Skyfield's built in functions with plotting using cartopy to build the map and projections.

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

ts = load.timescale(builtin=True)

TLE = """ISS (ZARYA)             
1 25544U 98067A   19203.81086311  .00000606  00000-0  18099-4 0  9996
2 25544  51.6423 184.5274 0006740 168.1171 264.4057 15.50995519180787"""

name, L1, L2 = TLE.splitlines()

sat = EarthSatellite(L1, L2)

minutes = np.arange(0, 200, 0.1) # about two orbits
times   = ts.utc(2019, 7, 23, 0, minutes)

geocentric = sat.at(times)
subsat = geocentric.subpoint()

fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())


plt.scatter(subsat.longitude.degrees, subsat.latitude.degrees, transform=ccrs.PlateCarree(),
  • $\begingroup$ Excellent! For additional funsies how about a way to construct a terrestrial subpoint for another object? Comments below this question as well in this chat suggest that the OP would eventually like to plot a trail of the sub-planet points across the Earth's surface. I can think of two ways; 1) use reverse_terra() in this answer (cont.) $\endgroup$
    – uhoh
    Commented Dec 3, 2019 at 7:02
  • $\begingroup$ ...or 2) construct a Geocentric object in Skyfield based on the planet's positions and then apply .subpoint() to it. Unfortunately the text of that question hasn't been updated to reflect the OP's intent yet. (note: the linked question is in Astronomy SE not here in Space SE], it's common to confuse the two) $\endgroup$
    – uhoh
    Commented Dec 3, 2019 at 7:09
  • $\begingroup$ fyi now that you are above 50 reputation points you can leave comments on any posts you like! $\endgroup$
    – uhoh
    Commented Dec 3, 2019 at 7:11
  • $\begingroup$ I've found that the line "ax.stock_img() results in the following python error: FileNotFoundError: [Errno 2] No such file or directory: '/usr/lib/python3/dist-packages/cartopy/data/raster/natural_earth/50-natural-earth-1-downsampled.png' - it seems the folder structure has changed and that image moved. Where else is the same file available? $\endgroup$
    – Puffin
    Commented Nov 22, 2020 at 23:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.