Runtime warning

[mworion]

Hi, during the usage of skyfield 1.39 I get regular a runtime warning:

/Users/mw/PycharmProjects/Envs/mw4/lib/python3.8/site-packages/skyfield/functions.py:95: RuntimeWarning: invalid value encountered in double_scalars lon = arctan2(y, x) % tau

Does somebody has some hints where to look as it is am internal function and I use skyfield only from api level. Many thanks in advance, Michel

[Bernmeister]

If you can supply a short working script which reproduces the warning, that would help to diagnose.

[brandon-rhodes]

Also, if you can tell NumPy to raise an exception rather than merely print a warning, which will make Python print a full traceback so we can see which part of Skyfield was calling that function:

https://stackoverflow.com/a/15934081/85360

[mworion]

Thanks for the hints. @brandon-rhodes: I thought the stack overflow would be the Swiss army knife for it. But I still did not get the exception to raise. What I found out so far is that the reason for this runtime warning is right as I see (just having a print statement - sorry for the simple bug tracking) put in is the values of r in the function are all 'nan'. And this happens, when I process satellites from the Celestak database: active.txt. The critical satellite is 'STARLINK-1914'. It's in the database, but Heavensabove do not show this satellite in their database. I found some code for reproducing it easy:

from skyfield.api import load, wgs84, EarthSatellite

def findSatUp(sat, loc, tStart, tEnd, alt):
    t, events = sat.find_events(loc, tStart, tEnd, altitude_degrees=alt)
    return t, events

def main():
    ts = load.timescale()
    tStart = ts.now()
    tEnd = ts.tt_jd(tStart.tt + 12 / 86400)
    loc = wgs84.latlon(latitude_degrees=49, longitude_degrees=-11)
    tle = ["STARLINK-1914",
           "1 47180U 20088BL  21303.19708368  .16584525  12000-4  30219-2 0  9999",
           "2 47180  53.0402 223.8709 0008872 210.0671 150.2394 16.31518727 52528"]
    sat = EarthSatellite(tle[1], tle[2],  name=tle[0])
    t, events = findSatUp(sat, loc, tStart, tEnd, 10)
    print(t, events)

if __name__ == "__main__":
    main()

Nothing found out so far. Thanks for the help. Michel

[brandon-rhodes]

I thought the stack overflow would be the Swiss army knife for it. But I still did not get the exception to raise.

Thanks very much for the working example code! I was able to switch NumPy into “exceptions mode” by putting this at the top of your script:

import numpy as np
np.seterr(all='raise')

The resulting exception on my system, for the record, is:

Traceback (most recent call last):
  File "tmp653.py", line 28, in <module>
    main()
  File "tmp653.py", line 23, in main
    t, events = findSatUp(sat, loc, tStart, tEnd, 10)
  File "tmp653.py", line 7, in findSatUp
    t, events = sat.find_events(loc, tStart, tEnd, altitude_degrees=alt)
  File "/home/brandon/skyfield/skyfield/sgp4lib.py", line 250, in find_events
    tmax, altitude = find_maxima(t0, t1, altitude_at, half_second, 12)
  File "/home/brandon/skyfield/skyfield/searchlib.py", line 163, in find_maxima
    left, right = _choose_brackets(y)
  File "/home/brandon/skyfield/skyfield/searchlib.py", line 184, in _choose_brackets
    indices = flatnonzero(dsd < 0)
FloatingPointError: invalid value encountered in less

The problem seems to be that the event-searching routine doesn't know what to do with nan values — which you are getting, it appears, by asking about a time too far from the satellite’s epoch. If you ask for its position on the epoch itself, you’ll get a non-NaN value:

    print(sat.epoch.utc_strftime())
    print(sat.at(sat.epoch).xyz)

Which gives:

2021-10-30 04:43:48 UTC
[-3.17399331e-05 -3.03775630e-05  2.21602777e-07] au

Drawing a simple graph of the satellite's x coordinate from its epoch out to your tEnd value, it looks like it devolves to nan after about one day:

(I've forced the NaNs there to zero, so they show a flat line.)

I’m not sure at this point that I can easily diagnose in any more detail, without going into the underlying SGP4 routine and adding printf() everywhere to see exactly where its math gets into trouble and can no longer predict the position of this satellite. I can only point out that this is always a danger with any satellite being propagated for times besides its actual epoch time, and that you might want to combine the NumPy code for turning on exceptions, with your own exception handler to skip that satellite or return a placeholder value, if you want to avoid unsightly error messages.

Or maybe you could minimally ask for the satellite's x,y,z position at tStart and tEnd, and skip the satellite if either is NaN? That wouldn't defend you from a few NaN's in the middle, but it would defend against this more common case where the satellite does a hard-fail a certain number of hours or days from the epoch.

Oh — and, recall that you can ask for the position’s .message to see why the SGP4 routine was upset. In this case:

mean eccentricity is outside the range 0.0 to 1.0

[mworion]

@brandon-rhodes, Many thanks for the explanations. I really appreciate your time investment as I take many insights out of them. I have to apologize about the missing stack trace. I did it exactly the same way in my overall application, but did not see anything and worked back and forth. For the code snippet I just did not try it again, but I should have done this.

As with the use of the satellite features of my application it will process many satellite data from various databases and time windows, I need to add some robustness for selecting them for further work. I will take some of the ideas back and try to find some solutions for it. Again, thanks for the help and closing the issue. Michel

[brandon-rhodes]

I'm glad that we at least worked out why the search routine wasn't working! The next time I'm in that code, I'll try to remember the NaN position, and see whether the search routine maybe could fail with a more informative error message if it’s given NaN coordinates to work with.

In the meantime, I hope you have success getting your satellite routines armored against this failure mode!

[mworion]

I'm working on that. Your idea checking the tEnd might have helped. We both observed a satellite problem:

From 04-Nov-21: 
https://orbit.ing-now.com/satellite/71914/2020-088bl/starlink-1914/ 

STARLINK-1914/2020-088BL Reentry 4 days ago!

So your chart showed the crash of one of these satellites.