[SOLVED] Python satellite tracking with spg4, pyephem

I'm trying to write a basic python scrip that will track a given satellite, defined with tle's, from a given location. I'm not a asto/orbital person but am trying to become smarter on it.

I am running into a problem where the different models I'm using are giving me very different position answers. I have tried using: pyEphem spg4 predict (exec system call from script)

The satellites I'm testing with are the ISS and directv10 (one fixed, one moving- with internet tracking available for verification):

0 Direct10
1 31862U 07032A   13099.15996183 -.00000126  00000-0  10000-3 0  1194
2 31862 000.0489 046.9646 0000388 001.7833 103.5813 01.00271667 21104
0 ISS
1 25544U 98067A   13112.50724749  .00016717  00000-0  10270-3 0  9148
2 25544  51.6465  24.5919 0009906 171.1474 188.9854 15.52429950 26067

I have modified the predict source to give me the eci location so I can use that to know the real location. I also have it give the az, el, range for using to verify observations. I'm using spg4 for getting the real location. For the observed location, I'm using PyEphem.

I'm getting the ECEF position from spg4 with:

def get_real(epoch, sv):
    satellite = twoline2rv(sv.tle1, sv.tle2, wgs84)

    #epoch = time.time()
    obsTime = datetime.datetime.utcfromtimestamp(epoch)
    position, velocity = satellite.propagate(obsTime.year,
                                             obsTime.month,
                                             obsTime.day,
                                             obsTime.hour,
                                             obsTime.minute,
                                             obsTime.second)


    x = position[0]
    y = position[1]
    z = position[2]

    x *= 1000
    y *= 1000
    z *= 1000

My code for pyephem based observations is:

def get_ob(epoch, sv, obsLoc):
    site = ephem.Observer()
    site.lon = str(obsLoc.lat)   # +E -104.77 here
    site.lat = str(obsLoc.lon)   # +N 38.95   here
    site.elevation = obsLoc.alt # meters    0 here
    #epoch = time.time()
    site.date = datetime.datetime.utcfromtimestamp(epoch)

    sat = ephem.readtle(sv.name,sv.tle1,sv.tle2)
    sat.compute(site)

    az       = degrees(sat.az)
    el       = degrees(sat.alt)
    #range in m
    range    = sat.range
    sat_lat  = degrees(sat.sublat)
    sat_long = degrees(sat.sublong)
    # elevation of sat in m
    sat_elev = sat.elevation

    #TODO: switch to using az,el,range for observed location calculation
    #x, y, z    = aer2ecef(az,el,range,38.95,-104.77,80 / 1000)
    x,y,z  = llh2ecef(sat_lat, sat_long, sat_elev)

The llh2ecef conversion:

def llh2ecef (flati,floni, altkmi ):
    #         lat,lon,height to xyz vector
    #
    # input:
    # flat      geodetic latitude in deg
    # flon      longitude in deg
    # altkm     altitude in km
    # output:
    # returns vector x 3 long ECEF in km

    dtr =  pi/180.0;

    flat  = float(flati);
    flon  = float(floni);
    altkm = float(altkmi);

    clat = cos(dtr*flat);
    slat = sin(dtr*flat);
    clon = cos(dtr*flon);
    slon = sin(dtr*flon);

    rrnrm  = radcur (flat);
    rn     = rrnrm[1];
    re     = rrnrm[0];

    ecc1    = ecc;
    esq1    = ecc1*ecc1

    x      =  (rn + altkm) * clat * clon;
    y      =  (rn + altkm) * clat * slon;
    z      =  ( (1-esq1)*rn + altkm ) * slat;

    return x,y,z

aer2ecef:

def aer2ecef(azimuthDeg, elevationDeg, slantRange, obs_lat, obs_long, obs_alt):

    #site ecef in meters
    sitex, sitey, sitez = llh2ecef(obs_lat,obs_long,obs_alt)

    #some needed calculations
    slat = sin(radians(obs_lat))
    slon = sin(radians(obs_long))
    clat = cos(radians(obs_lat))
    clon = cos(radians(obs_long))

    azRad = radians(azimuthDeg)
    elRad = radians(elevationDeg)

    # az,el,range to sez convertion
    south  = -slantRange * cos(elRad) * cos(azRad)
    east   =  slantRange * cos(elRad) * sin(azRad)
    zenith =  slantRange * sin(elRad)


    x = ( slat * clon * south) + (-slon * east) + (clat * clon * zenith) + sitex
    y = ( slat * slon * south) + ( clon * east) + (clat * slon * zenith) + sitey
    z = (-clat *        south) + ( slat * zenith) + sitez

    return x, y, z

When I compare and plot the locations on a 3D globe (using ecef positions), I'm getting answers all over the place. The predict eci position (converted to ecef) matchs what I see on ISS Tracking websites (http://www.n2yo.com/?s=25544)

The result from get_real() is way off in scale and location. The result from get_ob() are right in scale, but wrong in location on globe

example results:

predict based:

sv: ISS predict observed response    @ epoch: 1365630559.000000 : [111.485527, -69.072949, 12351.471383]
sv: ISS predict aer2ecef position(m) @ epoch: 1365630559.000000 : [4731598.706291642, 1844098.7384999825, -4521102.9225004213]
sv: ISS predict ecef position(m) @ epoch: 1365630559.000000 : [-3207559.6840419229, -3937040.5048992992, -4521102.9110000003]
sv: ISS predict ecef2llh(m)      @ epoch: 1365630559.000000 : [-41.67839724680753, -129.170165912171, 6792829.6884068651]
sv: Direct10 predict observed response    @ epoch: 1365630559.000000 : [39.692138, -49.219935, 46791.914833]
sv: Direct10 predict aer2ecef position(m) @ epoch: 1365630559.000000 : [28401835.38849232, 31161334.784188181, 3419.5400331273049]
sv: Direct10 predict ecef position(m) @ epoch: 1365630559.000000 : [-9348629.6463202238, -41113211.570621684, 3419.8620000000005]
sv: Direct10 predict ecef2llh(m)      @ epoch: 1365630559.000000 : [0.0046473273713214715, -102.81051792373036, 42156319.281573996]

python based:

sv: ISS ephem observed response    @ epoch: 1365630559.000000 : [344.067992722211, -72.38297754053431, 12587123.0][degrees(sat.az), degrees(sat.alt), sat.range]
sv: ISS ephem llh location(m)      @ epoch: 1365630559.000000 : [-41.678271938092195, -129.16682754513502, 421062.90625][degrees(sat.sublat0, degrees(sat.sublong), sat.elevation]
sv: ISS ephem xyz location(m)      @ epoch: 1365630559.000000 :[-201637.5647039332, -247524.53652043006, -284203.56557438202][llh2ecef(lat,long,elev)]
sv: ISS spg84 ecef position(m) @ epoch: 1365630559.000000 : [4031874.0758277094, 3087193.8810081254, -4521293.538866323]
sv: ISS spg84 ecef2llh(m)      @ epoch: 1365630559.000000 : [-41.68067424524357, 37.4411722245808, 6792812.8704163525]
sv: Direct10 ephem observed response    @ epoch: 1365630559.000000 : [320.8276456938389, -19.703680198781303, 43887572.0][degrees(sat.az), degrees(sat.alt), sat.range]
sv: Direct10 ephem llh location(m)      @ epoch: 1365630559.000000 : [0.004647324660923812, -102.8070784813048, 35784688.0][degrees(sat.sublat0, degrees(sat.sublong), sat.elevation]
sv: Direct10 ephem xyz location(m)      @ epoch: 1365630559.000000 :[-7933768.6901137345, -34900655.02490133, 2903.0498773286708][llh2ecef(lat,long,elev)]
sv: Direct10 spg84 ecef position(m) @ epoch: 1365630559.000000 : [18612307.532456037, 37832170.97306267, -14060.29781505302]
sv: Direct10 spg84 ecef2llh(m)      @ epoch: 1365630559.000000 : [-0.019106864351793953, 63.80418030988552, 42156299.077687643]

The az,el and range don't match between the two observations. The positions don't match for the "true" locations. (The lat and long does but height doesn't after a ecef2llh conversion.

When compared with the web based tracker, I notice the predict "true" llh locations match the website. For directv10, pyEphem matchs the azimuth and elevation- but not for the ISS

When I plot them on globe, the predict eci "true" location is in the right spot - matches tracker website). The spg84 ecef position (which I thought should be the same as predict, is on the other side of globe. The predict "observed" location, is close to the spg84 location. The pyEphem is completely off in altitude and not displayed (way too low, inside earth).

So my question is where am I using the python models wrong? My understanding was that spg84 propagate() call, should return the exec position of the satellite in meters. I would have thought at should match the predict position after the eci2efec conversion. I also would have expected that the match then the llh2ecef() when using sat.sublat,sat.sublong,sat.elevation.

As I said, I'm new to all things orbiting so I'm sure I an making a simple math error or soemthing. I have tried to google and search for answers, examples and tutorials as much as possible but nothing has helped so far (I have tried multiple ecef2llh and llh2ecef methods to try to work out those bugs.

Any suggestions, advice, pointers in the right direction would be greatly appreciated. I can post/send the complete cod eI'm using if that would be helpful for someone. I tried to make sure I posted the important parts here and didn't want to make this (already very) long post and longer.

Thanks for the help.

Aaron

UPDATE:

I found atleast part of the issue. spg84.propagate() returns the location in ECI, not ECEF. Quick run through eci2ecef and it lines up perfectly with the predict response.

I always seem to find solutions after posting for help ;)

Now need to figure out what is going on with the observer locations. This boils down to: How can I take the result from pyEphem.compute() and get the ecef position for the satellite? Prefer to do it with az,el, range values, not latitude, longitude, elevation.

I'm guessing the bug in my aer2ecef call.

Thanks.

UPDATE 2:

Got the observation to line up with the "true" position. Looks like I had a units issue. The working code:

az       = degrees(sat.az)
el       = degrees(sat.alt)
#range in km
range    = sat.range
sat_lat  = degrees(sat.sublat)
sat_long = degrees(sat.sublong)
# elevation of sat in km
sat_elev = sat.elevation


#x, y, z    = aer2ecef(az,el,range,obsLoc.lat,obsLoc.long,obsLoc.alt)
x,y,z  = llh2ecef(sat_lat, sat_long, sat_elev / 1000)

x *= 1000
y *= 1000
z *= 1000
return x,y,z

Now just need aer2ecef() method to return proper position...

Python satellite tracking with spg4, pyephem - positions not matching