Better way to get approximate ground track for a satellite using Skyfield?

Question

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.figure()
plt.plot(degs*ISSlon, degs*ISSlat, 'ok')
plt.show()

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())

ax.stock_img()

plt.scatter(subsat.longitude.degrees, subsat.latitude.degrees, transform=ccrs.PlateCarree(),
            color='red')
plt.show()