Comparación de resultados con Pyephem

Estoy tratando de calcular la posición de cada satélite de la constelación de GPS, usando archivos TLE y biblioteca pyephem , quiero saber la posición Elevación-Acimut. Pero cuando comparo los resultados con la página calsky con los resultados de mi código son diferentes, ¿por qué ?

Los datos iniciales son: Bogotá, Colombia
Longitud = -74.0817 Latitud = 4.60971 Alt : 2578.0
Fecha: 2017-04-18 Hora: 10:08:51.61

Resultados en CALSKY >45º

PRN 27/GPS BIIF-(39166 2013-023-A) az: 152,7° SSE h: 61,1°
PRN 16/GPS BIIR-(27663 2003-005-A) az: 341,8° NNW h: 52,5°
PRN 14/GPS BIIR-(26605 2000-071-A) az: 109,5° ESE h: 49,5°

Resultados en mi código > 45º

GPS BIIR-7 (PRN 18) - azim: 27,68 elev:
58,34 GPS BIIR-9 (PRN 21) - azim: 200,54 elev: 55,05
GPS BIIF-1 (PRN 25) - azim: 111,35 elev: 46,92

def GetPositionELAZ(Time,Longitude,Latitude,constelation,AllTheLines):
'''
- time: Time for look the satellite
- longitude: observer longitude
- latitude: observer latitude
- constelation: 
            1 - just GPS 
            2 - GPS - GLONASS
            3 - GPS - GLONASS - GALILEO
            
RETURN: (Tuble of lists...)
Satellite Name, Satellite Elevation, Satellite Azimuth, Type Satellite (1:GPS, 2: GLO, 3:GAL)
'''

contain_GPS = 'GPS'
contain_GLONASS = 'COSMOS'
contain_GALILEO = 'GSAT'

sat_alt, sat_az, sat_name, sat_type = [], [], [], []
original_name = []

longlat = str(Longitude) + "," + str(Latitude)
g = geocoder.elevation(longlat) #lon,lat
#print g.meters
import ephem
observer = ephem.Observer()
observer.long = Longitude
observer.lat = Latitude
observer.date = Time
observer.elevation = 2600. #float(str(g.meters))

GPS_list = 'https://celestrak.org/NORAD/elements/gps-ops.txt'
GLONASS_list = 'https://celestrak.org/NORAD/elements/glo-ops.txt'
GALILEO_list = 'https://celestrak.org/NORAD/elements/galileo.txt'

#'http://www.amsat.org/amsat/ftp/keps/current/nasabare.txt').readlines()
GPS2_list = 'http://www.tle.info/data/gps-ops.txt'
GLONASS2_list = 'http://www.tle.info/data/glo-ops.txt'

if constelation == 1:
    list_TLEs = [GPS_list]
if constelation == 2:
    list_TLEs = [GPS_list,GLONASS_list]
if constelation == 3:
    list_TLEs = [GPS_list,GLONASS_list,GALILEO_list]

for satlist in list_TLEs:
    tles = AllTheLines #urllib2.urlopen(satlist).readlines() 
    tles = [item.strip() for item in tles]
    tles = [(tles[i],tles[i+1],tles[i+2]) for i in xrange(0,len(tles)-2,3)]
    
    for tle in tles:

        try:
            sat = ephem.readtle(tle[0], tle[1], tle[2])
            rt, ra, tt, ta, st, sa = observer.next_pass(sat)

            if rt is not None and st is not None:
                sat.compute(observer)

                if Time >= ephem.localtime(st) and Time <= ephem.localtime(rt):
                    Title = tle[0]
                    sat_alt.append(np.rad2deg(sat.alt))
                    sat_az.append(np.rad2deg(sat.az))
                    original_name.append(tle[0])

                    text2 = Title.rsplit(')', 1)[0]
                    NamePRN = text2.rsplit('(', 1)[1]
                    sat_name.append(NamePRN)

                    if contain_GPS in Title:
                        sat_type.append(1)
                    elif contain_GLONASS in Title:
                        sat_type.append(2)
                    elif contain_GALILEO in Title:
                        sat_type.append(3)

        except ValueError as e:
            print e
return sat_name, sat_alt, sat_az, sat_type, original_name

Llamando a la función

GPS_list = 'https://celestrak.org/NORAD/elements/gps-ops.txt'
Downloaded = urllib2.urlopen(GPS_list).readlines()
TimeNow = datetime.datetime.now()
print Time2SecondsOfDay(TimeNow), 'time: ',TimeNow

names,elevs,azs,types,originalName = 
GetPositionELAZ(TimeNow,Longitude,Latitude,1,Downloaded)
Espero ver los mismos vehículos en el mismo tiempo, pero solo veo otros vehículos, pero no sé por qué no es lo mismo. para comparar algunos vehiculos solo muestro los satelites sobre 45º.
Esto probablemente no sea útil, pero tienes tu latitud y longitud invertidas. Déjame saber si eso ayuda/no ayuda.
gracias @barrycarter Tenía la latitud/longitud incorrecta, actualicé los resultados (y la pregunta), aunque todavía encuentro algunas diferencias en los resultados.
Si no me equivoco, la información de TLE te dará la órbita. A partir de eso, obtienes calcular la posición y la velocidad de la nave espacial en el marco ECI. Sin embargo, debe convertir eso al marco ECEF para calcular la información de altitud y elevación, y no veo que eso suceda en su código. Dicho esto, si hay una biblioteca que funciona bien, ¿por qué necesita volver a codificar esa función y no solo usar la biblioteca que funciona?
Además, ¿cuál es el código de sat.compute()? ¿Podría ocurrir allí la conversión del marco de referencia?
gracias @ChrisR, sobre el marco ECEF y ECI, cuando busca la elevación y el azimut, no importa el marco, pero necesita la elevación de la ubicación (uso el geocodificador de la biblioteca para esto) porque la elevación y el azimut son respecto a el horizonte
@mikesneider No, los marcos (casi) siempre importan. No se puede restar el vector del satélite expresado en el marco ECI con la posición de la estación que comúnmente se expresa en ECEF.
@ChrisR la respuesta literal es: "Para ser precisos, el marco de referencia de las coordenadas inerciales centradas en la Tierra (ECI) producidas por el modelo orbital SGP4/SDP4 es el ecuador verdadero, equinoccio medio (TEME) de época". fuente

Respuestas (1)

Descubrí que había algo mal en el código: la función observer.next_pass(sat), devuelve un vector de 5 valores

0 Tiempo
de subida 1 Acimut de subida
2 Tiempo de altitud
máxima 3 Altitud máxima
4 Tiempo
establecido 5 Acimut establecido

Entonces, la instrucción if Time >= ephem.localtime(st) and Time <= ephem.localtime(rt):que estaba usando era una comparación innecesaria porque ya calculaste la hora con la instrucción observer.date = Timey también la estaba comparando con la hora local, la mejor opción es tener todo en hora UTC.
y un consejo con satelites no uses datetime.now(), usa datetime.utcnow(), y luego puedes agregar el GTM de tu region.