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)