pytroll / pyorbital

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

Calculation method. #4

Closed sgongar closed 9 years ago

sgongar commented 9 years ago

I was testing PyOrbital comparing it with a well know, for me, like PyEphem but the results aren't similar.

The test code that I used was:

import ephem
from math import degrees
import pyorbital.orbital
import datetime

satellite_name = "HUMSAT"
line1 = "1 27848U 03031J   15160.96355163  .00000275  00000-0  14593-3 0  9996"
line2 = "2 27848  98.7047 169.4969 0009055 336.4898  80.2278 14.21451973619325" 
start_time = 1431832673
(lon, lat, ele) = (-8, 42, 0)

def sim_pyephem():  
    observer = ephem.Observer()
    observer.lon = lon
    observer.lat = lat
    observer.elevation = ele
    observer.date = (start_time + 2440587.5*86400)/86400 - 2415020

    satellite = ephem.readtle(satellite_name, line1, line2)
    satellite.compute(observer)

    alt1 = float(repr(satellite.alt))
    alt1 = degrees(alt1)
    az1 = float(repr(satellite.az))
    az1 = degrees(az1)

    print "PyEphem result:  ", start_time, alt1, az1

def sim_pyorbital():
    satellite = pyorbital.orbital.Orbital(satellite_name, line1=line1, line2=line2)
    time1 = datetime.datetime.utcfromtimestamp(start_time)
    az1, alt1 = satellite.get_observer_look(time1, lon, lat, ele)

    print "PyOrbital result:", start_time, alt1, az1    

sim_pyephem()
sim_pyorbital()

Could anybody help me? I don't know if I'm doing something wrong. Thanks!

mraspaud commented 9 years ago

@sgongar could you please provide some of the results also ? how much do pyephem and pyorbital differ ?

sgongar commented 9 years ago

The code I've shown you give the following output.

PyEphem result:   1431832673 -77.1037993892 104.170186164
PyOrbital result: 1431832673 -50.8534286102 163.696181969

I discovered the difference in values while performing a comparison project propagators. Here you have two pictures, one showing the azimuth and a the altitude of four different propagators.

alt az

The x-axis show the time in UNIX format.

You can find the source code for this project in this repo. https://github.com/satnet-project/propagators

mraspaud commented 9 years ago

@sgongar Ok, I looked at this a bit, and I have to start by admitting that I don't know anything about pyephem (except that array computations aren't supported, which is why we started pyorbital).

Now, I took your script and tested it with a few tles and satellites that we receive here, I the pyephem results don't make sense. For example, I attach the image of the reception overview we have here for NOAA-15 (elevation angle in green) and the results of your (adjusted) script: noaa15

PyEphem result:   2015-08-27 13:07:30 -85.3961544756 313.983768825
PyOrbital result: 2015-08-27 13:07:30 16.4159591386 36.6603352198

As you can see, pyephem predicts the satellite would be more or less on the other side of the earth at the time we actually receive it ! (note the time on the picture is local, so it would be 13:00 to 13:20 utc)

The maximum elevation of the pass was about 16.5 degrees, so I can't see any problem with pyorbital. I would be curious to hear if you find an explanation though, so you're very welcome to tell us about your findings.

I tried to change as little as possible to your script but I attach it here again for you to see

import ephem
from math import degrees
import pyorbital.orbital
import datetime

#satellite_name = "HUMSAT"
#line1 = "1 27848U 03031J   15160.96355163  .00000275  00000-0  14593-3 0  9996"
#line2 = "2 27848  98.7047 169.4969 0009055 336.4898  80.2278 14.21451973619325"
#start_time = 1431832673

satellite_name = "NOAA 15"
line1 = "1 25338U 98030A   15238.48036811  .00000066  00000-0  46773-4 0  9994"
line2 = "2 25338  98.7760 237.1745 0009991 330.5088  29.5527 14.25655225898859"
start_time = datetime.datetime(2015, 8, 27, 13, 7, 30)

#(lon, lat, ele) = (-8, 42, 0)
(lon, lat, ele) = (16.14, 58.58, 0)

def sim_pyephem():
    observer = ephem.Observer()
    observer.lon = lon
    observer.lat = lat
    observer.elevation = ele
    #observer.date = (start_time + 2440587.5 * 86400) / 86400 - 2415020
    observer.date = start_time
    satellite = ephem.readtle(satellite_name, line1, line2)
    satellite.compute(observer)

    alt1 = float(repr(satellite.alt))
    alt1 = degrees(alt1)
    az1 = float(repr(satellite.az))
    az1 = degrees(az1)

    print "PyEphem result:  ", start_time, alt1, az1

def sim_pyorbital():
    satellite = pyorbital.orbital.Orbital(satellite_name, line1=line1, line2=line2)
    #time1 = datetime.datetime.utcfromtimestamp(start_time)
    time1 = start_time
    az1, alt1 = satellite.get_observer_look(time1, lon, lat, ele)
    print "PyOrbital result:", start_time, alt1, az1

sim_pyephem()
sim_pyorbital()
mraspaud commented 9 years ago

@sgongar Just found the bug in the pyephem use, we need to write:

    observer.lon = str(lon)
    observer.lat = str(lat)

somehow pyephem doesn't like floats for lon and lat ??? :( The result is now:

PyEphem result:   2015-08-27 13:07:30 16.4698673721 36.6612384455
PyOrbital result: 2015-08-27 13:07:30 16.4159591386 36.6603352198
sgongar commented 9 years ago

Thanks @mraspaud , now it seems the problem is in my code because all propagators gave me the same value as PyEphem.

These days'll check my code but I think this issue in this repository can be considered closed.