diff --git a/pyorbital/orbital.py b/pyorbital/orbital.py index f7e701b3..27fa1f6e 100644 --- a/pyorbital/orbital.py +++ b/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) diff --git a/pyorbital/tests/TERRA.TLE b/pyorbital/tests/TERRA.TLE new file mode 100644 index 00000000..3344cea0 --- /dev/null +++ b/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 diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index 9f4135b8..6c0ea7ba 100644 --- a/pyorbital/tests/test_orbital.py +++ b/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)