diff --git a/level1c4pps/seviri2pps_lib.py b/level1c4pps/seviri2pps_lib.py index 9ea40da..5ddfb0f 100644 --- a/level1c4pps/seviri2pps_lib.py +++ b/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(' 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): diff --git a/level1c4pps/tests/test_seviri2pps.py b/level1c4pps/tests/test_seviri2pps.py index 59db6b1..519df61 100644 --- a/level1c4pps/tests/test_seviri2pps.py +++ b/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):