pytroll / pyorbital

Orbital and astronomy computations in python
http://pyorbital.readthedocs.org/
GNU General Public License v3.0
226 stars 77 forks source link

get_next_passes() bug? #5

Closed Box88 closed 9 years ago

Box88 commented 9 years ago

Hi Martin,

I am developing an application to predict pass conflict of NOAA 19 and NOAA 18 blind orbits with METOP-A and METOP-B passes over EUMETSAT Svalbard ground station. I have been used pyephem as TLE propagator and some days ago I found out also this nice library: pyorbital. I checked the prediction results of both libraries (pyephem and pyorbital) against the data given by the EUMETSAT flight dynamics team and results from pyorbital seem more accurate. Then I would like to integrate also pyorbital in my tool but I found some unexpected results using the get_next_passes function.

Here is the code for verification:

"""Pyorbital verification"""

from pyorbital import orbital
from datetime import datetime

def tle_propagator(line1, line2, line3, gs_info, start, stop):
    """Propagates TLEs using pyorbital package"""

    dt_format = '%Y-%m-%d %H:%M:%S'
    length = int((stop - start).days*24)
    gs_name = gs_info[0]
    lat = gs_info[1]
    lon = gs_info[2]
    alt = gs_info[3]/1000 # assuming that altitude must be given in km, am I right?

    orbit = orbital.Orbital(line1 ,line1 = line2, line2 = line3)
    visibility = orbit.get_next_passes(start,length,lon,lat,alt)

    rise_time = []
    rise_az = []
    max_el_time = []
    max_el = []
    set_time = []
    set_az = []
    duration = []
    sid = []
    antenna = []
    number = []

    for vis in visibility:
        rise_time.append(vis[0].strftime(dt_format))
        rise_az.append(orbit.get_observer_look(vis[0],lon,lat,alt)[0])
        max_el_time.append(vis[2].strftime(dt_format))
        max_el.append(orbit.get_observer_look(vis[2],lon,lat,alt)[1])
        set_time.append(vis[1].strftime(dt_format))
        set_az.append(orbit.get_observer_look(vis[1],lon,lat,alt)[0])
        duration.append((vis[1]-vis[0]).total_seconds())
        number.append(orbit.get_orbit_number(vis[0], tbus_style=True))
        sid.append(line1)
        antenna.append(gs_name)

    results =  [antenna, sid, rise_time, rise_az, max_el_time, max_el,
                set_time, set_az, duration, number]

    return results

if __name__ == '__main__':

    sat = 'METOP-A'
    # TLE from Celestrack updated at 18/08/2015
    if sat == 'METOP-B':         
        line1 = '1 38771U 12049A   15229.80404600  .00000046  00000-0  40891-4 0  9999'
        line2 = '2 38771  98.6922 288.8157 0001608  64.8585  28.9446 14.21489933151176'
    elif sat == 'METOP-A':         
        line1 = '1 29499U 06044A   15229.56308789  .00000049  00000-0  42323-4 0  9993'
        line2 = '2 29499  98.6704 287.8437 0001086  58.6288  78.6540 14.21498242458004'
    elif sat == 'NOAA 19':
        line1 = '1 33591U 09005A   15229.46280543  .00000108  00000-0  83803-4 0  9991'
        line2 = '2 33591  98.9943 179.5032 0014180   2.1481 357.9748 14.11967048336132'

    start = datetime(2015,8,20,0,0,0)
    stop = datetime(2015,9,10,0,0,0)
    gs_selected = 'CDA1'
    # ground stations
    if gs_selected == 'CDA1':
        lat = 78.228981740
        lon = 15.388229269
        alt = 490.819 # in m
    elif gs_selected == 'FAIRBANKS':
        lat = 64.973760002
        lon = -147.505823340
        alt = 397.033 # in m
    elif gs_selected == 'WALLOPS':
        lat = 37.947340276
        lon = -75.462919721
        alt = -21.250 # in m

    gs = gs_selected, lat, lon, alt

    results = tle_propagator(sat,line1,line2,gs,start,stop)
    passes = []
    for row in zip(*results):
        passes.append(row)
    print passes

If you run this code you will note that one pass has the following information:

('CDA1', 'METOP-A', '2015-09-01 03:04:50', 55.468140794607727, '2015-09-01 03:50:44', -79.594594264864824, '2015-09-01 04:48:11', 88.336125650096918, 6200.204852, 46008)

which is clearly wrong, but if you then change the start to datetime(2015,9,1,0,0,0) and stop to datetime(2015,9,2,0,0,0) then you get the following correct results;

[('CDA1', 'METOP-A', '2015-09-01 01:10:39', 292.84138407089591, '2015-09-01 01:16:20', 9.3752155556005459, '2015-09-01 01:22:01', 29.390424037000855, 682.417786, 46007),
('CDA1', 'METOP-A', '2015-09-01 02:53:44', 321.72323204016112, '2015-09-01 02:59:17', 8.7303638974740849, '2015-09-01 03:04:50', 55.468140697666151, 666.76633, 46008),
('CDA1', 'METOP-A', '2015-09-01 04:36:03', 343.2786315673265, '2015-09-01 04:42:07', 11.572023051734369, '2015-09-01 04:48:11', 88.336125650096918, 727.722306, 46009),
('CDA1', 'METOP-A', '2015-09-01 06:17:41', 359.7629881651298, '2015-09-01 06:24:30', 18.483175863740684, '2015-09-01 06:31:18', 124.45885193322324, 817.012133, 46010),
('CDA1', 'METOP-A', '2015-09-01 07:58:50', 14.171452278754581, '2015-09-01 08:06:15', 31.540118963785787, '2015-09-01 08:13:38', 160.64932840316101, 887.834728, 46011),
('CDA1', 'METOP-A', '2015-09-01 09:39:39', 28.668052730278006, '2015-09-01 09:47:22', 55.291433365215099, '2015-09-01 09:55:03', 195.37710788844643, 923.718758, 46012),
('CDA1', 'METOP-A', '2015-09-01 11:20:09', 44.866210724417016, '2015-09-01 11:27:55', 88.147082116548134, '2015-09-01 11:35:39', 227.85283058435132, 930.309171, 46013),
('CDA1', 'METOP-A', '2015-09-01 13:00:21', 64.099186057594977, '2015-09-01 13:08:03', 69.386741200588602, '2015-09-01 13:15:45', 257.29136007789896, 924.125122, 46014),
('CDA1', 'METOP-A', '2015-09-01 14:40:19', 87.279053370514859, '2015-09-01 14:47:59', 64.584814885569486, '2015-09-01 14:55:40', 282.89940850078244, 921.457021, 46015),
('CDA1', 'METOP-A', '2015-09-01 16:20:16', 114.56484795992259, '2015-09-01 16:27:59', 76.535293381903102, '2015-09-01 16:35:43', 304.36532841773203, 926.959947, 46016),
('CDA1', 'METOP-A', '2015-09-01 18:00:32', 145.37818714539699, '2015-09-01 18:08:16', 74.426529690613975, '2015-09-01 18:16:02', 322.19153008972489, 930.019627, 46017),
('CDA1', 'METOP-A', '2015-09-01 19:41:27', 178.88642476544609, '2015-09-01 19:49:02', 43.761765299394646, '2015-09-01 19:56:40', 337.50520175601247, 912.787052, 46018),
('CDA1', 'METOP-A', '2015-09-01 21:23:15', 214.34528904951168, '2015-09-01 21:30:25', 25.097361372111603, '2015-09-01 21:37:37', 351.75015409979113, 862.020291, 46019),
('CDA1', 'METOP-A', '2015-09-01 23:05:57', 250.78036364192184, '2015-09-01 23:12:26', 15.017357861950106, '2015-09-01 23:18:57', 6.7166803309648211, 779.819487, 46020)]

Any idea about this behavior? I also found the same problem over the same time period for METOP-B.

Thanks for your attention and congratulations for your project.

Marco Bosco

mraspaud commented 9 years ago

Hi Marco, Thanks for the but report. I'll have a look at it asap.

Best regards, Martin

mraspaud commented 9 years ago

@Box88 I just ran your script and can't reproduce the bug. I tested with both master and develop branches, but I don't see the line you mention... are you running on linux ?

Box88 commented 9 years ago

I am running on windows 7, python 2.7.10 32 bit, pyorbital 0.3.2 (installed using pip install) I copied and pasted the script again and I still get that line. That's weird...

mraspaud commented 9 years ago

ok, could you try using the code from github instead ? Master branch first, develop branch if master doesn't work ? It's possible that the latest changes in master didn't make it to pypi...

Box88 commented 9 years ago

Using the master branch I still get the bug but using the develop branch it is ok! Next week I will continue playing with pyorbital using the develop branch. Thank you very much!

mraspaud commented 9 years ago

@Box88 Ok, nice to hear. Fyi develop has now been merged to master, and v1.0.0 was just released.