Loading OMM and TLE satellite data

How are satellite orbits disseminated?

(Image by @fisadev)

• The only public source of orbital data for non-amateur Earth artificial satellites is the 18th Space Control Squadron (18 SPCS) from the United States Space Force (https://www.space-track.org/ and https://celestrak.org/NORAD/).

• The data is disseminated as general perturbations (GP) orbital data following the Simplified General Perturbations 4 (SGP4) propagation model.

• Such GP orbital data has traditionally been transmitted in Two-Line Element Sets or TLEs, although in May 2020 it started being published in CCSDS Orbit Mean-Elements Message (OMM) format as well.

• The advantage of SGP4 is that it’s an analytical model that includes the most important natural perturbations, and it’s much faster than an equivalent numerical model. Therefore, it’s good for initial orbit design based on low-precision requirements.

However, it turns out that GP data in general, and TLEs in particular, are poorly understood even by professionals ([1], [2], [3]). The core issue is that TLEs and OMMs contain Brouwer mean elements, which cannot be directly translated to osculating elements.

From “Spacetrack Report #3”:

The most important point to be noted is that not just any prediction model will suffice. The NORAD element sets are “mean” values obtained by removing periodic variations in a particular way. In order to obtain good predictions, these periodic variations must be reconstructed (by the prediction model) in exactly the same way they were removed by NORAD.

From the Celestrak TLE FAQ:

The elements in the two-line element sets are mean elements calculated to fit a set of observations using a specific model —the SGP4/SDP4 orbital model. Just as you shouldn’t expect the arithmetic and geometric means of a set of data to have the same value, you shouldn’t expect mean elements from different element sets —calculated using different orbital models— to have the same value. The short answer is that you cannot simply reformat the data unless you are willing to accept predictions with unpredictable errors.

From “Revisiting Spacetrack Report #3: Rev 3”:

Simply converting the orbital elements to an osculating state vector and propagating with a numerical propagator is equally invalid.

Therefore, the correct way of using GP data is:

1. Propagate the TLE or OMM using the SGP4 algorithm, which produces cartesian elements (position and velocity) in the True Equator Mean Equinox (TEME) reference frame.

2. Convert the resulting coordinates from TEME to the desired one (GCRS/ECI, ITRS/ECEF, other) (see note below).

“Revisiting Spacetrack Report #3: Rev 3” warns that there is an inherent ambiguity in the reference frame used by the 18 SPCS to generate their GP data (TEME of Date vs TEME of Epoch) and there is no confirmation from official sources of which one is used. In any case, if you need sub-kilometer precision, GP data with SGP4 is not enough for you anyway.

Loading general perturbations data

As explained in the Orbit Mean-Elements Messages (OMMs) support assessment deliverable of OpenSatCom, OMM input/output support in open source libraries is somewhat scattered. Luckily, python-sgp4 supports reading OMM in CSV and XML format, as well as usual TLE and 3LE formats. On the other hand, Astropy has accurate transformations from TEME to other reference frames.

# From https://github.com/poliastro/poliastro/blob/main/contrib/satgpio.py
"""
Author: Juan Luis Cano Rodríguez

Code to read GP data from Celestrak using the HTTP API and python-sgp4.

Requires some extra dependencies:

  $ pip install httpx sgp4

This is similar to https://gitlab.com/librespacefoundation/python-satellitetle,
but uses the XML API instead and returns a `Satrec` object from sgp4 directly.

"""

import io
import json
import xml.etree.ElementTree as ET

import httpx
from sgp4 import exporter, omm
from sgp4.api import Satrec


def _generate_url(catalog_number, international_designator, name):
    params = {
        "CATNR": catalog_number,
        "INTDES": international_designator,
        "NAME": name,
    }
    param_names = [
        param_name
        for param_name, param_value in params.items()
        if param_value is not None
    ]
    if len(param_names) != 1:
        raise ValueError(
            "Specify exactly one of catalog_number, international_designator, or name"
        )
    param_name = param_names[0]
    param_value = params[param_name]
    url = (
        "https://celestrak.org/NORAD/elements/gp.php?"
        f"{param_name}={param_value}"
        "&FORMAT=XML"
    )
    return url


def _segments_from_query(url):
    response = httpx.get(url)
    response.raise_for_status()

    if response.text == "No GP data found":
        raise ValueError(
            f"Query '{url}' did not return any results, try a different one"
        )
    tree = ET.parse(io.StringIO(response.text))
    root = tree.getroot()

    yield from omm.parse_xml(io.StringIO(response.text))


def load_gp_from_celestrak(
    *, catalog_number=None, international_designator=None, name=None
):
    """Load general perturbations orbital data from Celestrak.

    Returns
    -------
    Satrec
        Orbital data from specified object.

    Notes
    -----
    This uses the OMM XML format from Celestrak as described in [1]_.

    References
    ----------
    .. [1] Kelso, T.S. "A New Way to Obtain GP Data (aka TLEs)"
       https://celestrak.org/NORAD/documentation/gp-data-formats.php

    """
    # Assemble query, raise an error if malformed
    url = _generate_url(catalog_number, international_designator, name)

    # Make API call, raise an error if data is malformed
    for segment in _segments_from_query(url):
        # Initialize and return Satrec object
        sat = Satrec()
        omm.initialize(sat, segment)

        yield sat


def print_sat(sat, name):
    """Prints Satrec object in convenient form."""
    print(json.dumps(exporter.export_omm(sat, name), indent=2))

sat = list(load_gp_from_celestrak(name="ISS (Zarya)"))[0]
print_sat(sat, "ISS (Zarya)")

{
  "OBJECT_NAME": "ISS (Zarya)",
  "OBJECT_ID": "1998-067A",
  "CENTER_NAME": "EARTH",
  "REF_FRAME": "TEME",
  "TIME_SYSTEM": "UTC",
  "MEAN_ELEMENT_THEORY": "SGP4",
  "EPOCH": "2022-07-10T03:46:49.156896",
  "MEAN_MOTION": 15.498853529999998,
  "ECCENTRICITY": 0.0004765,
  "INCLINATION": 51.6428,
  "RA_OF_ASC_NODE": 222.824,
  "ARG_OF_PERICENTER": 353.1244,
  "MEAN_ANOMALY": 111.7704,
  "EPHEMERIS_TYPE": 0,
  "CLASSIFICATION_TYPE": "U",
  "NORAD_CAT_ID": 25544,
  "ELEMENT_SET_NO": 999,
  "REV_AT_EPOCH": 34874,
  "BSTAR": 0.00010529,
  "MEAN_MOTION_DOT": 5.534e-05,
  "MEAN_MOTION_DDOT": 0.0
}

The generator load_gp_from_celestrak might yield more than one satellite:

len(list(load_gp_from_celestrak(name="COSMOS 1408 DEB")))

692

Creating ephemerides from general perturbations data

Now that we know how to load GP data, let’s propagate it using the SGP4 algorithm:

from astropy import units as u
from astropy.time import Time

now = Time.now()
now.jd1, now.jd2

(2459771.0, 0.16455130069444446)

error, r, v = sat.sgp4(now.jd1, now.jd2)
assert error == 0

r << u.km

$[-651.64965,~-5125.4018,~4402.7695] \; \mathrm{km}$

v << (u.km / u.s)

$[6.5644841,~2.068562,~3.3756512] \; \mathrm{\frac{km}{s}}$

Note that Satrec also supports propagating for arrays of time values using the .sgp4_array method:

import numpy as np

from astropy.coordinates import CartesianRepresentation, CartesianDifferential

from poliastro.util import time_range

times = time_range(now, end=now + (1 << u.h), periods=3)

errors, rs, vs = sat.sgp4_array(times.jd1, times.jd2)
assert (errors == 0).all()

CartesianRepresentation(rs << u.km, xyz_axis=-1)

<CartesianRepresentation (x, y, z) in km
    [( -651.64964926, -5125.40181287,  4402.76952168),
     ( 5511.89886557,  3914.0454116 ,   718.6678296 ),
     (-4220.22649891,  1682.15539479, -5057.51706039)]>

CartesianDifferential(vs << (u.km / u.s), xyz_axis=-1)

<CartesianDifferential (d_x, d_y, d_z) in km / s
    [( 6.56448405,  2.068562  ,  3.37565123),
     (-2.24423769,  4.26596893, -5.95180601),
     (-4.5875843 , -5.82447553,  1.89634597)]>

Mixing all together, and leveraging poliastro Ephem objects, we can compute ephemerides from GP orbital data:

from warnings import warn

from astropy.coordinates import TEME, GCRS

from poliastro.ephem import Ephem
from poliastro.frames import Planes


def ephem_from_gp(sat, times):
    errors, rs, vs = sat.sgp4_array(times.jd1, times.jd2)
    if not (errors == 0).all():
        warn(
            "Some objects could not be propagated, "
            "proceeding with the rest",
            stacklevel=2,
        )
        rs = rs[errors == 0]
        vs = vs[errors == 0]
        times = times[errors == 0]

    cart_teme = CartesianRepresentation(
        rs << u.km,
        xyz_axis=-1,
        differentials=CartesianDifferential(
            vs << (u.km / u.s),
            xyz_axis=-1,
        ),
    )
    cart_gcrs = (
        TEME(cart_teme, obstime=times)
        .transform_to(GCRS(obstime=times))
        .cartesian
    )

    return Ephem(cart_gcrs, times, plane=Planes.EARTH_EQUATOR)

iss_ephem = ephem_from_gp(sat, time_range(now, end=now + (3 << u.h)))
iss_ephem

Ephemerides at 50 epochs from 2022-07-10 15:56:57.232380 (UTC) to 2022-07-10 18:56:57.232380 (UTC)

And plot it!

from poliastro.bodies import Earth
from poliastro.plotting import OrbitPlotter3D

plotter = OrbitPlotter3D()
plotter.set_attractor(Earth)
plotter.plot_ephem(iss_ephem, color="#333", label="ISS", trail=True)

plotter.show()

We can also plot the trajectory of the ISS for 90 minutes, plus the trajectories of the first 25 debris fragments:

from itertools import islice

epochs = time_range(now, end=now + (90 << u.minute))

iss_ephem = ephem_from_gp(sat, epochs)

plotter = OrbitPlotter3D()
plotter.set_attractor(Earth)
plotter.plot_ephem(iss_ephem, color="#333", label="ISS", trail=True)

for debris_fragment in islice(
    load_gp_from_celestrak(name="COSMOS 1408 DEB"), 25
):
    debris_ephem = ephem_from_gp(debris_fragment, epochs)
    plotter.plot_ephem(
        debris_ephem, color="#666", label=debris_fragment.satnum, trail=True
    )

plotter.show()

/tmp/ipykernel_2218/3231179520.py:14: UserWarning:

Some objects could not be propagated, proceeding with the rest