get_orbit_number not internally consistant

[travisrlh]

Code Sample, a minimal, complete, and verifiable piece of code

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyorbital.orbital

line1 = '1 41855U 16067H   20154.58226637 +.00000953 +00000-0 +87667-4 0  9999'
line2 = '2 41855 097.9270 261.6728 0008797 306.0913 053.9493 14.97197363194247'

orb = pyorbital.orbital.Orbital('sat', line1=line1, line2=line2)

iter_days = pd.date_range(start='2020-01-01', end='2022-01-01', freq='1D')
an_timestamps = pd.to_datetime([orb.get_last_an_time(timestamp) for timestamp in iter_days]).round(freq='10S').unique()
orbit_number = np.array([orb.get_orbit_number(an_timestamp, as_float=True) for an_timestamp in an_timestamps])
lonlatalt =  np.array([orb.get_lonlatalt(an_timestamp) for an_timestamp in an_timestamps])
plt.plot(orbit_number%1,'.-')
plt.plot(lonlatalt[:,1],'.-')
plt.show()

Problem description

when using a time from get_last_an_time, get_orbit_number should always return a value that is close to an integer value. Orbit number crosses from *.999 to *.001 only as it's passing across the equator.

Expected Output

lonlatalt[:,1] = a value that is epsilon from 0 degrees latitude
orbit_number %1 = an value that is epsilon from 0 or 1 of an orbit, depending on the an time just before or after the object passes the equator

Actual Result, Traceback if applicable

lonlatalt[:,1] = a value that is epsilon from 0, as expected.
orbit_number%1 = can be any value between 0 and 1. this is unexpected

if latitude is close to 0, the orbit_number should be close to an integer value. For most TLEs this is true but the occasional TLE (one in example above) it can drift between 0 and 1 unexpectedly

Versions of Python, package at hand and relevant dependencies

Python 3.8.5 (default, Jul 28 2020, 12:59:40)
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

pyorbital.version.version_json
('\n'
'{\n'
' "date": "2020-06-24T13:56:04+0200",\n'
' "dirty": false,\n'
' "error": null,\n'
' "full-revisionid": "3ff7a6b1631deea9aed6ef227017e74f3ebdc379",\n'
' "version": "1.6.0"\n'
'}\n')

[mraspaud]

First, thanks for the nice report, with code to reproduce, and even a plot!
My first reaction to this is that TLEs aren't supposed to be reliable more than a week ahead, so I'm not surprised the quality degrades after a couple of years. That's a limitation of the sgp4 model I think.
Would that explain the error?

[travisrlh]

I agree TLE data gets stale quickly and the coordinates themselves will not be accurate, but the orbit number and coordinate should agree with one another.

it looks like get_orbit_number doesn't use sgp4, it's using a polynomial to approximate. Here's another plot showing the same issue within a few days of TLE timestamp: where latitude < 0, the orbit_number %1 should be close to 1, where latitude is > 0, orbit_number%1 should be close to 0. For day 0 and 1 this is accurate, but for day 3 there is 1 second where the satellite has crossed the equator, latitude is >0, but orbit_number has not incremented. By day 9 there are 8 seconds where this happens as well.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyorbital.orbital

line1 = '1 44370U 19037F   20001.32119615  .00004603  00000-0  11727-3 0  9993'
line2 = '2 44370  45.0162 334.0694 0014940 280.2831  79.6372 15.40966114 28640'
orb = pyorbital.orbital.Orbital('sat', line1=line1, line2=line2)

tle_timestamp = pd.to_datetime('2020-01-01')
for days in range(10):
    an_timestamp = pd.to_datetime(orb.get_last_an_time(tle_timestamp + pd.to_timedelta(days*86400e9)))
    cross_an_timestamps = pd.date_range(start=an_timestamp-pd.to_timedelta(60e9), end=an_timestamp+pd.to_timedelta(60e9), freq='1S')
    orbit_number = np.array([orb.get_orbit_number(an_timestamp, as_float=True) for an_timestamp in cross_an_timestamps])
    lonlatalt =  np.array([orb.get_lonlatalt(an_timestamp) for an_timestamp in cross_an_timestamps])
    plt.plot(lonlatalt[:,1], orbit_number%1, '.-', alpha=0.5, label=days)

plt.title('crossing the ascending node')
plt.ylabel('get_orbit_number % 1')
plt.xlabel('latitude')
plt.legend(title='TLE age (days)')
plt.grid()
plt.show()

[mraspaud]

Really interesting! thanks for digging into this.

I agree that the orbit number and coordinate should be consistent. Do you feel like giving a shot at a PR to fix this?

[travisrlh]

I don't have enough knowledge of SGP4 internals to make a fix at this time, but when I do I can check back on this issue

[mraspaud]

Thanks anyway! I'll leave this open in case any of the maintainers has the time to fix it.