Add scanline acquisition time to hrit_jma

satpy/readers/hrit_jma.py

@@ -107,6 +107,15 @@
 }
 
 
+def mjd2datetime64(mjd):
+    """Convert Modified Julian Day (MJD) to datetime64."""
+
+    epoch = np.datetime64('1858-11-17 00:00')
+    day2usec = 24 * 3600 * 1E6
+    mjd_usec = (mjd * day2usec).astype(int).astype('timedelta64[us]')
+    return epoch + mjd_usec
+
+
 class HRITJMAFileHandler(HRITFileHandler):
     """JMA HRIT format reader."""
 
@@ -149,6 +158,7 @@ def __init__(self, filename, filename_info, filetype_info):
         if self.area_id not in AREA_NAMES:
             self.area_id = UNKNOWN_AREA
         self.area = self._get_area_def()
+        self.acq_time = self._get_acq_time()
 
     def _get_platform(self):
         """Get the platform name.
@@ -265,6 +275,10 @@ def get_dataset(self, key, info):
         # Calibrate and mask space pixels
         res = self._mask_space(self.calibrate(res, key.calibration))
 
+        # Add scanline acquisition time
+        res['acq_time'] = ('y', self.acq_time)
+        res['acq_time'].attrs['long_name'] = 'Scanline acquisition time'
+
         # Update attributes
         res.attrs.update(info)
         res.attrs['platform_name'] = self.platform
@@ -283,6 +297,42 @@ def _mask_space(self, data):
         geomask = get_geostationary_mask(area=self.area)
         return data.where(geomask)
 
+    def _get_acq_time(self):
+        """
+        Acquisition times for a subset of scanlines are stored in the header
+        as follows:
+
+        b'LINE:=1\rTIME:=54365.022558\rLINE:=21\rTIME:=54365.022664\r...'
+
+        Missing timestamps in between are computed using linear interpolation.
+        """
+        buf_b = np.frombuffer(self.mda['image_observation_time'],
+                              dtype=image_observation_time)
+
+        # Replace \r by \n before encoding, otherwise encoding will drop all
+        # elements except the last one
+        buf_s = b''.join(buf_b['times']).replace(b'\r', b'\n').decode()
+
+        # Split into key:=value pairs; then extract line number and timestamp
+        splits = buf_s.strip().split('\n')
+        lines_sparse = [int(s.split(':=')[1]) for s in splits[0::2]]
+        times_sparse = [float(s.split(':=')[1]) for s in splits[1::2]]
+
+        if self.platform == HIMAWARI8:
+            # Only 3 timestamps, and only the first and last are usable
+            # (the second equals the third).
+            lines_sparse = [lines_sparse[0], lines_sparse[-1]]
+            times_sparse = [times_sparse[0], times_sparse[-1]]
+
+        # Compute missing timestamps using linear interpolation.
+        lines = np.arange(lines_sparse[0], lines_sparse[-1]+1)
+        times = np.interp(lines, lines_sparse, times_sparse)
+
+        # Convert to np.datetime64
+        times64 = mjd2datetime64(times)
+
+        return times64
+
     @staticmethod
     def _interp(arr, cal):
         return np.interp(arr.ravel(), cal[:, 0], cal[:, 1]).reshape(arr.shape)
@@ -304,3 +354,13 @@ def calibrate(self, data, calibration):
         res = res.where(data < 65535)
         logger.debug("Calibration time " + str(datetime.now() - tic))
         return res
+
+    @property
+    def start_time(self):
+        """Get start time of the scan."""
+        return self.acq_time[0].astype(datetime)
+
+    @property
+    def end_time(self):
+        """Get end time of the scan."""
+        return self.acq_time[-1].astype(datetime)

satpy/tests/reader_tests/test_ahi_hrit.py

@@ -36,19 +36,40 @@ def _get_reader(self, mocked_init, mda, filename_info=None):
         HRITJMAFileHandler.mda = mda
         return HRITJMAFileHandler('filename', filename_info, {})
 
+    def _get_acq_time(self, nlines):
+        """Get sample header entry for scanline acquisition times.
+
+        Lines: 1, 21, 41, 61, ..., nlines
+        Times: 1970-01-01 00:00 + (1, 21, 41, 61, ..., nlines) hours
+
+        So the interpolated times are expected to be 1970-01-01 +
+        (1, 2, 3, 4, ..., nlines) hours. Note that there will be some
+        floating point inaccuracies, because timestamps are stored
+        with only 6 decimals precision.
+        """
+        mjd_1970 = 40587.0
+        lines_sparse = np.array(list(range(1, nlines, 20)) + [nlines])
+        times_sparse = mjd_1970 + lines_sparse / 24 / 3600
+        acq_time_s = ['LINE:={}\rTIME:={:.6f}\r'.format(l, t)
+                      for l, t in zip(lines_sparse, times_sparse)]
+        acq_time_b = ''.join(acq_time_s).encode()
+        return acq_time_b
+
     def _get_mda(self, loff=5500.0, coff=5500.0, nlines=11000, ncols=11000,
-                 segno=0, numseg=1, vis=True):
+                 segno=0, numseg=1, vis=True, platform='Himawari-8'):
         """Create metadata dict like HRITFileHandler would do it."""
         if vis:
             idf = b'$HALFTONE:=16\r_NAME:=VISIBLE\r_UNIT:=ALBEDO(%)\r' \
                   b'0:=-0.10\r1023:=100.00\r65535:=100.00\r'
         else:
             idf = b'$HALFTONE:=16\r_NAME:=INFRARED\r_UNIT:=KELVIN\r' \
                   b'0:=329.98\r1023:=130.02\r65535:=130.02\r'
-
+        proj_h8 = b'GEOS(140.70)                    '
+        proj_mtsat2 = b'GEOS(145.00)                    '
+        proj_name = proj_h8 if platform == 'Himawari-8' else proj_mtsat2
         return {'image_segm_seq_no': segno,
                 'total_no_image_segm': numseg,
-                'projection_name': b'GEOS(140.70)                    ',
+                'projection_name': proj_name,
                 'projection_parameters': {
                     'a': 6378169.00,
                     'b': 6356583.80,
@@ -60,7 +81,8 @@ def _get_mda(self, loff=5500.0, coff=5500.0, nlines=11000, ncols=11000,
                 'loff': loff,
                 'number_of_columns': ncols,
                 'number_of_lines': nlines,
-                'image_data_function': idf}
+                'image_data_function': idf,
+                'image_observation_time': self._get_acq_time(nlines)}
 
     def test_init(self):
         """Test creating the file handler."""
@@ -87,6 +109,9 @@ def test_init(self):
                                  [65535,  100]])
         self.assertTrue(np.all(reader.calibration_table == cal_expected))
 
+        # Check if scanline timestamps are there (dedicated test below)
+        self.assertIsInstance(reader.acq_time, np.ndarray)
+
         # Check platform
         self.assertEqual(reader.platform, HIMAWARI8)
 
@@ -271,3 +296,23 @@ def test_get_dataset(self, base_get_dataset):
         with mock.patch('logging.Logger.error') as log_mock:
             reader.get_dataset(key, {'units': '%', 'sensor': 'jami'})
             log_mock.assert_called()
+
+    def test_mjd2datetime64(self):
+        """Test conversion from modified julian day to datetime64"""
+        from satpy.readers.hrit_jma import mjd2datetime64
+        self.assertEqual(mjd2datetime64(0), np.datetime64('1858-11-17', 'us'))
+        self.assertEqual(mjd2datetime64(40587.5), np.datetime64('1970-01-01 12:00', 'us'))
+
+    def test_get_acq_time(self):
+        """Test computation of scanline acquisition times."""
+        dt_line = np.arange(1, 11000+1).astype('timedelta64[s]')
+        acq_time_exp = np.datetime64('1970-01-01', 'us') + dt_line
+
+        for platform in ['Himawari-8', 'MTSAT-2']:
+            # Results are not exactly identical because timestamps are stored in
+            # the header with only 6 decimals precision (max diff here: 45 msec).
+            mda = self._get_mda(platform=platform)
+            reader = self._get_reader(mda=mda)
+            np.testing.assert_allclose(reader.acq_time.astype(int),
+                                       acq_time_exp.astype(int),
+                                       atol=45000)