Searching along an orbit track

[gfireman]

I am trying to use skyfield.searchlib to search for specific solar angles along a EarthSatellite ground track. I've been able to find the ground track, and to find angles for a specific point, but not in putting the two together.

Here is my code for a simple demo:

#!/usr/bin/env python3
import numpy as np
from matplotlib import pyplot as plt
from skyfield.api import EarthSatellite, load, wgs84
from skyfield.searchlib import find_discrete, find_maxima

eph = load("de440s.bsp")
sun = eph["sun"]
earth = eph["earth"]
ts = load.timescale()

def fday(t, offset=0):  # time as fractional days
    s = np.array([dt.timestamp() for dt in t.utc_datetime()])
    return (s - offset.utc_datetime().timestamp()) / 86400

def altitude(t):  # solar altitude in degrees
    location = (earth + wgs84.latlon(lat, lon)).at(t)
    alt, _, _ = location.observe(sun).apparent().altaz()
    return alt.degrees

def horizon(t):  # when sun crosses nominal horizon
    return altitude(t) > 0

def rundemo():
    nresults = int(ndays / step_days + 0.5)

    # init plot
    plt.ion()
    fig, ax = plt.subplots(figsize=(5, 2))
    ax.grid()
    ax.axhline(0, color="b")
    ax.set(
        title="Solar Altitude",
        xlabel=f'days since {t[0].utc_strftime("%Y-%m-%d")}',
        ylabel="degrees",
    )

    # solar altitude in degrees
    alt = altitude(t)
    ax.plot(fday(t, offset=t[0]), alt)
    fig.tight_layout()
    fig.savefig(ofile)

    # find where sun crosses horizon
    horizon.step_days = step_days
    htimes, _ = find_discrete(t[0], t[-1], horizon, num=nresults)
    horalt = np.zeros_like(htimes)
    ax.plot(fday(htimes, offset=t[0]), horalt, "ko")
    fig.tight_layout()
    fig.savefig(ofile)

    # find where sun is at max altitude
    altitude.step_days = step_days
    mtimes, maxalt = find_maxima(t[0], t[-1], altitude, num=nresults)
    ax.plot(fday(mtimes, offset=t[0]), maxalt, "r.")
    fig.tight_layout()
    fig.savefig(ofile)

def times(ndays):  # define times in 10-minute intervals
    nvals = ndays * 1440 + 1
    return ts.utc(2024, 7, 10, 0, range(0, nvals, 10))

Specifying a single point works fine, whether or not coordinates are array-like:

# parameters
ofile = "plot1.png"
ndays = 7
step_days = 0.25
t = times(ndays)

# define observing position
lat = np.array([40.8939])
lon = np.array([-83.8917])

But when I search along a path, I get an error.

# parameters
ofile = "plot2.png"
ndays = 1
step_days = 1 / 30  # avg 14.5 orbits/day
t = times(ndays)

# define observing position as subsatellite coordinates
sat = EarthSatellite(
    "1 58928U 24025A   24192.40949872  .00000794  00000-0  15498-3 0  9996",
    "2 58928  98.0945 123.6525 0001014  97.0736 263.0583 14.65602003 22425",
)
lat, lon = wgs84.latlon_of(sat.at(t))
lat = lat.degrees
lon = lon.degrees

Traceback (most recent call last):
  File "./demo2.py", line 91, in <module>
    rundemo()
  File "./demo2.py", line 53, in rundemo
    htimes, _ = find_discrete(t[0], t[-1], horizon, num=nresults)
  File "/opt/homebrew/lib/python3.8/site-packages/skyfield/searchlib.py", line 44, in find_discrete
    return _find_discrete(ts, jd, f, epsilon, num)
  File "/opt/homebrew/lib/python3.8/site-packages/skyfield/searchlib.py", line 56, in _find_discrete
    y = f(t)
  File "./demo2.py", line 26, in horizon
    return altitude(t) > 0
  File "./demo2.py", line 20, in altitude
    location = (earth + wgs84.latlon(lat, lon)).at(t)
  File "/opt/homebrew/lib/python3.8/site-packages/skyfield/vectorlib.py", line 91, in at
    p, v, gcrs_position, message = self._at(t)
  File "/opt/homebrew/lib/python3.8/site-packages/skyfield/vectorlib.py", line 215, in _at
    p2, v2, _, message = vf._at(t)
  File "/opt/homebrew/lib/python3.8/site-packages/skyfield/toposlib.py", line 48, in _at
    r = mxv(RT, r)
  File "/opt/homebrew/lib/python3.8/site-packages/skyfield/functions.py", line 38, in mxv
    return einsum('ij...,j...->i...', M, v)
  File "<__array_function__ internals>", line 200, in einsum
  File "/opt/homebrew/lib/python3.8/site-packages/numpy/core/einsumfunc.py", line 1371, in einsum
    return c_einsum(*operands, **kwargs)
ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (3,3,32)->(3,32,3) (3,145)->(145,3)

This error seems to be related to issue #539.
I would be grateful if someone could enlighten me as to the correct method.

[brandon-rhodes]

Interesting! I've never tried combining $n$ longitudes and latitudes with $n$ times, but, happily, it looks like the operation works fine in Skyfield:

#!/usr/bin/env python3

import numpy as np
from skyfield.api import EarthSatellite, load, wgs84

ts = load.timescale()
t = ts.utc(2024, 7, 10, [1, 2, 3])

lat = np.array([40.0, 41.0, 42.0])
lon = np.array([-83.0, -84.0, -85.0])

track = wgs84.latlon(lat, lon)
print(track.at(t).xyz)  # multiple lat's and lon's seem to work fine

eph = load('de421.bsp')
earth = eph['earth']
location = (earth + track).at(t)  # and can be combined with multiple t's
print(location.xyz)

Which, happily, seems to work fine!

The problem seems to come in when you run routines like find_discrete() and find_maxima(), that do their work by asking about hundreds of points in between a start time and end time. What is happening in your particular traceback is that find_discrete() creates its own array of 145 times between t[0] and t[1] and wants to know the value of altitude(t) for those 145 times. But your altitude() routine takes those 145 times and tries to use .at(t) to combine them with the mere 32 latitude-longitude positions that you've provided. And since 145 ≠ 32, NumPy doesn't know how to proceed, since it's been given arrays of two different lengths.

So I think you need to upgrade your altitude() routine so that instead of being satisfied with your 32 latitude-longitude positions, it looks at the t it has been given, and computes new satellite latitude-longitude positions for that new array t. So, instead of trying to pre-compute lat and lon, try re-computing them in each call to altitude().

[gfireman]

That was it - thank you!
All I had to do was to replace the altitude function as:

def altitude(t):  # solar altitude in degrees
    lat, lon = wgs84.latlon_of(sat.at(t))
    lat = lat.degrees
    lon = lon.degrees
    location = (earth + wgs84.latlon(lat, lon)).at(t)
    alt, _, _ = location.observe(sun).apparent().altaz()
    return alt.degrees

[brandon-rhodes]

I'm glad we were able to make Skyfield work here! Thanks for circling back and sharing with future readers the code that got your script working.