Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Handle TLE archives with multiple entries #73

Open
wants to merge 17 commits into
base: main
Choose a base branch
from
Open

Handle TLE archives with multiple entries #73

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
a5cd067
added TLE archive parsing
espg Jan 26, 2021
15995c0
added fixed error handleing; proper digest for multiple times
espg Jan 26, 2021
222cc12
addressing None case for tle_file...
espg Feb 4, 2021
2934a94
logic for user supplied tle lines on init
espg Feb 4, 2021
9c8728e
fixed if statement checking for Nonetype
espg Feb 4, 2021
0cf4cea
multiline string fix
espg Feb 4, 2021
c9e2e02
type casting for numpy timestamp to datetime
espg Feb 4, 2021
0518107
flake8 and other typo fixes...
espg Feb 4, 2021
da12178
moved absolute value check to right place
espg Feb 4, 2021
74a83ed
error message cleanup
espg Feb 4, 2021
1925548
Fixing style errors.
stickler-ci Feb 4, 2021
6a01ae1
added unit tests, fixed failing idiom test, other feedback...
espg Mar 7, 2022
4f6ace6
Fixing style errors.
stickler-ci Mar 7, 2022
d220092
fixed tests, new method for inferring utc for tle's that are specifie…
espg Mar 7, 2022
c13f58d
Fixing style errors.
stickler-ci Mar 7, 2022
33fa641
cleaning up after bot fucked up my working directory by somehow pushi…
espg Mar 7, 2022
b9e9bc5
Fixing style errors.
stickler-ci Mar 7, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
153 changes: 152 additions & 1 deletion pyorbital/orbital.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,11 +160,159 @@ class Orbital(object):
def __init__(self, satellite, tle_file=None, line1=None, line2=None):
satellite = satellite.upper()
self.satellite_name = satellite
self.tle_file = tle_file
self.tle = tlefile.read(satellite, tle_file=tle_file,
line1=line1, line2=line2)
if self.tle_file:
self.utctime = None
else:
self.utctime = self.tle2datetime64(float(self._tle.line1[18:32]))
self.orbit_elements = OrbitElements(self.tle)
self.sgdp4 = _SGDP4(self.orbit_elements)

@property
def tle(self):
# using an archive; select appropriate TLE from file
if self.tle_file:
tle_data = self.read_tle_file(self.tle_file)
dates = self.tle2datetime64(
np.array([float(line[18:32]) for line in tle_data[::2]]))

# set utctime to arbitrary value inside archive if object init
if self.utctime is None:
sdate = dates[-1]
else:
sdate = np.datetime64(self.utctime) # .timestamp() #self.utcs[0]
# Find index "iindex" such that dates[iindex-1] < sdate <= dates[iindex]
# Notes:
# 1. If sdate < dates[0] then iindex = 0
# 2. If sdate > dates[-1] then iindex = len(dates), beyond the right boundary!
iindex = np.searchsorted(dates, sdate)

if iindex in (0, len(dates)):
if iindex == len(dates):
# Reset index if beyond the right boundary (see note 2. above)
iindex -= 1
elif abs(sdate - dates[iindex - 1]) < abs(sdate - dates[iindex]):
# Choose the closest of the two surrounding dates
iindex -= 1

# Make sure the TLE we found is within the threshold
delta_days = abs(sdate - dates[iindex]) / np.timedelta64(1, 'D')
tle_thresh = 7
if delta_days > tle_thresh:
raise IndexError(
"Can't find tle data for %s within +/- %d days around %s" %
(self.satellite_name, tle_thresh, sdate))

if delta_days > 3:
<<<<<<< HEAD
warnings.warn("Found TLE data for %s that is %f days apart" %
espg marked this conversation as resolved.
Show resolved Hide resolved
(sdate, int(delta_days)))
espg marked this conversation as resolved.
Show resolved Hide resolved
=======
espg marked this conversation as resolved.
Show resolved Hide resolved
logging.warning("Found TLE data for %s that is %f days apart",
espg marked this conversation as resolved.
Show resolved Hide resolved
espg marked this conversation as resolved.
Show resolved Hide resolved
sdate, delta_days)
else:
espg marked this conversation as resolved.
Show resolved Hide resolved
logging.debug("Found TLE data for %s that is %f days apart",
espg marked this conversation as resolved.
Show resolved Hide resolved
sdate, delta_days)
>>>>>>> 4f6ace6cfb3b6668adc90ecfd70da102f1860d8e

# Select TLE data
tle1 = tle_data[iindex * 2]
tle2 = tle_data[iindex * 2 + 1]
self._tle = tlefile.read(self.satellite_name, tle_file=None,
line1=tle1, line2=tle2)

# Not using TLE archive;
# Either using current celestrek TLE,
# or using user supplied Line 1 and line 2
else:
sat_time = self.tle2datetime64(float(self._tle.line1[18:32]))
# Object just created
if not self.utctime:
self.utctime = datetime.now()
#if not np.isclose(sat_time.astype(int),
# np.datetime64(datetime.now(), 'ms'),
# atol=720*7*60*1000):
# self.utctime = datetime.now()
#else:
# self.utctime = sat_time.astype(datetime)
mismatch = sat_time.astype(datetime) - self.utctime
if abs(mismatch.days) > 30:
raise IndexError("""
Current TLE from celestrek is %d days newer than
requested orbital parameter. Please use TLE archive
for accurate results""" % (mismatch.days))
if abs(mismatch.days) > 7:
warnings.warn("""
Current TLE from celestrek is %d days newer than
requested orbital parameter. Please use TLE archive
for accurate results""" % (mismatch.days))

return self._tle

@tle.setter
def tle(self, new_tle):
self._tle = new_tle

@property
def orbit_elements(self):
return OrbitElements(self.tle)

@orbit_elements.setter
def orbit_elements(self, new_orbit_elements):
self._orbit_elements = new_orbit_elements

@property
def sgdp4(self):
return _SGDP4(self.orbit_elements)

@sgdp4.setter
def sgdp4(self, new_sgdp4):
self._sgdp4 = _SGDP4(self.orbit_elements)

@property
def utctime(self):
return self._utctime

@utctime.setter
def utctime(self, utc_time):
if not isinstance(utc_time, datetime):
times = np.array(utc_time, dtype='datetime64[m]')
if times.max() - times.min() > np.timedelta64(3, 'D'):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, could we check instead that both min and max times are within 3 days of the epoch instead? or do you think it would be too restrictive?

raise ValueError(
"Dates must not exceed 3 days")
utctime = np.array(times.astype(float).mean(),
dtype='datetime64[m]').astype(datetime)
self._utctime = utctime.tolist()
else:
self._utctime = utc_time

@staticmethod
def tle2datetime64(times):
"""Convert TLE timestamps to numpy.datetime64.
Args:
times (float): TLE timestamps as %y%j.1234, e.g. 18001.25
"""
# Convert %y%j.12345 to %Y%j.12345 (valid for 1950-2049)
times = np.where(times > 50000, times + 1900000, times + 2000000)

# Convert float to datetime64
doys = (times % 1000).astype('int') - 1
years = (times // 1000).astype('int')
msecs = np.rint(24 * 3600 * 1000 * (times % 1))
times64 = (
years - 1970).astype('datetime64[Y]').astype('datetime64[ms]')
times64 += doys.astype('timedelta64[D]')
times64 += msecs.astype('timedelta64[ms]')

return times64

def read_tle_file(self, tle_filename):
"""Read TLE file."""
with open(tle_filename, 'r') as fp_:
return fp_.readlines()

def __str__(self):
return self.satellite_name + " " + str(self.tle)

Expand All @@ -174,6 +322,7 @@ def get_last_an_time(self, utc_time):
"""

# Propagate backwards to ascending node
self.utctime = utc_time
dt = np.timedelta64(10, 'm')
t_old = utc_time
t_new = t_old - dt
Expand Down Expand Up @@ -208,7 +357,8 @@ def get_position(self, utc_time, normalize=True):
"""Get the cartesian position and velocity from the satellite.
"""

kep = self._sgdp4.propagate(utc_time)
self.utctime = utc_time
kep = self.sgdp4.propagate(utc_time)
pos, vel = kep2xyz(kep)

if normalize:
Expand All @@ -221,6 +371,7 @@ def get_lonlatalt(self, utc_time):
"""Calculate sublon, sublat and altitude of satellite.
http://celestrak.com/columns/v02n03/
"""
self.utctime = utc_time
(pos_x, pos_y, pos_z), (vel_x, vel_y, vel_z) = self.get_position(
utc_time, normalize=True)

Expand Down
38 changes: 38 additions & 0 deletions pyorbital/tests/TERRA.TLE
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
1 25994U 99068A 15002.10622162 -.00010290 00000-0 -22790-2 0 9759
2 25994 98.2074 78.9859 0001342 97.3908 262.7439 14.57079636800035
1 25994U 99068A 15002.44957755 -.00011270 00000-0 -24979-2 0 9766
2 25994 98.2075 79.3243 0001324 97.4407 262.6918 14.57070694800084
1 25994U 99068A 15002.72426476 -.00011562 00000-0 -25635-2 0 9773
2 25994 98.2073 79.5948 0001388 99.3169 260.8205 14.57064107800128
1 25994U 99068A 15003.06761484 -.00006057 00000-0 -13365-2 0 9783
2 25994 98.2076 79.9333 0001481 95.1887 264.9645 14.57077442800172
1 25994U 99068A 15003.41096220 -.00000952 00000-0 -20131-3 0 9790
2 25994 98.2075 80.2718 0001425 95.6926 264.4504 14.57088545800221
1 25994U 99068A 15003.75431301 .00000617 00000-0 14694-3 0 9805
2 25994 98.2077 80.6095 0001466 92.3976 267.7414 14.57092512800279
1 25994U 99068A 15004.09766505 .00000916 00000-0 21346-3 0 9815
2 25994 98.2075 80.9478 0001466 91.5587 268.5799 14.57093702800327
1 25994U 99068A 15004.44101593 .00001079 00000-0 24956-3 0 9826
2 25994 98.2073 81.2861 0001448 91.1265 269.0094 14.57095042800372
1 25994U 99068A 15004.71569724 .00001055 00000-0 24423-3 0 9834
2 25994 98.2074 81.5566 0001441 89.4374 270.6982 14.57095542800416
1 25994U 99068A 15004.92170827 .00000907 00000-0 21141-3 0 9849
2 25994 98.2073 81.7597 0001439 87.7970 272.3361 14.57095246800445
1 25994U 99068A 15005.33374466 -.00002833 00000-0 00000+0 0 9859
2 25994 98.2069 82.1670 0001390 96.0356 264.1502 14.57082143800509
1 25994U 99068A 15005.40240164 .00000752 00000-0 17698-3 0 9867
2 25994 98.2073 82.2333 0001440 88.4447 271.6914 14.57095249800515
1 25994U 99068A 15005.67708292 .00000681 00000-0 16123-3 0 9871
2 25994 98.2072 82.5039 0001446 87.3999 272.7350 14.57095313800558
1 25994U 99068A 15006.08910599 .00000506 00000-0 12235-3 0 9889
2 25994 98.2072 82.9098 0001450 87.6000 272.5367 14.57094811800613
1 25994U 99068A 15006.22644661 .00000501 00000-0 12122-3 0 9890
2 25994 98.2073 83.0451 0001446 87.8459 272.2908 14.57094955800633
1 25994U 99068A 15006.36380590 -.00003936 00000-0 00000+0 0 9909
2 25994 98.2071 83.1828 0001316 94.3822 265.8029 14.57077694800658
1 25994U 99068A 15018.86171352 .00000465 00000-0 11317-3 0 272
2 25994 98.2045 95.4921 0001433 89.4794 270.6562 14.57110163802479
1 25994U 99068A 15019.20506188 .00000424 00000-0 10426-3 0 286
2 25994 98.2046 95.8304 0001428 91.1799 268.9572 14.57110302802520
1 25994U 99068A 15019.54840962 .00000410 00000-0 10107-3 0 296
2 25994 98.2044 96.1687 0001439 93.7627 266.3745 14.57110648802578
85 changes: 75 additions & 10 deletions pyorbital/tests/test_orbital.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,12 @@
from datetime import datetime, timedelta

import numpy as np
import os

from pyorbital import orbital

eps_deg = 10e-3
_DATAPATH = os.path.dirname(os.path.abspath(__file__))


class Test(unittest.TestCase):
Expand Down Expand Up @@ -119,32 +121,95 @@ def test_orbit_num_equator(self):
self.assertTrue(pos1[2] < 0)
self.assertTrue(pos2[2] > 0)

def test_error_four_day_interpolation(self):
"""Tests for error when time range exceeds three days"""
sat = orbital.Orbital("EOS-TERRA",
os.path.join(_DATAPATH, "./TERRA.TLE"))
espg marked this conversation as resolved.
Show resolved Hide resolved
# 720 minutes / 12 hours, times 5 to give 4.5 day search window
search=720*5
date=datetime(2015, 1, 25, 12)
time=np.array(date, dtype='datetime64[m]')
window=(time - search) + np.arange(search*2)
with self.assertRaises(Exception):
sat.get_lonlatalt(window)

def test_no_error_two_day_interpolation(self):
"""Tests for list of times, when list is under three days"""
<<<<<<< HEAD
sat = orbital.Orbital("EOS-TERRA",
espg marked this conversation as resolved.
Show resolved Hide resolved
os.path.join(_DATAPATH, "./TERRA.TLE"))
espg marked this conversation as resolved.
Show resolved Hide resolved
search=720*2
espg marked this conversation as resolved.
Show resolved Hide resolved
date = datetime(2015,1,25,12)
espg marked this conversation as resolved.
Show resolved Hide resolved
time = np.array(date, dtype='datetime64[m]')
window = (time - search) + np.arange(search*2)
self.assertEqual(len(sat.get_lonlatalt(window)[0]), 2880)

def test_warn_four_day_projection(self):
"""Tests for warning when TLE's are stale but still usable"""
sat = orbital.Orbital("EOS-TERRA",
os.path.join(_DATAPATH, "./TERRA.TLE"))
date = datetime(2015,1,26,12)
=======
sat=orbital.Orbital("EOS-TERRA",
os.path.join(_DATAPATH, "./TERRA.TLE")
search=720*2
date=datetime(2015, 1, 25, 12)
time=np.array(date, dtype='datetime64[m]')
window=(time - search) + np.arange(search*2)
self.assertEqual(sat.get_lonlatalt(window)[0], 2880)

def test_warn_four_day_projection(self):
"""Tests for warning when TLE's are stale but still usable"""
sat=orbital.Orbital("EOS-TERRA",
os.path.join(_DATAPATH, "./TERRA.TLE")
date=datetime(2015, 1, 25, 12)
>>>>>>> 4f6ace6cfb3b6668adc90ecfd70da102f1860d8e
with self.assertWarns(Warning):
sat.get_lonlatalt(date)

def test_error_ten_day_projection(self):
"""Tests for error on large TLE gap with no TLE file"""
sat=orbital.Orbital("EOS-TERRA")
# default behavior is to grab current TLE from celestrek
date=datetime(2011, 5, 1, 12)
with self.assertRaises(Exception):
sat.get_lonlatalt(date)

def test_error_ten_day_projection_file(self):
"""Tests for error on large TLE gap with TLE file"""
sat=orbital.Orbital("EOS-TERRA",
os.path.join(_DATAPATH, "./TERRA.TLE"))
# same as above, but with a local file
date=datetime(2011, 5, 1, 12)
with self.assertRaises(Exception):
sat.get_lonlatalt(date)

def test_get_next_passes_apogee(self):
"""Regression test #22."""
line1 = "1 24793U 97020B 18065.48735489 " \
line1="1 24793U 97020B 18065.48735489 " \
".00000075 00000-0 19863-4 0 9994"
line2 = "2 24793 86.3994 209.3241 0002020 " \
line2="2 24793 86.3994 209.3241 0002020 " \
"89.8714 270.2713 14.34246429 90794"

orb = orbital.Orbital('IRIDIUM 7 [+]', line1=line1, line2=line2)
d = datetime(2018, 3, 7, 3, 30, 15)
res = orb.get_next_passes(d, 1, 170.556, -43.368, 0.5, horizon=40)
orb=orbital.Orbital('IRIDIUM 7 [+]', line1=line1, line2=line2)
d=datetime(2018, 3, 7, 3, 30, 15)
res=orb.get_next_passes(d, 1, 170.556, -43.368, 0.5, horizon=40)
self.assertTrue(abs(
res[0][2] - datetime(2018, 3, 7, 3, 48, 13, 178439)) <
timedelta(seconds=0.01))

def test_get_next_passes_tricky(self):
""" Check issue #34 for reference """
line1 = "1 43125U 18004Q 18251.42128650 " \
line1="1 43125U 18004Q 18251.42128650 " \
"+.00001666 +00000-0 +73564-4 0 9991"

line2 = "2 43125 097.5269 314.3317 0010735 "\
line2="2 43125 097.5269 314.3317 0010735 "\
"157.6344 202.5362 15.23132245036381"

orb = orbital.Orbital('LEMUR-2-BROWNCOW', line1=line1, line2=line2)
d = datetime(2018, 9, 8)
orb=orbital.Orbital('LEMUR-2-BROWNCOW', line1=line1, line2=line2)
d=datetime(2018, 9, 8)

res = orb.get_next_passes(d, 72, -8.174163, 51.953319, 0.05, horizon=5)
res=orb.get_next_passes(d, 72, -8.174163, 51.953319, 0.05, horizon=5)

self.assertTrue(abs(
res[0][2] - datetime(2018, 9, 8, 9, 5, 46, 375248)) <
Expand Down