Skip to content

Instantly share code, notes, and snippets.

@RalucaNicola
Created April 19, 2022 19:49
Show Gist options
  • Save RalucaNicola/3df6506d70f6a39dc28bb5dd95aab1b6 to your computer and use it in GitHub Desktop.
Save RalucaNicola/3df6506d70f6a39dc28bb5dd95aab1b6 to your computer and use it in GitHub Desktop.
Orbital computations
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment