Open SHellmich opened 3 weeks ago
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)
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.
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?
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?
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
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!
@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
If you make finish = start + 6 hours, it works about as I expected when the satellite starts to descend below the culmination.
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.
@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
.
@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
reachculminate at least altitude_degrees above the horizon. A pass is defined by at least one culmination of the satellite.
Hi!
I just realized that find_events() misses the setting event when t0 is set to a time just after culmination but before setting.
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 withobserve()
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...