Efficient way to propagate satellite catalog over time

Question

Problem statement

I need to propagate the entire catalog of recent TLEs (need a free account to view) from space-track.org using skyfield or similar. There are typically 15k-16k TLEs in the list. I have it working, but it is very slow. On the order of hours using 46 cores on a server.

Not married to skyfield. If astropy or pyephem or something else is faster, I would happily accept an answer that shows what I am trying to do using that.

Minimal example

For my application, I am loading TLEs into a Pandas dataframe and doing my analysis there, so I will keep my example in the Pandas world. Minimal example follows.

Assuming the satellite catalog was saved as catalog.txt, set up the environment, then read TLEs, generate sf.sgp4lib.EarthSatellite objects, and load everything into a Pandas dataframe. We also offset positions to some observation site. I leave picking an observation site up to the reader (0, 0, 0 would be fine):

import skyfield as sf
import pandas as pd
from skyfield.api import load, Topos
from datetime import datetime, timezone, timedelta

with open('catalog.txt', 'r') as f:
    tle_list = [line.strip() for line in f.read().split('\n')
                if line is not '']
data = []
for i in range(0, len(tle), 2):  # every two lines
    temp = {}
    temp['tle1'] = tle_list[i]
    temp['tle2'] = tle_list[i+1]
    temp['earthsat'] = sf.sgp4lib.EarthSatellite(tle_list[i],
                                                 tle_list[i+1])
    data.append(temp)
df = pd.DataFrame(data=data)
site = Topos(latitude_degrees=site_lat,
             longitude_degrees=site_lon,
             elevation_m=site_alt)
df['earthsat'] = df.earthsat - site  # offset to site location

2.1 s ± 20 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Create an array of timezone-aware datetime objects over which to propagate all of the satellites. Here I chose every 10 minutes from 4 hours before midnight on the day that I write this to 4 hours after.

ts = load.timescale()
tz = timezone(timedelta(hours=-4))  # Eastern or whatever timezone
midnight = datetime(2018, 4, 4, 0, 0, 0, tzinfo=tz)  # midnight today
start = midnight - timedelta(hours=4)
end = midnight + timedelta(hours=4)
delta_time = timedelta(minutes=10)
# This is ugly, but I had issues using linspace or arange...
times = [start]
now = start
while now <= end:
    now += delta_time
    times.append(now)

189 ms ± 36.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Finally, calculate the astrometric positions for each time step for each satellite. This is what takes forever. Too long for me to run again just to get timing, but it is on the order of hours using 46 cores on a server.

df['astrometric'] = df.earthsat.apply(lambda x: [x.at(ts.utc(time)) for time in times])

Additional details

I found this discussion on date arrays in the documentation that suggests passing the whole array in at once: x.at(ts.utc(times)). So far that takes fewer cores and runs a bit faster, but still takes prohibitively long.

I got around this slowdown for a bit by creating generators for astrometrics (originally why I went away from passing the whole times array in at once), but eventually I actually need to evaluate things, so I am unable to avoid the heavy lifting forever.

In case the final use case lends itself to some specific speed-up, I ultimately need to get observation angles from site out of these objects, so [x.altaz() for x in row.astrometrics] type of thing.

Thoughts on a solution

My thought right now is that I am calculating the observation site's position at every time step throughout the night for every satellite in the catalog. I could be wrong, but if I am right then I imagine I would see a pretty decent speedup (maybe not enough) by calculating this once and then looking it up for each new satellite. Does anybody know how to do this?

Also, if there is a faster implementation of the orbit propagator out there, or a way to speed up skyfield's implementation, I would happily accept an answer that shows how to do what I am trying to do using that (thus including astropy and pyephem tags).

Thank you.

Answer 1

My best advice is to use the NASA SPICE Toolkit for this type of work. This let's you load in Two-Line Elements, and then use NASA's NAIF/SPICE kernels to do the rest of the work (you might also find some satellites in other formats, but TLE are fine).

If you were working in C, you would use the getelm_c method to read them and spkpos_c to get the position. Luckily, there is a python wrapper called spiceypy!

The getelm_c method is wrapped in spiceypy.spiceypy.getelm(frstyr, lineln, lines) to load in your TLE. You'd them want to use spiceypy.spiceypy.spkpos(targ, et, ref, abcorr, obs) to get your position relative to the reference body.

I suggest using the docs from spiceypy, specifically the position of Cassini example, to make sure you have all the correct files in place if you decide to go with NASA's SPICE kernels: https://spiceypy.readthedocs.io/en/master/exampleone.html