get_next_passes returns max-elevation-time time not between rise & fall time

[johnmccombs]

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

from pyorbital import tlefile
from pyorbital.orbital import Orbital
from datetime import datetime

TLE_FILE_NAME = 'iridium-tle.txt'

HORIZON = 40

# altitude (km) above seal-level and lat lon of the observer
ALTITUDE_ASL = 0.5
LON = 170.556
LAT = -43.368

tlefile.TLE_URLS = ('https://celestrak.com/NORAD/elements/iridium.txt', )
tlefile.fetch(TLE_FILE_NAME)

orb = Orbital('IRIDIUM 7 [+]', TLE_FILE_NAME)
d = datetime(2018, 3, 7, 3, 40, 15)
orb.get_next_passes(d, 1, LON, LAT, ALTITUDE_ASL, horizon=HORIZON)

Problem description

The max elevation time for a pass should lie between the rise and set times.

Expected Output

The documentation says get_next_passes returns " [(rise-time, fall-time, max-elevation-time), …]"

Max elevation time should lie between
003:46:14 and 03:50:12

Actual Result, Traceback if applicable

The max-elevation-time is 04:38:53, almost an hour after the set time.

[(datetime.datetime(2018, 3, 7, 3, 46, 14, 521170), datetime.datetime(2018, 3, 7, 3, 50, 12, 734428), datetime.datetime(2018, 3, 7, 4, 38, 53, 847556))]

Versions of Python, package at hand and relevant dependencies

Python 2.7.13 and 3.6.4. Pyorbital 1.2.0

Thank you for reporting an issue !

[mraspaud]

@johnmccombs thanks for spotting this. I believe I have a fix implemented in PR #23, would you mind checking if that works for you ? I tested against your minimal example, and it worked.

[johnmccombs]

Agreed that fixes my example, however I found some Iridium passes that go directly overhead, but are only 1.4 minutes long, which doesn't make sense. And also instances of the transit time well outside the rise/set times. Repro code below. I attached my current TLE file so we're using the same data.

I did a plot (attached) of duration vs transit elevation of a week's worth of calculations and you can see there are a number of high passes with short duration. The negative elevations are caused by the transit time be much later than the rise/set. I have a script that does all the passes for all the Iridium satellites for a week, if that would help.

iridium-tle.txt

from pyorbital import tlefile
from pyorbital.orbital import Orbital
from datetime import datetime

TLE_FILE_NAME = 'iridium-tle.txt'

HORIZON = 40

# altitude (km) above seal-level and lat lon of the observer
ALTITUDE_ASL = 0.5
LON = 170.556
LAT = -43.368

#tlefile.TLE_URLS = ('https://celestrak.com/NORAD/elements/iridium.txt', )
#tlefile.fetch(TLE_FILE_NAME)

def do_calc(d, sat_name):
    orb = Orbital(sat_name, TLE_FILE_NAME)
    passes = orb.get_next_passes(d, 1, LON, LAT, ALTITUDE_ASL, horizon=HORIZON)
    rise_az, rise_el = orb.get_observer_look(passes[0][0], LON, LAT, ALTITUDE_ASL)
    transit_az, transit_el = orb.get_observer_look(passes[0][2], LON, LAT, ALTITUDE_ASL)
    set_az, set_el = orb.get_observer_look(passes[0][1], LON, LAT, ALTITUDE_ASL)

    print(sat_name)
    print(passes[0][0], rise_az, rise_el)
    print(passes[0][2], transit_az, transit_el)
    print(passes[0][1], set_az, set_el)
    print()

do_calc(datetime(2018, 3, 9, 17, 35, 0), 'IRIDIUM 40 [-]')
do_calc(datetime(2018, 3, 8, 23, 21, 0), 'IRIDIUM 47 [+]')

[mraspaud]

Very interesting, thanks a lot for the investigation ! Looks like we need to check the horizon crossing times.

[mraspaud]

@johnmccombs Ok, I just had a look at this again, and I tested with multiple satellites, and fresh TLE data from celestrak, with passes over a period of one month. So far, I could only observe the problem you mention with Iridium 40, the others seem to have the expected behaviour.

Hence I looked further and saw that the orbit of Iridium 40 is much different from the other iridiums, and its orbit altitude seem to oscillate between 250 and 500 km (the other iridiums are stable around 750km). Hence, making a similar plot to yours, with the altitudes giving the color to the points, I get the following:

So to me, it seems logical that lower orbits have smaller pass durations. Or am I missing something here ?

[johnmccombs]

You're right about the short ones. I took a look at Iridium 40 and the perigee is almost overhead at my chosen location, so a brief overhead pass is reasonable.

The other issue though at the right-hand end of my graph with negative transit elevations. I added some code to brute-force the transit. The script below shows the brute-forced values on the second line of output and the get_next_passes() value on the third line. Looks like the transit calculation in get_next_passes() occasionally gets lost and returns the wrong time. The script is setup to show an example of this. Hope that helps.

from pyorbital import tlefile
from pyorbital.orbital import Orbital
from datetime import datetime
from datetime import timedelta
import re

TLE_FILE_NAME = 'iridium-tle.txt'

HORIZON = 40

# altitude (km) above sea-level and lat lon of the observer
ALTITUDE_ASL = 0.5
LON = 170.556
LAT = -43.368

SATELLITE_LIST = []

tlefile.TLE_URLS = ('https://celestrak.com/NORAD/elements/iridium.txt', )
tlefile.fetch(TLE_FILE_NAME)


def get_satellite_names():
    # read the element file to get a list of names for the satellites
    with open(TLE_FILE_NAME, 'rt') as f:
        for l in f:
            if not re.match('^[12] ', l) and not l.startswith('DUMMY MASS'):
                name = l.strip()
                SATELLITE_LIST.append(name)


def timerange(start_date, end_date):
    for n in range(int ((end_date - start_date).seconds)):
        yield start_date + timedelta(0, n)


def find_transit(orb, t0, t1):

    max_el = -999
    max_az = -999
    t_max  = -999

    # hack to handle case where SV just touches the horizon
    if t0 + timedelta(0, 2) < t1:

        for t in timerange(t0, t1):
            az, el = orb.get_observer_look(t, LON, LAT, ALTITUDE_ASL)
            if el > max_el:
                max_el = el
                max_az = az
                t_max = t

        return t_max, max_az, max_el

    else:
        return t0, orb.get_observer_look(t0, LON, LAT, ALTITUDE_ASL)
  

def do_calc(d, sat_name, duration):
    orb = Orbital(sat_name, TLE_FILE_NAME)
    passes = orb.get_next_passes(d, duration, LON, LAT, ALTITUDE_ASL, horizon=HORIZON)

    for p in passes:
        rise_az, rise_el = orb.get_observer_look(p[0], LON, LAT, ALTITUDE_ASL)
        transit_time, transit_az, transit_el = find_transit(orb, p[0], p[1])
        set_az, set_el = orb.get_observer_look(p[1], LON, LAT, ALTITUDE_ASL)

        tr_az, tr_el = orb.get_observer_look(p[2], LON, LAT, ALTITUDE_ASL)

        alert = ' <<<<<<<<<<<<<' if tr_el < 0.01 or not (p[2] > p[0] < p[1]) else ''

        print(sat_name)
        print(p[0], rise_az, rise_el)
        print(transit_time, transit_az, transit_el)
        print(p[2], tr_az, tr_el, alert)
        print(p[1], set_az, set_el)
        print()

do_calc(datetime(2018, 3, 12, 14, 30, 0), 'IRIDIUM 7 [+]',  2)

# get_satellite_names()
# for s in SATELLITE_LIST:
#     try:
#         do_calc(datetime(2018, 3, 9, 17, 35, 0), s, 96)
#     except:
#         pass

[grypoB]

I was able to reproduce the bug with this short script:

from datetime import datetime, timedelta
import pyorbital.orbital
import numpy as np
from scipy import optimize

tle_lines = [
    "",
    "1 44083U 19018F   21029.41783953  .00001194  00000-0  55122-4 0  9996",
    "2 44083  97.3931  93.0075 0014683  99.5364 260.7532 15.22064008101741",
]
gs = (46.517862, 6.5602297, 0.390)

sat = pyorbital.orbital.Orbital("", line1=tle_lines[1], line2=tle_lines[2])
t0 = datetime(2020, 7, 19, 8)

passes = sat.get_next_passes(t0, 2, gs[1], gs[0], gs[2])

for c in passes:
    az, el = sat.get_observer_look(c[2], gs[1], gs[0], gs[2])
    print([t.isoformat() for t in c], f"Max el: {el:.4}")

    if c[2] > c[1] or c[2] < c[0]: # max-elevation time not between rise and fall time
        c = list(c)
        c[2] = c[0] + (c[1] - c[0]) / 2
        print("Max-elevation time not within bounds, try mid-point approx")
        c = tuple(c)
        az, el = sat.get_observer_look(c[2], gs[1], gs[0], gs[2])
        print("correction", [t.isoformat() for t in c], f"Max el: {el:.4}")

Which gives the following output:

['2020-07-19T08:05:14.306864', '2020-07-19T08:13:14.786040', '2020-07-19T08:09:15.019480'] Max el: 6.819
['2020-07-19T09:37:50.372532', '2020-07-19T09:49:20.394733', '2020-07-19T11:16:43.327782'] Max el: 8.588
Max-elevation time not within bounds, try mid-point approx
correction ['2020-07-19T09:37:50.372532', '2020-07-19T09:49:20.394733', '2020-07-19T09:43:35.383632'] Max el: 89.22

I was able to pin-point the issue in the function to the line 429-432 of orbital.py:

                    highest = utc_time + \
                        timedelta(minutes=get_max_parab(
                            elevation_inv,
                            max(risemins, middle - 1), min(fallmins, middle + 1),
                            tol=tol / 60.0
                        ))

And simply replacing line 432 by the following fixes the issue:

max(risemins, middle - 2), min(fallmins, middle + 2)

I'm not familiar enough with the codebase to know if this will generate new bugs, but it fixes mine at least (where the max-elev time for a high elevation pass (>89 deg) is not between rise and fall time)

@mraspaud let me know if you want me to create a PR for that fix

[mraspaud]

Thanks for finding a minimal example, that'll help writing a failing unit test!

Hmm, it's been a while, can you remind me what the unit of middle is?

[grypoB]

From my understanding of the function, `middle` should be in minutes (where 0 is the `utc_time` given in the function input). It also looks like an index, but I think that's because you have an array of elevation `elev` computed at 1 minute intervals And now that I look at the code more closely I also think the tolerance given to the get_max_parab function should be tol*60, instead of tol/60. Get_max_parab works on minutes, but tol of the get_next_passes is in seconds according to the doc

…

On Mon, Feb 1, 2021, 7:01 PM Martin Raspaud ***@***.***> wrote: Thanks for finding a minimal example, that'll help writing a failing unit test! Hmm, it's been a while, can you remind me what the unit of middle is? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#22 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACOQ4UGCDFJI5RHAY2Z4L4TS43T7TANCNFSM4ET6HZIA> .

[dagnic]

Hi !

Is there some update about this issue ?

I have some errors of this kind (example below). @grypoB how did you come to the conclusion that line 432 should be replaced by the following?

max(risemins, middle - 2), min(fallmins, middle + 2)

If it works, is it possible to implement it?

Thanks !

import pyorbital.orbital
import datetime
orb = pyorbital.orbital.Orbital("NOAA 18")
t = datetime.datetime(2021,3,9,22)
orb.get_next_passes(t, 1, -15.6335, 27.762, 0.)

which gives

Out[1]:
[(datetime.datetime(2021, 3, 9, 22, 17, 42, 127720),
  datetime.datetime(2021, 3, 9, 22, 33, 32, 873936),
  datetime.datetime(2021, 3, 9, 19, 56, 29, 161898))]

[dagnic]

So @grypoB I tested your solution, replacing line 432 by:

max(risemins, middle - 2), min(fallmins, middle + 2)

but, alone, it does not work for my example I gave above

Out[1]:
[(datetime.datetime(2021, 3, 9, 22, 17, 42, 63769),
  datetime.datetime(2021, 3, 9, 22, 33, 32, 820920),
  datetime.datetime(2021, 3, 10, 0, 6, 42, 220473))]

In order to make it work, I have to change also line 433 with tol=tol * 60.0 as you mentioned before.

Out[1]:
[(datetime.datetime(2021, 3, 9, 22, 17, 42, 63769),
  datetime.datetime(2021, 3, 9, 22, 33, 32, 820920),
  datetime.datetime(2021, 3, 9, 22, 25, 38, 125490))]

Is it possible to implement it, or at least detect that the maxtime is not within rise and falltime, then perform mid-poind approximation and warn?

[mraspaud]

@dagnic could you provide the TLE lines for that example?

[dagnic]

'1 28654U 05018A   21083.16603416  .00000102  00000-0  79268-4 0  9999'
'2 28654  99.0035 147.6583 0014816 159.4931 200.6838 14.12591533816498'

[mraspaud]

thanks!

[mraspaud]

@grypoB @dagnic I think I have a fix in #76, could you check if it works for you?

[dagnic]

@mraspaud I clone the branch, its looks good. Thanks.

[grypoB]

@mraspaud same here, it does fix the issue in my example. Thanks!

[mraspaud]

@dagnic @grypoB Thanks for the feedback! I'll merge my PR then.