Handle TLE archives with multiple entries

pyorbital/orbital.py

@@ -160,11 +160,145 @@ 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:
+ warnings.warn("Found TLE data for %s that is %f days apart" %
+ (sdate, int(delta_days)))
+
+ # 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()
+ 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'):
+ 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)
 
@@ -174,6 +308,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
@@ -208,7 +343,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:
@@ -221,6 +357,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)
 

pyorbital/tests/TERRA.TLE

@@ -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

pyorbital/tests/test_orbital.py

@@ -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):
@@ -119,12 +121,59 @@ 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"))
+ # 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"""
+ 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(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)
+ 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 " \
-  ".00000075 00000-0 19863-4 0 9994"
+ ".00000075 00000-0 19863-4 0 9994"
  line2 = "2 24793 86.3994 209.3241 0002020 " \
-  "89.8714 270.2713 14.34246429 90794"
+ "89.8714 270.2713 14.34246429 90794"
 
  orb = orbital.Orbital('IRIDIUM 7 [+]', line1=line1, line2=line2)
  d = datetime(2018, 3, 7, 3, 30, 15)