# How can I plot a satellite's orbit in 3D from a TLE using Python and Skyfield?

I have obtained a Two Line Element (TLE) of a satellite in Earth orbit from Celestrak at https://www.celestrak.com/NORAD/elements/ and I would like to use it to calculate an orbit.

How can I calculate the orbit from the TLE, and then plot it in 3D, using Python and the Skyfield package?

### Code updated for python3*

def makecubelimits(axis, centers=None, hw=None):
lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
if centers == None:
centers = [0.5*sum(pair) for pair in lims]

if hw == None:
widths  = [pair[1] - pair[0] for pair in lims]
hw      = 0.5*max(widths)
ax.set_xlim(centers[0]-hw, centers[0]+hw)
ax.set_ylim(centers[1]-hw, centers[1]+hw)
ax.set_zlim(centers[2]-hw, centers[2]+hw)
print("hw was None so set to:", hw)
else:
try:
hwx, hwy, hwz = hw
print("ok hw requested: ", hwx, hwy, hwz)

ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
except:
print("nope hw requested: ", hw)
ax.set_xlim(centers[0]-hw, centers[0]+hw)
ax.set_ylim(centers[1]-hw, centers[1]+hw)
ax.set_zlim(centers[2]-hw, centers[2]+hw)

return centers, hw

TLE = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()

from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]

earth   = planets['earth']

hours = np.arange(0, 3, 0.01)

time = ts.utc(2018, 2, 7, hours)

print(Rpos.shape)

re = 6378.

theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in [np.cos, np.sin, np.zeros_like]]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
cph, sph = [f(phi) for f in [np.cos, np.sin]]
lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
lon0[1]*cph + lon0[0]*sph,
lon0[2]) )
lons.append(lon)

lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
cph, sph = [f(phi) for f in [np.cos, np.sin]]
lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
lats.append(lat)

if True:
fig = plt.figure(figsize=[10, 8])  # [12, 10]

ax  = fig.add_subplot(1, 1, 1, projection='3d')

x, y, z = Rpos
ax.plot(x, y, z)
for x, y, z in lons:
ax.plot(x, y, z, '-k')
for x, y, z in lats:
ax.plot(x, y, z, '-k')

centers, hw = makecubelimits(ax)

print("centers are: ", centers)
print("hw is:       ", hw)

plt.show()

if True:
plt.figure()
plt.xlabel('hours', fontsize=14)
plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
plt.show()

• I've accepted this answer since we should be using Python 3 now. Some relevant details are still helpful in this answer with the old Python 2. – uhoh Dec 27 '19 at 10:32

SGP4 is the standard procedure that TLEs are intended to work with. All of the following are extremely helpful, but the most important point would be use a standard, recent SGP4 package that is recommended, do not try to use the elements in a TLE yourself. This is becuase the TLE and the SGP4 package are built for each other.

One interesting point is that for high altitude orbits with periods longer than 225 minutes, a modern SGP4 implementation will switch to an algorithm that is called SDP4. See the question “Deep space” corrections in SGP4; how does it account for the Sun's and Moon's gravity? for more on that.

The easiest to access SGP4 that I know is in the Python package Skyfield. You can find SGP4 routines in many languages in many places. I would recommend you choose something that is used widely and well maintained, and not just any random code on the internet.

Skyfield's SGP4 is based on :https://github.com/brandon-rhodes/python-sgp4

Here's a plot of the low earth orbit of the SpaceX F9 upper stage known as Roadster generated using Skyfield and Roadster's TLE, of course before the second burn that put it in heliocentric orbit.

def makecubelimits(axis, centers=None, hw=None):
lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
if centers == None:
centers = [0.5*sum(pair) for pair in lims]

if hw == None:
widths  = [pair[1] - pair[0] for pair in lims]
hw      = 0.5*max(widths)
ax.set_xlim(centers[0]-hw, centers[0]+hw)
ax.set_ylim(centers[1]-hw, centers[1]+hw)
ax.set_zlim(centers[2]-hw, centers[2]+hw)
print("hw was None so set to:", hw)
else:
try:
hwx, hwy, hwz = hw
print("ok hw requested: ", hwx, hwy, hwz)

ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
except:
print("nope hw requested: ", hw)
ax.set_xlim(centers[0]-hw, centers[0]+hw)
ax.set_ylim(centers[1]-hw, centers[1]+hw)
ax.set_zlim(centers[2]-hw, centers[2]+hw)

return centers, hw

TLE = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()

from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

earth   = planets['earth']

hours = np.arange(0, 3, 0.01)

time = ts.utc(2018, 2, 7, hours)

print(Rpos.shape)

re = 6378.

theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
cph, sph = [f(phi) for f in (np.cos, np.sin)]
lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
lon0[1]*cph + lon0[0]*sph,
lon0[2]) )
lons.append(lon)

lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
cph, sph = [f(phi) for f in (np.cos, np.sin)]
lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
lats.append(lat)

if True:
fig = plt.figure(figsize=[10, 8])  # [12, 10]

ax  = fig.add_subplot(1, 1, 1, projection='3d')

x, y, z = Rpos
ax.plot(x, y, z)
for x, y, z in lons:
ax.plot(x, y, z, '-k')
for x, y, z in lats:
ax.plot(x, y, z, '-k')

centers, hw = makecubelimits(ax)

print("centers are: ", centers)
print("hw is:       ", hw)

plt.show()

if True:
plt.figure()

• @EricDuminil that's great! It looks like you can delete one of the two occurrences of load() in your pastebin script, and that you've omitted the definition of load load = Loader(path) Also I'm not sure how to use the comments at the top, but they look intriguing. Usually all my Skyfield scripts need adding the extra parenthesis around tuples like here[f*np.pi for f in 0.5, 1, 2] Did anything else actually require changing for 2 to 3? – uhoh May 14 '18 at 14:11
• @uhoh: My goal was to get a working Python3 script with the least amount of changes to your code. The comments at first are just indications on how to install the required packages in an Anaconda environment (anaconda.com/download). There's also a link to download the JPL ephemeris. load is defined directly in from skyfield.api import load. The only required changes for Python3 were parens around print arguments and tuples. – Eric Duminil May 14 '18 at 15:17