satellite.find_events is missing events

[SHellmich]

Hi!

I just realized that find_events() misses the setting event when t0 is set to a time just after culmination but before setting.

from datetime import timedelta

from astropy.time import Time
from skyfield.api import EarthSatellite, load, wgs84

ts = load.timescale()

line1 = "1 39191U 13031D   21365.68950013 -.00000013  00000-0  00000-0 0  9995"
line2 = "2 39191 000.0397 004.0913 0002586 278.0623 077.8173 05.00115674155430"
satellite = EarthSatellite(line1, line2, ts=ts)

observer = wgs84.latlon(-24.626331, -70.403964, 2369.34)

t_1 = ts.from_astropy(Time("2022-01-02T00:19:53.630", format="isot", scale="utc"))
t_2 = ts.from_astropy(Time("2022-01-02T05:19:53.630", format="isot", scale="utc"))

t, events = satellite.find_events(observer, t_1, t_2, altitude_degrees=30)
print(t, events)

for t_i, event_i in zip(t, events):
    if event_i == 1:
        t, events = satellite.find_events(
            observer, t_i + timedelta(seconds=0.1), t_2, altitude_degrees=30
        )
        print(t, events)

Am I doing something wrong?

Some wider context of what I'm up to:
I'm currently working on a script that searches for potentially observed satellites in astronomical images. The goal is to compute the precise time when an object enters and leaves the field. Since I need to check for a large number of satellites, I'm trying to filter the candidates before I search for the precise observation time. I tried with find_minima() but it turned out that for satellites below the horizon, the (satellite - observer).at(t).altaz() returns nan values, so I tried with observe() but that made the program crash (It seems there is a memory issue, I'll try to make an example to reproduce this and open another issue). As a workaround, I'm trying to reduce the number of candidates by searching those that are around the pointed altitude within the observation time, but this unfortunately misses objects that have already culminated at the start of the exposure...

[EndlessDex]

Uh. altaz does not produce NaN when below horizon?

import skyfield.api as sf

ts = sf.load.timescale()

line1 = "1 39191U 13031D   21365.68950013 -.00000013  00000-0  00000-0 0  9995"
line2 = "2 39191 000.0397 004.0913 0002586 278.0623 077.8173 05.00115674155430"
satellite = sf.EarthSatellite(line1, line2, ts=ts)

observer = sf.wgs84.latlon(-24.626331, -70.403964, 2369.34)

SPAN_HOURS = 5
RESOLUTION_MINUTES = 30
timespan = ts.utc(2022, 1, 2, 0, range(0, 60 * SPAN_HOURS, RESOLUTION_MINUTES))
topocentric = (satellite - observer).at(timespan)
alt, az, dist = topocentric.altaz()

positions = list(zip(timespan, alt.degrees, az.degrees))
print(*positions, sep='\n')

(<Time tt=2459581.500800741>, -71.23754483418848, 154.90586068966513)
(<Time tt=2459581.521634074>, -68.52840940036559, 219.3380915094899)
(<Time tt=2459581.5424674074>, -52.28736970797745, 249.93179972242018)
(<Time tt=2459581.563300741>, -31.99605286574774, 265.2578019863132)
(<Time tt=2459581.584134074>, -8.681496693825755, 278.0790936068631)
(<Time tt=2459581.6049674074>, 18.56151577764979, 295.53690498335936)
(<Time tt=2459581.625800741>, 44.96484143086254, 334.90575913142794)
(<Time tt=2459581.646634074>, 39.3190793071871, 39.52207633404662)
(<Time tt=2459581.6674674074>, 10.974153945616925, 70.08777838124449)
(<Time tt=2459581.688300741>, -15.16793472645499, 85.39693494677331)

[EndlessDex]

Regarding your usage of find_events:

The function specifically looks for local maxima. You cut out the only local maxima in the timespan. So it is working correctly when it says that there are no events.

[brandon-rhodes]

I tried with find_minima() but it turned out that for satellites below the horizon, the (satellite - observer).at(t).altaz() returns nan values…

Any chance you could share a small example where nan is returned, for me to investigate? Or is it not likely that you could quickly re-create that case?

[SHellmich]

I'll try to create one. It indeed works with the code I shared (sorry for not testing for that). Either there is some other bug in my code or the issue is more difficult to reproduce.

Thanks for the comment on find_events! Is there a method that can find the setting below a certain altitude after culmination would the solution be to simply test if the satellite is up and search for the next setting manually?

[EndlessDex]

Regarding your stated goal, this is how I would do it. Modify the variables to meet your needs and BAM!

I find 1% of orbit period hits the sweetspot of fast and effective, but you can definitely adjust depending on the size of your FOV

import skyfield.api as sf
from skyfield.searchlib import find_discrete
import numpy as np

# Setup
ts = sf.load.timescale()

line1 = "1 39191U 13031D   21365.68950013 -.00000013  00000-0  00000-0 0  9995"
line2 = "2 39191 000.0397 004.0913 0002586 278.0623 077.8173 05.00115674155430"
satellite = sf.EarthSatellite(line1, line2, ts=ts)

observer = sf.wgs84.latlon(-24.626331, -70.403964, 2369.34)
# Setup End

# Variables
ALT_CENTER = 45
ALT_RADIUS = 5
AZ_CENTER = 340
AZ_RADIUS = 5
PERCENT_OF_PERIOD = 0.01
T_START = ts.utc(2022, 1, 2, 0)
T_STOP = ts.utc(2022, 1, 2, 5)
# Variables End

# Algorithm
def satellite_in_view(t):
    alt, az, _ = (satellite - observer).at(t).altaz()
    results = np.array([])
    for alt_i, az_i in zip(alt.degrees, az.degrees):
        in_alt_view = (ALT_CENTER - ALT_RADIUS) < alt_i and alt_i < (ALT_CENTER + ALT_RADIUS)
        in_az_view = (AZ_CENTER - AZ_RADIUS) < az_i and az_i < (AZ_CENTER + AZ_RADIUS)
        results = np.append(results, in_alt_view and in_az_view)
    return results

rad_per_min_TO_rev_per_day = (60.0 * 24.0) / (2 * np.pi)
period_of_orbit = 1.0 / (satellite.model.no_kozai * rad_per_min_TO_rev_per_day)
satellite_in_view.step_days = period_of_orbit * PERCENT_OF_PERIOD

t, values = find_discrete(T_START, T_STOP, satellite_in_view)
# Algorithm End

# Output
for t_i, v_i in zip(t, values):
    alt, az, _ = (satellite - observer).at(t_i).altaz()
    print(t_i.utc_iso(), "ENTERING" if v_i else "EXITING ", alt.degrees, az.degrees)
# Output End

2022-01-02T02:57:31Z ENTERING 43.42015241974112 330.0000009069733
2022-01-02T03:06:50Z EXITING  47.79899187616711 350.0000073975163

[SHellmich]

Thanks for the example! I used find_minima but find_discrete is indeed much more elegant.

I'm having a hard time to reproduce the memory issue with a minimal example. I really don't know what caused it but in my actual program, at some point, tracemalloc showed that memory consumption in skyfield/iokit.py increased. However, I think this was caused by calling load.timescale() and load("de421.bsp") multiple times. After cleaning the code a bit, skyfield memory consumption did not increase anymore and it turned out that rather astropy/io/fits/hdu/image.py seems to cause the problems. I think that altaz() returning nan values might have been a result of memory was already corrupted and lead to unpredictable behavior just before the program eventually crashed.

I'll redesign my code to use find_discrete, as you suggested and report back here (or at astropy) if there are still problems. Thanks a lot for the support! And thanks for skyfield! I makes things much easier and elegant!

[lpmrfentazis]

@EndlessDex
@brandon-rhodes
Hi, I'm facing a similar problem, but I don't seem to understand your comment correctly.
The way I see it, find_events should provide any events that occur between t0 and t1. In my case, if the satellite is below horizon at t0, but before t1 it rises above horizon (in this example at t1 the culmination). But I am observing empty list.

I apologize for almost the same question that a colleague asked, but I seem to have lost the gist of the answer due to translation.

Here is my working example:

from skyfield.api import load, Topos, EarthSatellite, N, E
from datetime import datetime, timedelta
import pytz

# ELEKTRO-L 2 TLE from Celestrack 26.08.2024 8:34 UTC
line1 = "1 41105U 15074A   24238.84268576 -.00000128  00000+0  00000+0 0  9993"
line2 = "2 41105   5.0153  78.6193 0002491 153.4123  31.4253  1.00270890 31881"


satellite = EarthSatellite(line1=line1, line2=line2, name="ELEKTRO-L 2")
observer = Topos(55.671429 * N, 37.62539 * E, elevation_m=180)

ts = load.timescale()
start = datetime.fromtimestamp(1724661467.610997, tz=pytz.utc)  # 26.08.2024 8:38 UTC
finish = start + timedelta(hours=5)

horizon = 14

diff = satellite - observer

times, events = satellite.find_events(observer, 
                                      t0=ts.from_datetime(start), 
                                      t1=ts.from_datetime(finish), 
                                      altitude_degrees=horizon)

# events empty
print(events)

# but

pos = diff.at(ts.from_datetime(start)).altaz()
print(f"finish el={pos[0]}")

pos = diff.at(ts.from_datetime(finish)).altaz()
print(f"finish el={pos[0]}")

Python 3.11.3
Skyfield version: 1.46

[lpmrfentazis]

If you make finish = start + 6 hours, it works about as I expected when the satellite starts to descend below the culmination.

[lpmrfentazis]

[SHellmich]

This behavior is a bit unexpected but still intended. If I understand correctly, find_events() will do a find_maxima() to search for culminations and then find the rising and setting events around these maxima. If the search interval does not include a culmination, there is no maximum to be found by find_maxima(), so nothing will be returned.

[EndlessDex]

@lpmrfentazis

The comment by @SHellmich is correct.

find_events specifically targets culminations over a certain altitude. Here is a much reduced version of the implementation of find_events.

Find Events

def find_events(self, topos, t0, t1, altitude_degrees=0.0):
        # First, we find the moments of maximum altitude over the time
        # period.  Some of these maxima will be negative, meaning the
        # satellite failed to crest the horizon.
        def altitude_at(t):
            return at(t).altaz()[0].degrees
        tmax, altitude = find_maxima(t0, t1, altitude_at)

        # Next, filter out the maxima that are not high enough.
        keepers = altitude >= altitude_degrees
        jdmax = tmax.tt[keepers]

        # Finally, find the rising and setting that bracket each maximum
        # altitude.  We guess that the satellite will be back below the
        # horizon in between each pair of adjancent maxima.
        def below_horizon_at(t):
            return at(t).altaz()[0].degrees < altitude_degrees
        trs, rs = find_discrete(****, below_horizon_at)

@brandon-rhodes potentially the description of find_events could be clarified here?

Perhaps:

Searches between t0 and t1, which should each be a Skyfield Time object, for passes of this satellite above the location topos that reach culminate at least altitude_degrees above the horizon. A pass is defined by at least one culmination of the satellite.