Memory issue with satellite.at()?

[lucabaldini]

I have an application where I am trying to calculate the position of a satellite on a very fine time grid (say millions of points).

The simple call to
geocentric = satellite.at(t)
seems to run into performance issue, where for large n scales much worse than linearly, and by the time n is O(1 M) the thing is essentially requiring >10 GB or RAM and bringing my terminal down.

I can make a more detailed report, including all the related information, but I want to make a quick check first of whether I am doing something wrong and/or you have any immediate insight into this.

Thanks in advance!

[brandon-rhodes]

Good question! On the one hand, with a quick script I cannot confirm performance worse than linear — in fact I show slightly better than linear, with the cost per point generated dropping as the size of the time vector grows larger:

from sgp4.api import accelerated
from skyfield.api import EarthSatellite, load
from time import time

assert accelerated

ts = load.timescale()

times = []
for n in 10, 100, 1000, 10000, 100000:
    t = ts.utc(2020, 5, 12, 19, 16, range(n))
    times.append(t)

for t in times:
    t0 = time()
    sat = EarthSatellite(
        '1 25544U 98067A   20128.53256469  .00000758  00000-0  21651-4 0  9992',
        '2 25544  51.6442 189.6875 0001104 261.0974 196.9979 15.49353594225718',
    )
    p = sat.at(t)
    duration = time() - t0
    size, = t.shape
    print(duration, size, duration / size)

But, on the other hand, my laptop runs out of memory trying to generate a mere 100k positions!

0.008205175399780273 10 0.0008205175399780273
0.035399675369262695 100 0.00035399675369262695
0.3386518955230713 1000 0.0003386518955230713
2.879323959350586 10000 0.0002879323959350586
zsh: killed     python tmp55.py

It should only take about 4.8 MB to hold 100,000 points × 3 coordinates (x,y,z) × 8 bytes per floating point number × 2 for both position and velocity, so there’s definitely something unexpected happening with memory allocation.

Let me go look at the sgp4 module’s code for what might be going wrong, and I’ll report back!

[lucabaldini]

Thanks---I appreciate.

(Note this is not a real blocker for me, although I could definitely used some speed-up.)

I know anecdotal evidence isn't really useful for a bug report, but I am also under the impression that when doing the same thing twice the second is much faster. Does it make sense?

[brandon-rhodes]

Yes, it makes perfect sense given the results of my investigation! Here’s what I found:

The cost here is not in computing the satellite positions with sgp4 which is running fine. The problem is in computing the high-precision orientation of the Earth that allows those satellite coordinates to be rotated into the Solar System’s frame of reference that Skyfield uses internally. The expensive but highly accurate IAU2000A nutation model needs more work — the way I have implemented it in Python must create an intermediate NumPy array so gargantuan that it runs our machines out of memory. In my case I can't even compute when t is length 100,000, as we saw above.

The reason you are seeing speedup is probably that you are re-using the same Time object. The Time object caches the results of the nutation calculation, so it only has to be run the first time you use that Time object to translate coordinates to or from the Earth's reference frame.

So we should think about three things: immediately getting you a workaround; second, thinking about whether the workaround should kick in automatically for Earth satellite coordinates; and, third, looking into the IAU2000A implementation to see if those huge intermediate arrays can be avoided.

So, first:

Here's a script that measures the degradation in coordinates if you switch to the less accurate but much less time-and-memory-intensive IAU2000B:

from skyfield.api import EarthSatellite, load
from skyfield.functions import length_of
from skyfield.nutationlib import iau2000b

sat = EarthSatellite(
    '1 25544U 98067A   20128.53256469  .00000758  00000-0  21651-4 0  9992',
    '2 25544  51.6442 189.6875 0001104 261.0974 196.9979 15.49353594225718',
)
h = range(10000)

ts = load.timescale()
t = ts.utc(2020, 5, 12, h)
p1 = sat.at(t)

t = ts.utc(2020, 5, 12, h)
t._nutation_angles = iau2000b(t.tt)
p2 = sat.at(t)

delta = length_of(p1.position.m - p2.position.m)

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(h, delta)
ax.set(xlabel='Hours', ylabel='Meters difference in position',
       title='Satellite position with iau2000a vs iau2000b')
ax.grid()
fig.savefig('tmp.png')

And the result:

Given that satellite coordinates are generally only accurate to a kilometer or so, this tiny 1mm or 2mm difference in position is negligible. I encourage you to immediately update your script create time objects where you set t._nutation_angles = iau2000b(t.tt) like in the script above, and see if your computation time (and memory footprint) improves.

Second:

Delivering a quicker and less accurate ICRF matrix in cases where Earth satellite objects are involved will be tricky. The actual coordinate transform routine just asks the time for its M matrix which in turn is build from N, P, and B (all this is in timelib.py), and the computation for N asks the time for its _earth_tilt attribute, which asks the time for its _nutation_angles attribute which is what costs so much to compute.

So would we invent a new rotation matrix M_low_accuracy that asks for N_low_accuracy which asks for _earth_tilt_low_accuracy which asks for _nutation_angles_low_accuracy which calls the iau2000b() routine? It might almost be worth it given the speedup it would provide for the many users that call Skyfield for Earth satellite positions. But I will have to think through whether the extra complexity could be trimmed back somehow!

Third:

I am not sure why those intermediate matrices get so large. I'll have to debug the iau2000a() function sometime with lots of print statements.

[lucabaldini]

Wow! The
t._nutation_angles = iau2000b(t.tt)
really saved my day.

I haven't made extensive tests, yet, but at the first glance the time/memory footprint is well below the level at which I would notice it.

I don't know if mine is really a peculiar use case, but it might be worth noting this in the docs.

Thanks!

[brandon-rhodes]

@lucabaldini — I spent yesterday looking at the nutation routines, and I found several spots where very large intermediate results were being produced, and figured out how to use numpy.dot() to go ahead and sum the results as they’re produced so the full intermediate result never gets stored.

Here's how to try out the new version if you would like to let me know ahead of the next release if you see a difference in your own calculations:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

[brandon-rhodes]

The new code has now been released in Skyfield 1.21. I'm going to close the issue for now, but please re-open if you run into further problems!