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)