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, ZZpos) # 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()