Interpolate broken timestamps

level1c4pps/seviri2pps_lib.py

@@ -35,21 +35,24 @@
 from trollsift.parser import globify, Parser
 from pyorbital.astronomy import get_alt_az, sun_zenith_angle
 from pyorbital.orbital import get_observer_look
-
+from level1c4pps import dt64_to_datetime
 from level1c4pps.calibration_coefs import get_calibration, CalibrationData
 from level1c4pps import (make_azidiff_angle,
+                         convert_angles,
+                         apply_sunz_correction,
                          get_encoding,
                          compose_filename,
                          update_angle_attributes,
                          fix_sun_earth_distance_correction_factor)
 
-
+import logging
 try:
     FileNotFoundError
 except NameError:
     # Python 2
     FileNotFoundError = IOError
 
+logger = logging.getLogger("seviri2pps")
 
 class UnexpectedSatpyVersion(Exception):
     """Exception if unexpected satpy version."""
@@ -115,7 +118,7 @@ def load_and_calibrate(filenames, rotate,
     _check_is_seviri_data(scn_)
     _load_bands(scn_, rotate)
     _update_scene_attrs(scn_, {'image_rotated': rotate})
-
+    scn_.attrs["filename_starttime"] = info['start_time']
     return scn_
 
 
@@ -267,6 +270,29 @@ def set_attrs(scene):
                      'platform_name', 'sensor', 'georef_offset_corrected']:
             scene[band].attrs.pop(attr, None)
 
+def interpolate_nats(acq_times, filename_starttime):
+    """Interpolate NaTs."""
+    # out_times.interpolate_na(dim = "y")  => can't cast array data from dtype('<M8[ns]') to dtype('float64')
+
+    update = np.isnat(acq_times.values)
+
+    too_late = [ind for (ind, acq_time_i) in enumerate(acq_times.values) if not np.isnat(acq_time_i)
+                and dt64_to_datetime(acq_time_i) - filename_starttime > timedelta(seconds=60*15)]
+    too_early = [ind for (ind, acq_time_i) in enumerate(acq_times.values) if not np.isnat(acq_time_i)
+                 and dt64_to_datetime(acq_time_i) - filename_starttime < timedelta(seconds=0) ]
+    update[too_late] = True
+    update[too_early] = True
+    x_val = np.array(list(range(0, len(acq_times))))
+    x_val_ok = x_val[~update]
+    y_val_ok = acq_times.values[~update].astype('float64')
+    index_to_update = [ind for ind in np.where(update)[0] if ind > 52 and ind < 3660]
+
+    if len(index_to_update) > 0:
+        logger.info(f'Interpolating timestamp at index {index_to_update}')
+        interp = np.interp(x_val, x_val_ok, y_val_ok)
+        acq_times[index_to_update] = dt64_to_datetime(interp[index_to_update])
+    return acq_times
+
 
 def get_mean_acq_time(scene):
     """Compute mean scanline acquisition time over all bands."""
@@ -284,7 +310,8 @@ def get_mean_acq_time(scene):
 
     # Compute average over all bands (skip NaNs)
     acq_times = xr.concat(acq_times, 'bands')
-    return acq_times.mean(dim='bands', skipna=True).astype(dtype)
+    out_times = acq_times.mean(dim='bands', skipna=True).astype(dtype)
+    return interpolate_nats(out_times, scene.attrs["filename_starttime"])
 
 
 def update_coords(scene):

level1c4pps/tests/test_seviri2pps.py

@@ -229,25 +229,34 @@ def test_set_attrs(self):
 
     def test_get_mean_acq_time(self):
         """Test computation of mean scanline acquisition time."""
+        scene = get_fake_scene()
         seviri2pps.BANDNAMES = ['VIS006', 'IR_108']
         vis006 = xr.DataArray(
-            data=[0, 0, 0],
+            data=list(range(56)),
             dims=('y', ),
-            coords={'acq_time': ('y', [None,
-                                       None,
-                                       dt.datetime(2009, 7, 1, 12, 1, 0)])})
+            coords={'acq_time': ('y', [None] * 52 + [None,
+                                                     np.nan,
+                                                     dt.datetime(2008, 7, 1, 12, 1, 0),
+                                                     dt.datetime(2009, 7, 1, 12, 1, 0)])})
         ir_108 = xr.DataArray(
-            data=[0, 0, 0],
+            data=list(range(56)),
             dims=('y', ),
-            coords={'acq_time': ('y', [None,
-                                       dt.datetime(2009, 7, 1, 12, 0, 30),
-                                       dt.datetime(2009, 7, 1, 12, 1, 30)])})
-        scene = {'VIS006': vis006, 'IR_108': ir_108}
+            coords={'acq_time': ('y', [None] * 52 + [dt.datetime(2009, 7, 1, 12, 0, 30),
+                                                      np.nan,
+                                                      dt.datetime(2008, 7, 1, 12, 1, 30),
+                                                      dt.datetime(2009, 7, 1, 12, 1, 30)])})
+        scene['VIS006']  = vis006
+        scene['IR_108'] = ir_108
+        scene.attrs["filename_starttime"] = dt.datetime(2009, 7, 1, 12, 0, 30)
+
         acq_exp = np.array(['NaT',
                             '2009-07-01 12:00:30',
+                            '2009-07-01 12:00:45',
+                            '2009-07-01 12:01:00',
                             '2009-07-01 12:01:15'], dtype='datetime64[s]')
         acq = seviri2pps.get_mean_acq_time(scene)
-        np.testing.assert_array_equal(acq, acq_exp)
+        print (acq)
+        np.testing.assert_array_equal(acq[-5:], acq_exp)
 
     @mock.patch('level1c4pps.seviri2pps_lib.get_mean_acq_time')
     def test_update_coords(self, get_mean_acq_time):