Code for plotting orbital parameters

tle_swap.py

#!/usr/bin/env python3
from __future__ import print_function, division
from spacetrack import SpaceTrackClient
import os
import json
import ephem
from sgp4.earth_gravity import wgs84
from sgp4.io import twoline2rv
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from astropy.time import Time

class twolineelement:
    """TLE class"""

    def __init__(self, tle0, tle1, tle2):
        """Define a TLE"""

        self.tle0 = tle0
        self.tle1 = tle1
        self.tle2 = tle2
        if tle0[:2] == "0 ":
            self.name = tle0[2:].rstrip()
        else:
            self.name = tle0.rstrip()
        if tle1.split(" ")[1] == "":
            self.id = int(tle1.split(" ")[2][:4])
        else:
            self.id = int(tle1.split(" ")[1][:5])
        self.desig = tle1.split(" ")[2]

def download_tles_from_spacetrack(settings, intldes_query, fname):
    # Login
    st = SpaceTrackClient(settings["username"], settings["password"])

    # Find objects
    satcat = st.satcat(intldes=intldes_query)

    # Get NORAD IDs
    norad_cat_ids = sorted([int(s["NORAD_CAT_ID"]) for s in satcat])

    # Download TLEs
    with open(fname, "w") as fp:
        for norad_cat_id in norad_cat_ids:
            lines = st.tle(iter_lines=True, norad_cat_id=norad_cat_id, orderby="epoch desc", format="3le")
            for line in lines:
                fp.write("%s\n" % line)

    return

def read_tles(fname):
    with open(fname, "r") as fp:
        lines = fp.readlines()

        tles = []
        for i in range(0, len(lines), 3):
            tles.append(twolineelement(lines[i], lines[i+1], lines[i+2]))


    return tles

if __name__ == "__main__":
    # Query
    intldes_query = "~~2010-058"
    
    # Set login
    with open("login.json", "r") as fp:
        settings = json.load(fp)

    # Download TLEs if absent
    if not os.path.exists("database.txt"):
        download_tles_from_spacetrack(settings, intldes_query, "database.txt")

    # Load TLEs
    tles = read_tles("database.txt")

    # Read TLEs
    sats = []
    for tle in tles:
        sats.append(twoline2rv(tle.tle1, tle.tle2, wgs84))

    # Extract parameters
    norad_cat_ids = [sat.satnum for sat in sats]
    a = np.array([sat.a-1 for sat in sats])*wgs84.radiusearthkm
    raan = np.array([sat.nodeo for sat in sats])*180.0/np.pi
    incl = np.array([sat.inclo for sat in sats])*180.0/np.pi
    argp = np.array([sat.argpo for sat in sats])*180.0/np.pi
    ma = np.array([sat.mo for sat in sats])*180.0/np.pi
    ecc = np.array([sat.ecco for sat in sats])
    q = (1.0 - ecc) * a
    Q = (1.0 + ecc) * a
    
    t = Time([sat.epoch for sat in sats])
    names = np.array([tle.name for tle in tles])
    desigs = np.array([tle.desig for tle in tles])

    
    fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(15, 10), sharex=True)

    for norad_cat_id in np.unique(norad_cat_ids):
        c = norad_cat_ids==norad_cat_id

        name, desig = names[c][0], desigs[c][0]
        ax1.plot_date(t[c].datetime, a[c], linestyle="-", marker=None, label="[%d/%s] %s" % (norad_cat_id, desig, name), alpha=0.9)

    ax1.legend()
    ax1.grid()
    #ax1.set_xlabel("Date (UTC)")
    ax1.set_ylabel("Orbital altitude (km)")

    for norad_cat_id in np.unique(norad_cat_ids):
        c = norad_cat_ids==norad_cat_id

        name, desig = names[c][0], desigs[c][0]
        ax2.plot_date(t[c].datetime, ecc[c], linestyle="-", marker=None, label="[%d/%s] %s" % (norad_cat_id, desig, name), alpha=0.9)

    #ax2.legend()
    ax2.grid()
    #ax2.set_xlabel("Date (UTC)")
    ax2.set_ylabel("Eccentricity")

    for norad_cat_id in np.unique(norad_cat_ids):
        c = norad_cat_ids==norad_cat_id

        name, desig = names[c][0], desigs[c][0]
        ax3.plot_date(t[c].datetime, incl[c], linestyle="-", marker=None, label="[%d/%s] %s" % (norad_cat_id, desig, name), alpha=0.9)

    #ax3.legend()
    ax3.grid()
    #ax3.set_xlabel("Date (UTC)")
    ax3.set_ylabel("Inclination (deg)")    

    for norad_cat_id in np.unique(norad_cat_ids):
        c = norad_cat_ids==norad_cat_id

        name, desig = names[c][0], desigs[c][0]
        ax4.plot_date(t[c].datetime, argp[c], linestyle="-", marker=None, label="[%d/%s] %s" % (norad_cat_id, desig, name), alpha=0.9)

    #ax4.legend()
    ax4.grid()
    ax4.set_xlabel("Date (UTC)")
    ax4.set_ylabel("Argument of perigee (deg)")    
    
    plt.tight_layout()
    plt.show()
    #plt.savefig("tle_swap.png")