Refactoring

Signed-off-by: Adam.Dybbroe <[email protected]>

pyorbital/sno_utils.py

@@ -107,6 +107,11 @@ def __init__(self, platform_id, calipso_id, time_window, sno_min_thr, arc_len_mi
 
         self._check_platforms()
 
+        self.tle_id_platform_one = get_satellite_catalogue_number_from_name(self.calipso_id)
+        self.tle_id_platform_other = get_satellite_catalogue_number_from_name(self.platform_id)
+
+        self._minthr_step = 20  # min less than half an orbit
+
     def _check_platforms(self):
         # Check if satellite is supported:
         if OSCAR_NAMES.get(self.platform_id, self.platform_id).upper() not in SATELLITES.keys():
@@ -160,21 +165,23 @@ def get_snos_within_time_window(self):
             if filename_other.exists():
                 break
 
-        TLE_ID_CALIPSO = get_satellite_catalogue_number_from_name(self.calipso_id)
-        TLE_ID_OTHER = get_satellite_catalogue_number_from_name(self.platform_id)
-
-        minthr_step = 20  # min less than half an orbit probably
-        dtime = timedelta(seconds=60 * minthr_step * 2.0)
-        timestep_double = timedelta(seconds=60 * minthr_step * 2.0)
-        # make sure the two sat pass the SNO in the same step. We need and overlap of at least half minthr minutes.
-        timestep_plus_30s = timedelta(seconds=60 * minthr_step * 1.0 + (self.sno_minute_threshold*0.5)*60 + 30)
+        dtime = timedelta(seconds=60 * self._minthr_step * 2.0)
+        timestep_double = timedelta(seconds=60 * self._minthr_step * 2.0)
 
         calipso_obj = dt.datetime(1970, 1, 1).replace(tzinfo=timezone.utc)
         other_obj = dt.datetime(1970, 1, 1).replace(tzinfo=timezone.utc)
 
+        tle_finder_one = TwoLineElementsFinder(self.tle_id_platform_one, str(filename_calipso))
+        tle_finder_one.populate_tle_buffer()
+        self._tle_buffer_calipso = tle_finder_one.tle_buffer
+
+        tle_finder_other = TwoLineElementsFinder(self.tle_id_platform_other, str(filename_other))
+        tle_finder_other.populate_tle_buffer()
+        self._tle_buffer_other = tle_finder_other.tle_buffer
+
         tobj = self.time_start
         tle_calipso = None
-        tle_the_other_one = None
+        tle_the_other = None
         tobj_tmp = self.time_start
         i = 0
         t_diff = timedelta(days=1)
@@ -185,72 +192,34 @@ def get_snos_within_time_window(self):
                 message = time.time() - tic, "seconds", dtime
                 print(message)
 
-            if not tle_calipso or calipso_obj - tobj > t_diff or tobj - calipso_obj > t_diff:
-                tle_finder = TwoLineElementsFinder(TLE_ID_CALIPSO, str(filename_calipso),
-                                                   tle_buffer=self._tle_buffer_calipso)
-                tle_calipso = tle_finder.get_tle_archive(tobj)
-                self._tle_buffer_calipso = tle_finder.tle_buffer
-                # tle_calipso = get_tle_archive(tobj, str(filename_calipso), TLE_ID_CALIPSO, TLE_BUFFER_CALIPSO)
+            if not tle_calipso or abs(calipso_obj - tobj) > t_diff:
+                tle_calipso = tle_finder_one.get_best_tle_from_archive(tobj)
+                self._tle_buffer_calipso = tle_finder_one.tle_buffer
 
             if tle_calipso:
                 calipso_obj = get_datetime_from_tle(tle_calipso)
 
-            if not tle_the_other_one or other_obj - tobj > t_diff or tobj - other_obj > t_diff:
-                tle_finder = TwoLineElementsFinder(TLE_ID_OTHER, str(filename_other),
-                                                   tle_buffer=self._tle_buffer_other)
-                tle_the_other_one = tle_finder.get_tle_archive(tobj)
-                self._tle_buffer_other = tle_finder.tle_buffer
-                # tle_the_other_one = get_tle_archive(tobj, filename_other, TLE_ID_OTHER, TLE_BUFFER_OTHER)
-
-            if tle_the_other_one:
-                other_obj = get_datetime_from_tle(tle_the_other_one)
-
-            calipso = Orbital(self.calipso_id,
-                              line1=tle_calipso.line1,
-                              line2=tle_calipso.line2)
-            the_other_one = Orbital(self.platform_id,
-                                    line1=tle_the_other_one.line1,
-                                    line2=tle_the_other_one.line2)
-
-            got_intersection_acurate = False
-            arc_calipso_vector = get_arc_vector(tobj, timestep_plus_30s, calipso, self.arc_len_min)
-            arc_the_other_one_vector = get_arc_vector(tobj, timestep_plus_30s, the_other_one, self.arc_len_min)
-            # Approximate tracks with one arc each self.arc_len_min minutes.
-            # For each pair of arcs check if they intersect.
-            # There is atmost one intersection. Quit when we find it.
-            for arc_calipso in arc_calipso_vector:
-                for arc_the_other_one in arc_the_other_one_vector:
-                    if arc_calipso.intersects(arc_the_other_one):
-                        got_intersection_acurate = True
-                    if got_intersection_acurate:
-                        break
-                if got_intersection_acurate:
-                    break
-
-            if got_intersection_acurate:
-                sno = get_sno_point(calipso, the_other_one,
-                                    arc_calipso, arc_the_other_one,
+            if not tle_the_other or abs(other_obj - tobj) > t_diff:
+                tle_the_other = tle_finder_other.get_best_tle_from_archive(tobj)
+                self._tle_buffer_other = tle_finder_other.tle_buffer
+
+            if tle_the_other:
+                other_obj = get_datetime_from_tle(tle_the_other)
+
+            calipso = Orbital(self.calipso_id, line1=tle_calipso.line1, line2=tle_calipso.line2)
+            the_other = Orbital(self.platform_id, line1=tle_the_other.line1, line2=tle_the_other.line2)
+
+            retv = self.search_intersection(tobj, calipso, the_other)
+            if retv:
+                arc_calipso, arc_the_other = retv
+                sno = get_sno_point(calipso, the_other,
+                                    arc_calipso, arc_the_other,
                                     tobj, self.sno_minute_threshold, self.station)
 
                 if sno:
-                    # For debugging:
-                    create_geojson_line('./calipso_arc_%d.geojson' % i, arc_calipso)
+                    self.write_output_on_screen(sno, arc_calipso, arc_the_other, i)
                     results.append(sno)
 
-                    seconds_a = int(sno['satAdatetime'].strftime("%S")) + \
-                        float(sno['satAdatetime'].strftime("%f"))/1000000.
-                    seconds_b = int(sno['satBdatetime'].strftime("%S")) + \
-                        float(sno['satBdatetime'].strftime("%f"))/1000000.
-
-                    print("  " +
-                          str(sno['satBdatetime'].strftime("%Y%m%d %H:%M")) +
-                          "%5.1fs" % seconds_b + " "*5 +
-                          str(sno['satAdatetime'].strftime("%Y%m%d %H:%M")) +
-                          "%5.1fs" % seconds_a +
-                          " "*6 + "(%7.2f, %7.2f)" % (sno['sno_latitude'], sno['sno_longitude']) +
-                          "   " + "%4.1f min" % (sno['minutes_diff']) + "   " + str(sno['within_local_reception_area'])
-                          )
-
             tobj = tobj + timestep_double
             if tobj - tobj_tmp > timedelta(days=1):
                 tobj_tmp = tobj
@@ -259,6 +228,42 @@ def get_snos_within_time_window(self):
         print(str(results[0]))
         return pd.DataFrame(results)
 
+    def search_intersection(self, dtime, orb_obj1, orb_obj2):
+        """Get two pieces of a great arc circle and check if there is an intersection."""
+        # make sure the two sat pass the SNO in the same step. We need and overlap of at least half minthr minutes.
+        timestep_plus_30s = timedelta(seconds=60 * self._minthr_step * 1.0 + (self.sno_minute_threshold*0.5)*60 + 30)
+
+        arc_orb_obj1_vector = get_arc_vector(dtime, timestep_plus_30s, orb_obj1, self.arc_len_min)
+        arc_orb_obj2_vector = get_arc_vector(dtime, timestep_plus_30s, orb_obj2, self.arc_len_min)
+        # Approximate tracks with one arc each self.arc_len_min minutes.
+        # For each pair of arcs check if they intersect.
+        # There is atmost one intersection. Quit when we find it.
+        for arc_orb_obj1 in arc_orb_obj1_vector:
+            for arc_orb_obj2 in arc_orb_obj2_vector:
+                if arc_orb_obj1.intersects(arc_orb_obj2):
+                    return arc_orb_obj1, arc_orb_obj2
+
+        return
+
+    def write_output_on_screen(self, sno, arc_calipso, arc_the_other, idx):
+        """Write output results on stdout."""
+        # For debugging:
+        create_geojson_line('./calipso_arc_%d.geojson' % idx, arc_calipso)
+
+        seconds_a = int(sno['satAdatetime'].strftime("%S")) + \
+            float(sno['satAdatetime'].strftime("%f"))/1000000.
+        seconds_b = int(sno['satBdatetime'].strftime("%S")) + \
+            float(sno['satBdatetime'].strftime("%f"))/1000000.
+
+        print("  " +
+              str(sno['satBdatetime'].strftime("%Y%m%d %H:%M")) +
+              "%5.1fs" % seconds_b + " "*5 +
+              str(sno['satAdatetime'].strftime("%Y%m%d %H:%M")) +
+              "%5.1fs" % seconds_a +
+              " "*6 + "(%7.2f, %7.2f)" % (sno['sno_latitude'], sno['sno_longitude']) +
+              "   " + "%4.1f min" % (sno['minutes_diff']) + "   " + str(sno['within_local_reception_area'])
+              )
+
 
 def get_satellite_catalogue_number_from_name(satname):
     """Return the satellite catalogue number from the platform name."""
@@ -303,24 +308,21 @@ def get_arc_vector(timeobj, delta_t, sat, arc_len_min):
     return arcs
 
 
-def get_sno_point(calipso, the_other_one, arc_calipso, arc_the_other_one, tobj, minthr, station):
+def get_sno_point(calipso, the_other, arc_calipso, arc_the_other, tobj, minthr, station):
     """Get the SNO point if there is any.
 
     If the two sub-satellite tracks of the overpasses intersects
     get the sub-satellite position and time where they cross,
     and determine if the time deviation is smaller than the require threshold:
     """
     import math
-    intersect = arc_calipso.intersection(arc_the_other_one)
+    intersect = arc_calipso.intersection(arc_the_other)
     point = (math.degrees(intersect.lon),
              math.degrees(intersect.lat))
-    nextp = the_other_one.get_next_passes(tobj - timedelta(seconds=60*60),
-                                          # SNO around tobj check for passes between +- one hour.
-                                          # So that wanted pass for sure is next pass!
-                                          2,  # Number of hours to find overpasses
-                                          point[0],
-                                          point[1],
-                                          0)
+    # SNO around tobj check for passes between +- one hour.
+    # So that wanted pass for sure is next pass!
+    nextp = the_other.get_next_passes(tobj - timedelta(seconds=60*60), 2,  # Number of hours to find overpasses
+                                      point[0], point[1], 0)
 
     minthr_step = 20  # min less than half an orbit probably
     dtime = timedelta(seconds=60 * minthr_step * 2.0)
@@ -347,8 +349,8 @@ def get_sno_point(calipso, the_other_one, arc_calipso, arc_the_other_one, tobj,
     # Get observer look from Norrkoping to the satellite when it is
     # in zenith over the SNO point:
 
-    azi, elev = the_other_one.get_observer_look(maxt,
-                                                station['lon'], station['lat'], station['alt'])
+    azi, elev = the_other.get_observer_look(maxt,
+                                            station['lon'], station['lat'], station['alt'])
     isNorrk = (elev > 0.0)
 
     tdelta = (maxt_calipso - maxt)

pyorbital/tests/conftest.py

@@ -69,6 +69,38 @@ def fake_yamlconfig_file(tmp_path):
 2 37849  98.7741 305.5505 0001426 126.5558   7.5558 14.19528199113326
 """
 
+TEST_TLE_FILE2_CALIPSO = """1 29108U 06016B   13365.56860824  .00000817  00000-0  19114-3 0  9995
+2 29108 098.2183 306.3131 0001262 081.5858 278.5464 14.57147103408355
+1 29108U 06016B   14001.73596105  .00000831  00000-0  19433-3 0  9990
+2 29108 098.2189 307.4643 0001241 081.3032 278.8307 14.57149241408523
+1 29108U 06016B   14002.83464415  .00000945  00000-0  21969-3 0  9993
+2 29108 098.2193 308.5492 0001225 081.5021 278.6342 14.57151825408689
+1 29108U 06016B   14004.00199336 -.00003567  00000-0 -78180-3 0  9993
+2 29108 098.2205 309.7016 0001599 074.2235 285.7542 14.57124395408850
+1 29108U 06016B   14004.34536874 -.00007900  00000-0 -17461-2 0  9991
+2 29108 098.2198 310.0372 0001753 016.4256 343.5560 14.57094406408902
+1 29108U 06016B   14004.75741906 -.00011581  00000-0 -25675-2 0  9996
+2 29108 098.2201 310.4470 0001148 111.5107 248.5630 14.57068459408967
+1 29108U 06016B   14005.10078897 -.00014693  00000-0 -32632-2 0  9996
+2 29108 098.2210 310.7843 0001695 087.8266 272.2525 14.57045430409010
+1 29108U 06016B   14006.26824747 -.00016004  00000-0 -35589-2 0  9999
+2 29108 098.2212 311.9389 0000985 074.8347 285.2961 14.57006491409187
+1 29108U 06016B   14006.61162114 -.00010721  00000-0 -23779-2 0  9999
+2 29108 098.2215 312.2750 0001785 054.7782 305.4877 14.57027072409230
+1 29108U 06016B   14006.95496482 -.00007161  00000-0 -15838-2 0  9995
+2 29108 098.2196 312.6152 0001053 142.1452 218.0627 14.57036401409285
+1 29108U 06016B   14007.36698875  .00001082  00000-0  25039-3 0  9990
+2 29108 098.2203 313.0220 0001243 089.7863 270.4344 14.57065316409340
+1 29108U 06016B   14008.46571865  .00001229  00000-0  28299-3 0  9995
+2 29108 098.2204 314.1066 0001330 084.9063 275.2275 14.57068253409504
+1 29108U 06016B   14009.56446154  .00001213  00000-0  27957-3 0  9992
+2 29108 098.2198 315.1909 0001273 090.7695 269.3617 14.57071150409666
+1 29108U 06016B   14010.73187492  .00001212  00000-0  27917-3 0  9995
+2 29108 098.2194 316.3431 0001281 087.5323 272.5997 14.57074056409830
+1 29108U 06016B   14011.83061533  .00001150  00000-0  26558-3 0  9995
+2 29108 098.2191 317.4273 0001267 090.5743 269.5612 14.57076467409992
+"""
+
 
 @pytest.fixture
 def fake_tle_file1_calipso(tmp_path):
@@ -80,6 +112,16 @@ def fake_tle_file1_calipso(tmp_path):
     yield file_path
 
 
+@pytest.fixture
+def fake_tle_file2_calipso(tmp_path):
+    """Write fake TLE file for Calipso."""
+    file_path = tmp_path / 'TLE_calipso.txt'
+    with open(file_path, 'w') as fpt:
+        fpt.write(TEST_TLE_FILE2_CALIPSO)
+
+    yield file_path
+
+
 @pytest.fixture
 def fake_tle_file1_snpp(tmp_path):
     """Write fake TLE file for Suomi-NPP."""

pyorbital/tests/test_tle_archive.py

@@ -24,8 +24,6 @@
 import datetime as dt
 from datetime import timezone
 import pyorbital
-# from pyorbital.tle_archive import populate_tle_buffer
-# from pyorbital.tle_archive import get_tle_archive
 from pyorbital.tle_archive import TwoLineElementsFinder
 
 
@@ -54,14 +52,14 @@ def test_populate_tle_buffer(fake_tle_file1_calipso):
         assert tlelines[idx*2+1].strip() == tleobj.line2
 
 
-def test_get_tle_archive(fake_tle_file1_calipso):
+def test_get_best_tle_from_archive(fake_tle_file1_calipso):
     """Test getting all the TLEs from file with many TLEs."""
     tle_filename = str(fake_tle_file1_calipso)
     tleid_calipso = '29108'
     dtobj = dt.datetime(2014, 1, 2, tzinfo=timezone.utc)
 
     tle_finder = TwoLineElementsFinder(tleid_calipso, tle_filename)
-    tleobj = tle_finder.get_tle_archive(dtobj)
+    tleobj = tle_finder.get_best_tle_from_archive(dtobj)
 
     ts = (tleobj.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
     dtime_valid = dt.datetime.fromtimestamp(ts, tz=timezone.utc)
@@ -71,7 +69,7 @@ def test_get_tle_archive(fake_tle_file1_calipso):
 
     dtobj = dt.datetime(2014, 1, 1, 12, tzinfo=timezone.utc)
     tle_finder = TwoLineElementsFinder(tleid_calipso, tle_filename)
-    tleobj = tle_finder.get_tle_archive(dtobj)
+    tleobj = tle_finder.get_best_tle_from_archive(dtobj)
 
     ts = (tleobj.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
     dtime_valid = dt.datetime.fromtimestamp(ts, tz=timezone.utc)
@@ -81,10 +79,26 @@ def test_get_tle_archive(fake_tle_file1_calipso):
 
     dtobj = dt.datetime(2014, 1, 1, 16, tzinfo=timezone.utc)
     tle_finder = TwoLineElementsFinder(tleid_calipso, tle_filename)
-    tleobj = tle_finder.get_tle_archive(dtobj)
+    tleobj = tle_finder.get_best_tle_from_archive(dtobj)
 
     ts = (tleobj.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
     dtime_valid = dt.datetime.fromtimestamp(ts, tz=timezone.utc)
     expected = dt.datetime(2014, 1, 1, 17, 39, 47, 34720, tzinfo=timezone.utc)
 
     assert abs(expected - dtime_valid).total_seconds() < 0.001
+
+
+# def test_get_best_tle_from_archive(fake_tle_file2_calipso):
+#     """Test getting all the TLEs from file with many TLEs."""
+#     tle_filename = str(fake_tle_file2_calipso)
+#     tleid_calipso = '29108'
+#     dtobj = dt.datetime(2014, 1, 6, 12, tzinfo=timezone.utc)
+
+#     tle_finder = TwoLineElementsFinder(tleid_calipso, tle_filename)
+#     tleobj = tle_finder.get_best_tle_from_archive(dtobj)
+
+#     ts = (tleobj.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
+#     dtime_valid = dt.datetime.fromtimestamp(ts, tz=timezone.utc)
+#     expected = dt.datetime(2014, 1, 1, 17, 39, 47, 34720, tzinfo=timezone.utc)
+
+#     assert abs(expected - dtime_valid).total_seconds() < 0.001