Add scanline timestamps to seviri_l1b_hrit (#752)

Add scanline timestamps to seviri_l1b_hrit

satpy/readers/seviri_base.py

@@ -24,7 +24,6 @@
 """Utilities and eventually also base classes for MSG HRIT/Native data reading
 """
 
-from datetime import datetime, timedelta
 import numpy as np
 from numpy.polynomial.chebyshev import Chebyshev
 import dask.array as da
@@ -193,11 +192,30 @@
 
 
 def get_cds_time(days, msecs):
-    """Get the datetime object of the time since epoch given in days and
-    milliseconds of day
+    """Compute timestamp given the days since epoch and milliseconds of the day
+
+    1958-01-01 00:00 is interpreted as fill value and will be replaced by NaT (Not a Time).
+
+    Args:
+        days (int, either scalar or numpy.ndarray):
+            Days since 1958-01-01
+        msecs (int, either scalar or numpy.ndarray):
+            Milliseconds of the day
+
+    Returns:
+        numpy.datetime64: Timestamp(s)
     """
-    return datetime(1958, 1, 1) + timedelta(days=float(days),
-                                            milliseconds=float(msecs))
+    if np.isscalar(days):
+        days = np.array([days], dtype='int64')
+        msecs = np.array([msecs], dtype='int64')
+
+    time = np.datetime64('1958-01-01').astype('datetime64[ms]') + \
+        days.astype('timedelta64[D]') + msecs.astype('timedelta64[ms]')
+    time[time == np.datetime64('1958-01-01 00:00')] = np.datetime64("NaT")
+
+    if len(time) == 1:
+        return time[0]
+    return time
 
 
 def dec10216(inbuf):

satpy/readers/seviri_l1b_hrit.py

@@ -23,14 +23,108 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 """SEVIRI HRIT format reader
-******************************
+============================
+
+Introduction
+------------
+
+The ``seviri_l1b_hrit`` reader reads and calibrates MSG-SEVIRI L1.5 image data in HRIT format. The format is explained
+in the `MSG Level 1.5 Image Format Description`_. The files are usually named as
+follows:
+
+.. code-block:: none
+
+    H-000-MSG4__-MSG4________-_________-PRO______-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000001___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000002___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000003___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000004___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000005___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000006___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000007___-201903011200-__
+    H-000-MSG4__-MSG4________-IR_108___-000008___-201903011200-__
+    H-000-MSG4__-MSG4________-_________-EPI______-201903011200-__
+
+Each image is decomposed into 24 segments (files) for the high-resolution-visible (HRV) channel and 8 segments for other
+visible (VIS) and infrared (IR) channels. Additionally there is one prologue and one epilogue file for the entire scan
+which contain global metadata valid for all channels.
+
+Example
+-------
+Here is an example how to read the data in satpy:
+
+.. code-block:: python
+
+    from satpy import Scene
+    import glob
+
+    filenames = glob.glob('data/H-000-MSG4__-MSG4________-*201903011200*')
+    scn = Scene(filenames=filenames, reader='seviri_l1b_hrit')
+    scn.load(['VIS006', 'IR_108'])
+    print(scn['IR_108'])
+
+Output:
+
+.. code-block:: none
+
+    <xarray.DataArray 'reshape-5b8fc7364b289af7dec1f45b88ad2056' (y: 3712, x: 3712)>
+    dask.array<shape=(3712, 3712), dtype=float32, chunksize=(464, 3712)>
+    Coordinates:
+        acq_time  (y) datetime64[ns] NaT NaT NaT NaT NaT NaT ... NaT NaT NaT NaT NaT
+      * x         (x) float64 5.566e+06 5.563e+06 5.56e+06 ... -5.566e+06 -5.569e+06
+      * y         (y) float64 -5.566e+06 -5.563e+06 ... 5.566e+06 5.569e+06
+    Attributes:
+        satellite_longitude:      0.0
+        satellite_latitude:       0.0
+        satellite_altitude:       35785831.0
+        georef_offset_corrected:  True
+        wavelength:               (9.8, 10.8, 11.8)
+        units:                    K
+        standard_name:            brightness_temperature
+        sensor:                   seviri
+        navigation:               {'satellite_nominal_longitude': 0.0, 'satellite...
+        projection:               {'satellite_longitude': 0.0, 'satellite_latitud...
+        platform_name:            Meteosat-11
+        start_time:               2019-03-01 12:00:09.716000
+        end_time:                 2019-03-01 12:12:42.946000
+        area:                     Area ID: some_area_name\\nDescription: On-the-fl...
+        name:                     IR_108
+        resolution:               3000.403165817
+        calibration:              brightness_temperature
+        polarization:             None
+        level:                    None
+        modifiers:                ()
+        ancillary_variables:      []
+
+
+* The ``projection`` attribute specifies the projection parameters which are used to, for example, compute lat/lon
+  coordinates.
+* The ``navigation`` attribute holds the actual position of the satellite required for computing viewing
+  angles etc.
+* You can choose between nominal and GSICS calibration coefficients or even specify your own coefficients, see
+  :class:`HRITMSGFileHandler`.
+* The ``acq_time`` coordinate provides the acquisition time for each scanline. Use a ``MultiIndex`` to enable selection
+  by acquisition time:
+
+  .. code-block:: python
+
+      import pandas as pd
+      mi = pd.MultiIndex.from_arrays([scn['IR_108']['y'].data, scn['IR_108']['acq_time'].data],
+                                     names=('y_coord', 'time'))
+      scn['IR_108']['y'] = mi
+      scn['IR_108'].sel(time=np.datetime64('2019-03-01T12:06:13.052000000'))
+
 
 References:
-    - MSG Level 1.5 Image Data Format Description
-    - Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent
-      Spectral Blackbody Radiance
+    - `MSG Level 1.5 Image Format Description`_
+    - `Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance`_
 
+.. _MSG Level 1.5 Image Format Description: http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=
+    PDF_TEN_05105_MSG_IMG_DATA&RevisionSelectionMethod=LatestReleased&Rendition=Web
 
+.. _Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent Spectral Blackbody Radiance:
+    https://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_MSG_SEVIRI_RAD_CALIB&
+    RevisionSelectionMethod=LatestReleased&Rendition=Web
 """
 
 import logging
@@ -47,7 +141,7 @@
                                      annotation_header, base_hdr_map,
                                      image_data_function)
 
-from satpy.readers.seviri_base import SEVIRICalibrationHandler, chebyshev
+from satpy.readers.seviri_base import SEVIRICalibrationHandler, chebyshev, get_cds_time
 from satpy.readers.seviri_base import (CHANNEL_NAMES, VIS_CHANNELS, CALIB, SATNUM)
 
 from satpy.readers.seviri_l1b_native_hdr import (hrit_prologue, hrit_epilogue,
@@ -529,6 +623,10 @@ def get_dataset(self, key, info):
         res.attrs['navigation'] = self.mda['navigation_parameters'].copy()
         res.attrs['georef_offset_corrected'] = self.mda['offset_corrected']
 
+        # Add scanline timestamps as additional y-coordinate
+        res['acq_time'] = ('y', self._get_timestamps())
+        res['acq_time'].attrs['long_name'] = 'Mean scanline acquisition time'
+
         return res
 
     def calibrate(self, data, calibration):
@@ -577,6 +675,11 @@ def calibrate(self, data, calibration):
         logger.debug("Calibration time " + str(datetime.now() - tic))
         return res
 
+    def _get_timestamps(self):
+        """Read scanline timestamps from the segment header"""
+        tline = self.mda['image_segment_line_quality']['line_mean_acquisition']
+        return get_cds_time(days=tline['days'], msecs=tline['milliseconds'])
+
 
 def show(data, negate=False):
     """Show the stretched data.

satpy/tests/reader_tests/test_seviri_base.py

@@ -25,47 +25,56 @@
 
 import sys
 import numpy as np
-from satpy.readers.seviri_base import dec10216, chebyshev
+from satpy.readers.seviri_base import dec10216, chebyshev, get_cds_time
 
 if sys.version_info < (2, 7):
     import unittest2 as unittest
 else:
     import unittest
 
 
-class TestDec10216(unittest.TestCase):
-    """Test the dec10216 function."""
+def chebyshev4(c, x, domain):
+    """Evaluate 4th order Chebyshev polynomial"""
+    start_x, end_x = domain
+    t = (x - 0.5 * (end_x + start_x)) / (0.5 * (end_x - start_x))
+    return c[0] + c[1]*t + c[2]*(2*t**2 - 1) + c[3]*(4*t**3 - 3*t) - 0.5*c[0]
 
+
+class SeviriBaseTest(unittest.TestCase):
     def test_dec10216(self):
+        """Test the dec10216 function."""
         res = dec10216(np.array([255, 255, 255, 255, 255], dtype=np.uint8))
         exp = (np.ones((4, )) * 1023).astype(np.uint16)
         self.assertTrue(np.all(res == exp))
         res = dec10216(np.array([1, 1, 1, 1, 1], dtype=np.uint8))
         exp = np.array([4,  16,  64, 257], dtype=np.uint16)
         self.assertTrue(np.all(res == exp))
 
-
-class TestChebyshev(unittest.TestCase):
-    def chebyshev4(self, c, x, domain):
-        """Evaluate 4th order Chebyshev polynomial"""
-        start_x, end_x = domain
-        t = (x - 0.5 * (end_x + start_x)) / (0.5 * (end_x - start_x))
-        return c[0] + c[1]*t + c[2]*(2*t**2 - 1) + c[3]*(4*t**3 - 3*t) - 0.5*c[0]
-
     def test_chebyshev(self):
         coefs = [1, 2, 3, 4]
         time = 123
         domain = [120, 130]
         res = chebyshev(coefs=[1, 2, 3, 4], time=time, domain=domain)
-        exp = self.chebyshev4(coefs, time, domain)
+        exp = chebyshev4(coefs, time, domain)
         self.assertTrue(np.allclose(res, exp))
 
+    def test_get_cds_time(self):
+        # Scalar
+        self.assertEqual(get_cds_time(days=21246, msecs=12*3600*1000),
+                         np.datetime64('2016-03-03 12:00'))
+
+        # Array
+        days = np.array([21246, 21247, 21248])
+        msecs = np.array([12*3600*1000, 13*3600*1000 + 1, 14*3600*1000 + 2])
+        expected = np.array([np.datetime64('2016-03-03 12:00:00.000'),
+                             np.datetime64('2016-03-04 13:00:00.001'),
+                             np.datetime64('2016-03-05 14:00:00.002')])
+        self.assertTrue(np.all(get_cds_time(days=days, msecs=msecs) == expected))
+
 
 def suite():
     """The test suite for test_scene."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
-    tests = [TestDec10216, TestChebyshev]
-    for test in tests:
-        mysuite.addTest(loader.loadTestsFromTestCase(test))
+    mysuite.addTest(loader.loadTestsFromTestCase(SeviriBaseTest))
     return mysuite

satpy/tests/reader_tests/test_seviri_l1b_hrit.py

@@ -112,6 +112,11 @@ def setUp(self, fromfile):
                     self.reader.mda['navigation_parameters']['satellite_actual_latitude'] = -0.5
                     self.reader.mda['navigation_parameters']['satellite_actual_altitude'] = 35783328
 
+                    tline = np.zeros(nlines, dtype=[('days', '>u2'), ('milliseconds', '>u4')])
+                    tline['days'][1:-1] = 21246 * np.ones(nlines-2)  # 2016-03-03
+                    tline['milliseconds'][1:-1] = np.arange(nlines-2)
+                    self.reader.mda['image_segment_line_quality'] = {'line_mean_acquisition': tline}
+
     def test_get_xy_from_linecol(self):
         """Test get_xy_from_linecol."""
         x__, y__ = self.reader.get_xy_from_linecol(0, 0, (10, 10), (5, 5))
@@ -239,13 +244,17 @@ def get_header_patched(self):
         self.assertRaises(ValueError, HRITMSGFileHandler, filename=None, filename_info=None,
                           filetype_info=None, prologue=pro, epilogue=epi, calib_mode='invalid')
 
+    @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler._get_timestamps')
     @mock.patch('satpy.readers.seviri_l1b_hrit.HRITFileHandler.get_dataset')
     @mock.patch('satpy.readers.seviri_l1b_hrit.HRITMSGFileHandler.calibrate')
-    def test_get_dataset(self, calibrate, parent_get_dataset):
+    def test_get_dataset(self, calibrate, parent_get_dataset, _get_timestamps):
         key = mock.MagicMock(calibration='calibration')
         info = {'units': 'units', 'wavelength': 'wavelength', 'standard_name': 'standard_name'}
+        timestamps = np.array([1, 2, 3], dtype='datetime64[ns]')
+
         parent_get_dataset.return_value = mock.MagicMock()
-        calibrate.return_value = mock.MagicMock(attrs={})
+        calibrate.return_value = xr.DataArray(data=np.zeros((3, 3)), dims=('y', 'x'))
+        _get_timestamps.return_value = timestamps
 
         res = self.reader.get_dataset(key, info)
 
@@ -267,9 +276,29 @@ def test_get_dataset(self, calibrate, parent_get_dataset):
             'navigation': self.reader.mda['navigation_parameters'],
             'georef_offset_corrected': self.reader.mda['offset_corrected']
         })
-
         self.assertDictEqual(attrs_exp, res.attrs)
 
+        # Test timestamps
+        self.assertTrue(np.all(res['acq_time'] == timestamps))
+        self.assertEqual(res['acq_time'].attrs['long_name'], 'Mean scanline acquisition time')
+
+    def test_get_timestamps(self):
+        tline = self.reader._get_timestamps()
+
+        # First and last scanline have invalid timestamps (space)
+        self.assertTrue(np.isnat(tline[0]))
+        self.assertTrue(np.isnat(tline[-1]))
+
+        # Test remaining lines
+        year = tline.astype('datetime64[Y]').astype(int) + 1970
+        month = tline.astype('datetime64[M]').astype(int) % 12 + 1
+        day = (tline.astype('datetime64[D]') - tline.astype('datetime64[M]') + 1).astype(int)
+        msec = (tline - tline.astype('datetime64[D]')).astype(int)
+        self.assertTrue(np.all(year[1:-1] == 2016))
+        self.assertTrue(np.all(month[1:-1] == 3))
+        self.assertTrue(np.all(day[1:-1] == 3))
+        self.assertTrue(np.all(msec[1:-1] == np.arange(len(tline) - 2)))
+
 
 class TestHRITMSGPrologueFileHandler(unittest.TestCase):
     """Test the HRIT prologue file handler."""