Vectorized find_events computation

[jvarada]

Is there a way to call the find_events method in some sort of vectorized fashion to compute the visibility windows of a satellite at multiple positions (I'm imagining something like passing a list of Topos objects for a given start/end time and altitude)?

Edit: Also interested in computing find_events at the same location/times/altitude for multiple satellites.

[brandon-rhodes]

There are not built-in routines for either kind of vectorization, I don't think, and I'm not sure how large any speed increase would be — the underlying satellite position computation doesn't use vector math, but (hopefully on modern machines at least) is a C++ routine that is called n times if n positions are needed.

Your second case — same location, many satellites — is interesting, though, because the underlying C++ does at least have the ability to produce NumPy arrays from a loop over n satellite objects rather than having to create dozens or hundreds of separate Python float objects. Let me go see what I can put together!

[brandon-rhodes]

So: again, your first situation is not likely to support a vectorized speedup, because the number of computations necessary to do the computation is large compared to the number of satellites, so even if there were a way to somehow run it as vectors (and there's probably not? feel free to suggest a possible vectorization if you can think of one and we can consider it!), you'd only eliminate one iteration per satellite which I think is many fewer than the number of positions each satellite needs to be computed for to get rise and set times.

But for asking at time t which satellites are above the horizon can benefit directly from the underlying C++ library's ability to consider a vector of satellites. And better yet, if we want to be low-level we can skip converting all the vectors into Skyfield's base ICRF reference frame and can instead leave all the vectors in an Earth-fixed reference frame. See if, given a date and TLE file, this script produces the altitude in degrees (ignoring atmospheric refraction) of all the satellites — and, if you have a few minutes, try comparing its runtime to what you've been doing already:

# Generate Earth location.

import numpy as np
from skyfield import api
from skyfield.functions import length_of, dots
from skyfield.sgp4lib import TEME_to_ITRF
from sgp4.api import SatrecArray

ts = api.load.timescale(builtin=True)
t = ts.utc(2020, 4, 28, 11, 44)
topos = api.Topos('40.8939 N', '83.8917 W')
p = topos.itrf_xyz().km

# Build satellite array.

sats = api.load.tle_file('stations.txt')
print(sats[0])

sat_array = SatrecArray([sat.model for sat in sats])

jd = np.array([t._utc_float()])
fr = jd * 0.0
e, r, v = sat_array.sgp4(jd, fr)

# r shape = v shape = n satellites × 1 time × 3 coordinates

r = r[:,0]
v = v[:,0]

# r shape = v shape = n satellites × 3 coordinates

r = r.swapaxes(0, 1)
v = v.swapaxes(0, 1)

# r shape = v shape = 3 coordinates × n satellites

# Convert to ITRF.

r, v = TEME_to_ITRF(t.ut1, r, v)

# Offset the satellite vectors so they are centered on our Earth
# position rather than on the Earth's center.

r = r - p[:,None]

# Since the dot product of two vectors reflects not only the angle
# between them but also their magnitudes, we first turn our vectors into
# unit vectors.

r_unit = r / length_of(r)
p_unit = p / length_of(p)

zenith_angle_cosine = dots(r_unit, p_unit[:,None])
zenith_angle = np.arccos(zenith_angle_cosine)
altitude = api.tau / 4.0 - zenith_angle
altitude_degrees = api.Angle(radians=altitude).degrees

# Sort satellites by their altitude to check our result.

for i in np.argsort(altitude_degrees):
    sat = sats[i]
    alt = altitude_degrees[i]
    print('%6.2f  %5d  %s' % (alt, sat.model.satnum, sat.name))

[brandon-rhodes]

@jvarada — I wanted to check back and see whether the example code worked for your use case or not; if it did, we could think about whether similar code might one day become a feature of Skyfield.

[jvarada]

@brandon-rhodes Apologies for the super late followup to this... Another priority came up...

Last time I had tried the code above I had some success as it was written (although I could never quite get it to work over multiple satellites and multiple times). However the code fails to run with the current versions of the APIs (skyfield 1.39, sgp4 2.20):

import numpy as np
from skyfield import api
from skyfield.functions import length_of, dots
from skyfield.sgp4lib import TEME_to_ITRF
from sgp4.api import SatrecArray, jday
from datetime import datetime, timedelta, timezone

ts = api.load.timescale(builtin=True)
t = ts.utc(2020, 4, 28, 11, 44)

topos = api.Topos('40.8939 N', '83.8917 W')
p = topos.itrf_xyz().km

# Build satellite array.
sats = api.load.tle('http://celestrak.com/NORAD/elements/active.txt', reload=True)

# sat_array = SatrecArray([sat.model for sat in sats.values()])

# TEMP - select first 5 for testing
sat_array = [sat.model for sat in sats.values()][:5]
sat_array = SatrecArray(sat_array)

# Create time range vector of interest
start = datetime.now(timezone.utc)
start = start.replace(microsecond=0)
end = datetime.now(timezone.utc) + timedelta(hours=24)
end = end.replace(microsecond=0)
total_secs = (end - start).total_seconds()
timepoints = [start + timedelta(seconds=x) for x in range(int(total_secs) + 1)]

julian_timepoints = [jday(t.year, t.month, t.day, t.hour, t.minute, t.second) for t in timepoints]
jd, fr = zip(*julian_timepoints)

# TEMP - Select first timepoint for testing
jd_sample, fr_sample = np.array([jd[0]]), np.array([fr[0]])

# Propagate sats
e, r, v = sat_array.sgp4(jd_sample, fr_sample)

# r shape = v shape = n satellites × 1 time × 3 coordinates

r = r[:,0]
v = v[:,0]

# r shape = v shape = n satellites × 3 coordinates

r = r.swapaxes(0, 1)
v = v.swapaxes(0, 1)

# r shape = v shape = 3 coordinates × n satellites

# Convert to ITRF.

r, v = TEME_to_ITRF(t.ut1, r, v)

# Offset the satellite vectors so they are centered on our Earth
# position rather than on the Earth's center.

r = r - p[:,None]

# Since the dot product of two vectors reflects not only the angle
# between them but also their magnitudes, we first turn our vectors into
# unit vectors.

r_unit = r / length_of(r)
p_unit = p / length_of(p)

zenith_angle_cosine = dots(r_unit, p_unit[:,None])
zenith_angle = np.arccos(zenith_angle_cosine)
altitude = api.tau / 4.0 - zenith_angle
altitude_degrees = api.Angle(radians=altitude).degrees

# Sort satellites by their altitude to check our result.

for i in np.argsort(altitude_degrees):
    sat = sats[i]
    alt = altitude_degrees[i]
    print('%6.2f  %5d  %s' % (alt, sat.model.satnum, sat.name))

The error I am getting is:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-13-a5a772cece94> in <module>
     52 # Convert to ITRF.
     53 
---> 54 r, v = TEME_to_ITRF(t.ut1, r, v)
     55 
     56 # Offset the satellite vectors so they are centered on our Earth

~/anaconda3/envs/scripts/lib/python3.7/site-packages/skyfield/sgp4lib.py in TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp, yp, fraction_ut1)
    338     else:
    339         rPEF = mxv(R, rTEME)
--> 340         vPEF = mxv(R, vTEME) + _cross(angular_velocity, rPEF)
    341 
    342     if xp == 0.0 and yp == 0.0:

~/anaconda3/envs/scripts/lib/python3.7/site-packages/skyfield/sgp4lib.py in _cross(a, b)
    324 def _cross(a, b):
    325     # Nearly 4x speedup over numpy cross(). TODO: Maybe move to .functions?
--> 326     return a[_cross120] * b[_cross201] - a[_cross201] * b[_cross120]
    327 
    328 def TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp=0.0, yp=0.0, fraction_ut1=0.0):

ValueError: operands could not be broadcast together with shapes (3,) (3,5)

An associated question, it seems like the Topos object is being deprecated and the preferred way to create geographical coords is to use something like api.wgs84.latlon? I noticed that the itrs coordinates are ever so slightly off from the itrf coordinates (I assume they're different implementations of that coordinate system). Is there a way to go from itrs to teme in order to avoid calling TEME_to_ITRF in the first place?