Merge pull request #907 from sfinkens/fix-satpos-seviri-hrit

Handle bad orbit coefficients in SEVIRI HRIT header

satpy/readers/seviri_l1b_hrit.py

@@ -335,16 +335,28 @@ def _find_orbit_coefs(self):
         Returns: Corresponding index in the coefficient list.
 
         """
-        # Find index of interval enclosing the nominal timestamp of the scan
         time = np.datetime64(self.prologue['ImageAcquisition']['PlannedAcquisitionTime']['TrueRepeatCycleStart'])
         intervals_tstart = self.prologue['SatelliteStatus']['Orbit']['OrbitPolynomial']['StartTime'][0].astype(
             'datetime64[us]')
         intervals_tend = self.prologue['SatelliteStatus']['Orbit']['OrbitPolynomial']['EndTime'][0].astype(
             'datetime64[us]')
         try:
-            return np.where(np.logical_and(time >= intervals_tstart, time < intervals_tend))[0][0]
-        except IndexError:
-            raise NoValidOrbitParams('Unable to find orbit coefficients valid for {}'.format(time))
+            # Find index of interval enclosing the nominal timestamp of the scan. If there are
+            # multiple enclosing intervals, use the most recent one.
+            enclosing = np.where(np.logical_and(time >= intervals_tstart, time < intervals_tend))[0]
+            most_recent = np.argmax(intervals_tstart[enclosing])
+            return enclosing[most_recent]
+        except ValueError:
+            # No enclosing interval. Instead, find the interval whose centre is closest to the scan's timestamp
+            # (but not more than 6 hours apart)
+            intervals_centre = intervals_tstart + 0.5 * (intervals_tend - intervals_tstart)
+            diffs_us = (time - intervals_centre).astype('i8')
+            closest_match = np.argmin(np.fabs(diffs_us))
+            if abs(intervals_centre[closest_match] - time) < np.timedelta64(6, 'h'):
+                logger.warning('No orbit coefficients valid for {}. Using closest match.'.format(time))
+                return closest_match
+            else:
+                raise NoValidOrbitParams('Unable to find orbit coefficients valid for {}'.format(time))
 
     def get_earth_radii(self):
         """Get earth radii from prologue.
@@ -513,9 +525,10 @@ def _get_header(self):
         self.mda['orbital_parameters']['satellite_nominal_longitude'] = self.prologue['SatelliteStatus'][
             'SatelliteDefinition']['NominalLongitude']
         self.mda['orbital_parameters']['satellite_nominal_latitude'] = 0.0
-        self.mda['orbital_parameters']['satellite_actual_longitude'] = actual_lon
-        self.mda['orbital_parameters']['satellite_actual_latitude'] = actual_lat
-        self.mda['orbital_parameters']['satellite_actual_altitude'] = actual_alt
+        if actual_lon is not None:
+            self.mda['orbital_parameters']['satellite_actual_longitude'] = actual_lon
+            self.mda['orbital_parameters']['satellite_actual_latitude'] = actual_lat
+            self.mda['orbital_parameters']['satellite_actual_altitude'] = actual_alt
 
         # Misc
         self.platform_id = self.prologue["SatelliteStatus"][

satpy/tests/reader_tests/test_seviri_l1b_hrit.py

@@ -41,7 +41,7 @@
 
 def new_get_hd(instance, hdr_info):
     """Generate some metadata."""
-    instance.mda = {'spectral_channel_id': 'bla'}
+    instance.mda = {'spectral_channel_id': 1}
     instance.mda.setdefault('number_of_bits_per_pixel', 10)
 
     instance.mda['projection_parameters'] = {'a': 6378169.00,
@@ -534,6 +534,18 @@ def test_get_timestamps(self):
         self.assertTrue(np.all(day[1:-1] == 3))
         self.assertTrue(np.all(msec[1:-1] == np.arange(len(tline) - 2)))
 
+    def test_get_header(self):
+        # Make sure that the actual satellite position is only included if available
+        self.reader.mda['orbital_parameters'] = {}
+        self.reader.prologue_.get_satpos.return_value = 1, 2, 3
+        self.reader._get_header()
+        self.assertIn('satellite_actual_longitude', self.reader.mda['orbital_parameters'])
+
+        self.reader.mda['orbital_parameters'] = {}
+        self.reader.prologue_.get_satpos.return_value = None, None, None
+        self.reader._get_header()
+        self.assertNotIn('satellite_actual_longitude', self.reader.mda['orbital_parameters'])
+
 
 class TestHRITMSGPrologueFileHandler(unittest.TestCase):
     """Test the HRIT prologue file handler."""
@@ -560,9 +572,11 @@ def setUp(self, *mocks):
                 'Orbit': {
                     'OrbitPolynomial': {
                         'StartTime': np.array([
-                            [datetime(2006, 1, 1, 6), datetime(2006, 1, 1, 12), datetime(2006, 1, 1, 18)]]),
+                            [datetime(2006, 1, 1, 6), datetime(2006, 1, 1, 12), datetime(2006, 1, 1, 18),
+                             datetime(1958, 1, 1, 0)]]),
                         'EndTime': np.array([
-                            [datetime(2006, 1, 1, 12), datetime(2006, 1, 1, 18), datetime(2006, 1, 2, 0)]]),
+                            [datetime(2006, 1, 1, 12), datetime(2006, 1, 1, 18), datetime(2006, 1, 2, 0),
+                             datetime(1958, 1, 1, 0)]]),
                         'X': [np.zeros(8),
                               [8.41607082e+04, 2.94319260e+00, 9.86748617e-01, -2.70135453e-01,
                                -3.84364650e-02, 8.48718433e-03, 7.70548174e-04, -1.44262718e-04],
@@ -598,11 +612,60 @@ def init_patched(self, *args, **kwargs):
 
     def test_find_orbit_coefs(self):
         """Test identification of orbit coefficients."""
+        # Contiguous validity intervals (that's the norm)
         self.assertEqual(self.reader._find_orbit_coefs(), 1)
 
-        # No interval enclosing the given timestamp
-        self.reader.prologue['ImageAcquisition']['PlannedAcquisitionTime'][
-            'TrueRepeatCycleStart'] = datetime(2000, 1, 1)
+        # No interval enclosing the given timestamp ...
+        # a) closest interval should be selected (if not too far away)
+        self.reader.prologue['SatelliteStatus'] = {
+            'Orbit': {
+                'OrbitPolynomial': {
+                    'StartTime': np.array([
+                        [datetime(2006, 1, 1, 10), datetime(2006, 1, 1, 13)]]),
+                    'EndTime': np.array([
+                        [datetime(2006, 1, 1, 12), datetime(2006, 1, 1, 18)]])
+                }
+            }
+        }
+        self.assertEqual(self.reader._find_orbit_coefs(), 0)
+
+        # b) closest interval too far away
+        self.reader.prologue['SatelliteStatus'] = {
+            'Orbit': {
+                'OrbitPolynomial': {
+                    'StartTime': np.array([
+                        [datetime(2006, 1, 1, 0), datetime(2006, 1, 1, 18)]]),
+                    'EndTime': np.array([
+                        [datetime(2006, 1, 1, 4), datetime(2006, 1, 1, 22)]])
+                }
+            }
+        }
+        self.assertRaises(NoValidOrbitParams, self.reader._find_orbit_coefs)
+
+        # Overlapping intervals -> most recent interval should be selected
+        self.reader.prologue['SatelliteStatus'] = {
+            'Orbit': {
+                'OrbitPolynomial': {
+                    'StartTime': np.array([
+                        [datetime(2006, 1, 1, 6), datetime(2006, 1, 1, 10)]]),
+                    'EndTime': np.array([
+                        [datetime(2006, 1, 1, 13), datetime(2006, 1, 1, 18)]])
+                }
+            }
+        }
+        self.assertEqual(self.reader._find_orbit_coefs(), 1)
+
+        # No valid coefficients at all
+        self.reader.prologue['SatelliteStatus'] = {
+            'Orbit': {
+                'OrbitPolynomial': {
+                    'StartTime': np.array([
+                        [datetime(1958, 1, 1, 0), datetime(1958, 1, 1)]]),
+                    'EndTime': np.array([
+                        [datetime(1958, 1, 1, 0), datetime(1958, 1, 1)]])
+                }
+            }
+        }
         self.assertRaises(NoValidOrbitParams, self.reader._find_orbit_coefs)
 
     @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGPrologueFileHandler._find_orbit_coefs')