I tried tracking a satellite - or rather confirming the data I got from the AMSAT pass prediction web page - using python, more precisly the sgp4 and navpy packages. Now I don't really now what I'm doing here, but I thought this would work:

from sgp4.api import Satrec
from sgp4.api import jday
import navpy
import numpy as np


s = '1 25544U 98067A   21336.70775771  .00004247  00000-0  85987-4 0  9994'
t = '2 25544  51.6437 229.4568 0004362 271.6918 182.8759 15.48732282314711'
satellite = Satrec.twoline2rv(s, t)

jd, fr = jday(2021, 12, 8, 9, 4, 0)
e, r, v = satellite.sgp4(jd, fr) # e = error, r = position vector, v = speed vector

lat, lon, alt = 48.1, 16.3, 250 
r0 = navpy.lla2ecef(lat, lon, alt, latlon_unit='deg', alt_unit='m', model='wgs84')

x, y, z = r - r0
az = np.arctan2(x,y)
el = np.arctan2(z,y)
az_deg = 180*az/np.pi
el_deg = 180*el/np.pi

What I expected from this code is, that it would give me a simillar output to what the AMSAT page those, i.e., give me elevation and azimuth of the satellite, relative to my position (r0) at a certain time. However, the results don't match at all and I suspect it has something to do with the used coordinate systems, but I can't point my finger at it. Anyone can help me with this?


It's the coordinate system. You need an additional tool to make the conversion. It's right there on the home page (https://pypi.org/project/sgp4/):

Note that the SGP4 propagator returns raw x,y,z Cartesian coordinates in a “True Equator Mean Equinox” (TEME) reference frame that’s centered on the Earth but does not rotate with it — an “Earth centered inertial” (ECI) reference frame. The SGP4 propagator itself does not implement the math to convert these positions into more official ECI frames like J2000 or the ICRF; nor to convert positions into any Earth-centered Earth-fixed (ECEF) frames like the ITRS; nor to convert them to latitudes and longitudes through an Earth ellipsoid like WGS84.

For conversions into other coordinate frames, look for a comprehensive astronomy library that is built atop this one, like the Skyfield library

Note carefully that TEME is an ECI frame, not the only ECI frame. There are many different options for ECI, each of which tries to be "inertial" in a slightly different way, though none of them succeed exactly. Similarly, WGS84 is an ECEF frame, not the only ECEF frame. In fact, SGP4 is so old, that it actually uses WGS72 internally.

  • $\begingroup$ It is not obvious to me how to use Skyfield to make my code work, can you explain further? (As I commented below, M2FKXY's code works for me, but I'd be interested in wrapping my head around this. $\endgroup$
    – RadaD
    Dec 14 '21 at 12:19

What are the results when using skyfield? Here is some sample code to print an azimuth and elevation at a certain time, and avoid any reference frame and coordinates manipulations.

from skyfield.api import load, wgs84

ts = load.timescale()

satUrl = 'http://celestrak.com/NORAD/elements/active.txt'
sats = load.tle_file(satUrl)
print('Loaded', len(sats), 'satellites')

>> Loaded 4848 satellites

satNames = {sat.name: sat for sat in sats}
sat = satNames["SENTINEL-1A"]

>> SENTINEL-1A catalog #39634 epoch 2021-12-08 08:16:29 UTC

myPos = wgs84.latlon(+60, +25)

t0 = ts.utc(2021, 12, 10, 0)
t1 = ts.utc(2021, 12, 10, 12)

t, events = sat.find_events(myPos, t0, t1, altitude_degrees = 20)

for ti, event in zip(t, events):
    name = ("rise above 20°", "culminate", "set below 20°")[event]
    print(ti.utc_strftime('%Y %b %d %H:%M:%S'), name)
>> 2021 Dec 10 05:11:52 rise above 20°
>> 2021 Dec 10 05:15:11 culminate
>> 2021 Dec 10 05:18:30 set below 20°
>> 2021 Dec 10 06:50:38 rise above 20°
>> 2021 Dec 10 06:52:21 culminate
>> 2021 Dec 10 06:54:05 set below 20°

t = ts.utc(2021, 12, 10, 5, 15, 11)

diff = sat - myPos
topocentric = diff.at(t)
alt, az, dist = topocentric.altaz()

print("Altitude ", alt, "; azimuth ", az, " at ", t.utc_strftime(), )      
>> Altitude  85deg 49' 06.4" ; azimuth  105deg 50' 27.2"  at  2021-12-10 05:15:11 UTC
  • $\begingroup$ Thanks, I get the same AOS and LOS data as with the mentioned pass prediction tool. Is it correct that I would have to correct for a observer altitude (my altitude) myself? The code does not take this into account, right? $\endgroup$
    – RadaD
    Dec 14 '21 at 12:18
  • 1
    $\begingroup$ @RadaD, to account for your altitude construct myPos = wgs84.latlon(lat, lon, alt) where alt (altitude relative to ellipsoid) is an extra parameters (units in meters). It's not obvious at all in the Skyfield documentation that you can do this. $\endgroup$ Dec 17 '21 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.