Orbital computations

main.py

import pandas as pd
from skyfield.api import load
from datetime import timedelta
from math import isnan, pi, sqrt, atan2, sin, cos
from geojson import FeatureCollection, Feature, LineString, dump
from arcgis.gis import GIS

#constants
ts = load.timescale()
now = ts.now()
deg2rad = pi / 180
rad2deg = 180 / pi


# finds the greenwich sidereal time
def gstime(date):
    # julian date in ut1 (days from 4713 bc)
    jdut1 = date.tt
    tut1 = (jdut1 - 2451545.0) / 36525.0

    temp = (-6.2e-6 * tut1 * tut1 * tut1) + (0.093104 * tut1 * tut1) + (((876600.0 * 3600) + 8640184.812866) * tut1) + 67310.54841
    temp = ((temp * deg2rad) / 240.0) % (pi * 2)

    if temp < 0.0:
        temp += pi * 2
    return temp


# converts a Earth-Centered Earth-Fixed coordinate system to Earth-Centered Inertial
def ecf_to_eci(ecf, gmst):
    x = (ecf[0] * cos(gmst)) - (ecf[1] * sin(gmst))
    y = (ecf[0] * (sin(gmst))) + (ecf[1] * cos(gmst))
    z = ecf[2]
    return [x, y, z]


# converts a Earth-Centered Inertial to geodetic
def eci_to_geodetic(eci, gmst):
    a = 6378.137
    b = 6356.7523142
    R = sqrt((eci[0] * eci[0]) + (eci[1] * eci[1]))
    f = (a - b) / a
    e2 = ((2 * f) - (f * f))

    longitude = atan2(eci[1], eci[0]) - gmst
    while longitude < -pi:
        longitude += pi * 2
    while longitude > pi:
        longitude -= pi * 2

    kmax = 20
    k = 0
    latitude = atan2(eci[2], sqrt((eci[0] * eci[0]) + (eci[1] * eci[1])))
    while k < kmax:
        C = 1 / sqrt(1 - (e2 * (sin(latitude) * sin(latitude))))
        latitude = atan2(eci[2] + (a * C * e2 * sin(latitude)), R)
        k += 1
    height = (R / cos(latitude)) - (a * C)
    coordinates = dict({'longitude': longitude, 'latitude': latitude, 'height': height})
    return coordinates


def get_location(satellite, date):
    position = satellite.at(date).xyz.km
    gmst = gstime(date)
    eci = ecf_to_eci(position, gmst)
    pos = eci_to_geodetic(eci, gmst)
    if isnan(pos['longitude']) or isnan(pos['latitude']):
        return None
    else:
        return pos['longitude'] * rad2deg, pos['latitude'] * rad2deg, pos['height'] * 1000


def get_orbit(satellite, period, start):
    segments = 200 if period > 1000 else 100 if period > 400 else 50
    milliseconds = (period * 60000) / segments
    vertices = []
    for index in range(0, segments + 1):
        date = start + timedelta(milliseconds=index * milliseconds)
        sat_location = get_location(satellite, date)
        if not sat_location:
            continue
        vertices.append(sat_location)
    return vertices


def get_attributes(norad):
    sat_attr = sat_meta.loc[norad]
    return {
        'norad': norad,
        'name': sat_attr['Name of Satellite, Alternate Names'],
        'official_name': sat_attr['Current Official Name of Satellite'],
        'country_operator': sat_attr['Country of Operator/Owner'],
        'operator': sat_attr['Operator/Owner'],
        'purpose': sat_attr['Purpose'],
        'orbit_class': sat_attr['Class of Orbit'],
        'orbit_type': sat_attr['Type of Orbit'] if not pd.isna(sat_attr['Type of Orbit']) else '',
        'perigee': sat_attr['Perigee (km)'],
        'apogee': sat_attr['Apogee (km)'],
        'launch_date': sat_attr['Date of Launch'],
        'launch_site': sat_attr['Launch Site'] if not pd.isna(sat_attr['Launch Site']) else '',
        'launch_vehicle': sat_attr['Launch Vehicle'],
        'cospar': sat_attr['COSPAR Number']
    }


# setting a threshold of 5 for checking similarity between orbital parameters
def are_similar(value1, value2):
    return abs(value1 - value2) < 5


def get_features(sat_tle, sat_meta):
    orbit_features = []
    unique_orbit_features = []
    for i in range(0, len(sat_tle) - 1):
        sat = sat_tle[i]
        norad = sat.model.satnum
        if norad in sat_meta.index:
        # if norad == 48865:
            period = (2 * pi) / sat.model.no
            inclination = sat.model.inclo * 180 / pi
            eccentricity = sat.model.ecco * 180 / pi
            perigee_argument = sat.model.argpo * 180 / pi
            node = sat.model.nodeo * 180 / pi
            display = 1
            for j in range(0, len(unique_orbit_features)):
                params = unique_orbit_features[j]
                if are_similar(params['inclination'], inclination) and are_similar(params['perigee_argument'], perigee_argument) and are_similar(params['node'], node):
                    display = 0
                    break
                # make an exception for Landsat 9
                if norad == 39084:
                    display = 1
            if display == 1:
                unique_orbit_features.append({
                    'inclination': inclination,
                    'perigee_argument': perigee_argument,
                    'node': node
                })
            if isinstance(period, float) and not isnan(period):
                orbit_coordinates = get_orbit(sat, period, now)
                print(norad)
                if len(orbit_coordinates) > 0:
                    orbit_properties = get_attributes(norad)
                    orbit_properties['period'] = period
                    orbit_properties['inclination'] = inclination
                    orbit_properties['eccentricity'] = eccentricity
                    orbit_properties['display'] = display
                    orbit_feature = Feature(geometry=LineString(orbit_coordinates), properties=orbit_properties)
                    orbit_features.append(orbit_feature)
    print('Total orbits: ', len(orbit_features))
    print('Unique orbits: ', len(unique_orbit_features))
    return orbit_features


if __name__ == '__main__':
    # read data
    sat_meta = pd.read_csv('./updated/sat_metadata_012022.csv')
    # check the Norad numbers which are not unique
    # in the current dataset there were about 10 - I manually corrected them
    print('Non-unique NORAD ids: ', sat_meta[sat_meta.duplicated(subset=['NORAD Number'], keep="first")]['NORAD Number'])
    sat_meta.set_index('NORAD Number', inplace=True)
    sat_meta.fillna('')
    # converts tle file into an iterable list of satrec information
    sat_tle = load.tle_file('./updated/norad-tle.txt')

    # generate orbits
    features = get_features(sat_tle, sat_meta)

    # save orbits as geojson
    with open('satellite_orbits.geojson', 'w') as f:
        dump(FeatureCollection(features), f)