Working with Earth Satellites Using Astropy Coordinates

This document discusses Two-Line Element ephemerides and the True Equator, Mean Equinox frame. For satellite ephemerides given directly in geocentric ITRS coordinates (e.g. ILRS ephemeris format) please see the example transform to AltAz below starting with the geocentric ITRS coordinate frame.

Satellite data is normally provided in the Two-Line Element (TLE) format (see here for a definition). These datasets are designed to be used in combination with a theory for orbital propagation model to predict the positions of satellites.

The history of such models is discussed in detail in Vallado et al (2006) who also provide a reference implementation of the SGP4 orbital propagation code, designed to be compatible with the TLE sets provided by the United States Department of Defense, which are available from a source like Celestrak.

The output coordinate frame of the SGP4 model is the True Equator, Mean Equinox frame (TEME), which is one of the frames built-in to astropy.coordinates. TEME is an Earth-centered inertial frame (i.e., it does not rotate with respect to the stars). Several definitions exist; astropy uses the implementation described in Vallado et al (2006).

Finding TEME Coordinates from TLE Data

There is currently no support in astropy.coordinates for computing satellite orbits from TLE orbital element sets. Full support for handling TLE files is available in the Skyfield library, but some advice for dealing with satellite data in astropy is below.

You will need some external library to compute the position and velocity of the satellite from the TLE orbital elements. The SGP4 library can do this. An example of using this library to find the TEME coordinates of a satellite is:

>>> from sgp4.api import Satrec
>>> from sgp4.api import SGP4_ERRORS
>>> s = '1 25544U 98067A   19343.69339541  .00001764  00000-0  38792-4 0  9991'
>>> t = '2 25544  51.6439 211.2001 0007417  17.6667  85.6398 15.50103472202482'
>>> satellite = Satrec.twoline2rv(s, t)

The satellite object has a method, satellite.sgp4, that will try to compute the TEME position and velocity at a given time:

>>> from astropy.time import Time
>>> t = Time(2458827.362605, format='jd')
>>> error_code, teme_p, teme_v = satellite.sgp4(t.jd1, t.jd2)  # in km and km/s
>>> if error_code != 0:
...     raise RuntimeError(SGP4_ERRORS[error_code])

Now that we have the position and velocity in kilometers and kilometers per second, we can create a position in the TEME reference frame:

>>> from astropy.coordinates import TEME, CartesianDifferential, CartesianRepresentation
>>> from astropy import units as u
>>> teme_p = CartesianRepresentation(teme_p*u.km)
>>> teme_v = CartesianDifferential(teme_v*u.km/u.s)
>>> teme = TEME(teme_p.with_differentials(teme_v), obstime=t)

Note how we are careful to set the observed time of the TEME frame to the time at which we calculated satellite position.

Transforming TEME to Other Coordinate Systems

Once you have satellite positions in TEME coordinates they can be transformed into any astropy.coordinates frame.

For example, to find the overhead latitude, longitude, and height of the satellite:

>>> from astropy.coordinates import ITRS
>>> itrs_geo = teme.transform_to(ITRS(obstime=t))
>>> location = itrs_geo.earth_location
>>> location.geodetic  
GeodeticLocation(lon=<Longitude 160.34199789 deg>, lat=<Latitude -24.6609379 deg>, height=<Quantity 420.17927591 km>)

Or, if you want to find the altitude and azimuth of the satellite from a particular location:

Note

In this example, the intermediate step of manually setting up a topocentric ITRS frame is necessary in order to avoid the change in stellar aberration that would occur if a direct transform from geocentric to topocentric coordinates using transform_to was used. Please see the documentation of the ITRS frame for more details.

>>> from astropy.coordinates import EarthLocation, AltAz
>>> siding_spring = EarthLocation.of_site('aao')  
>>> topo_itrs_repr = itrs_geo.cartesian.without_differentials() - siding_spring.get_itrs(t).cartesian
>>> itrs_topo = ITRS(topo_itrs_repr, obstime = t, location=siding_spring)
>>> aa = itrs_topo.transform_to(AltAz(obstime=t, location=siding_spring))
>>> aa.alt  
<Latitude 10.94799670 deg>
>>> aa.az  
<Longitude 59.28803392 deg>

For a stationary observer, velocity in the ITRS is independent of location, so if you want to carry the velocity to the topocentric frame, you can do so as follows:

>>> itrs_geo_p = itrs_geo.cartesian.without_differentials()
>>> itrs_geo_v = itrs_geo.cartesian.differentials['s']
>>> topo_itrs_p = itrs_geo_p - siding_spring.get_itrs(t).cartesian
>>> topo_itrs_repr = topo_itrs_p.with_differentials(itrs_geo_v)
>>> itrs_topo = ITRS(topo_itrs_repr, obstime = t, location=siding_spring)