Merge pull request #473 from djhoese/bugfix-ahi-hsd-time

Change combine_metadata to average any 'time' fields

satpy/dataset.py

@@ -24,9 +24,11 @@
 """Dataset objects.
 """
 
+import sys
 import logging
 import numbers
 from collections import namedtuple
+from datetime import datetime
 
 import numpy as np
 
@@ -46,20 +48,54 @@ def id(self):
         return DatasetID.from_dict(self.attrs)
 
 
-def combine_metadata(*metadata_objects):
+def average_datetimes(dt_list):
+    """Average a series of datetime objects.
+
+    .. note::
+
+        This function assumes all datetime objects are naive and in the same
+        time zone (UTC).
+
+    Args:
+        dt_list (iterable): Datetime objects to average
+
+    Returns: Average datetime as a datetime object
+
+    """
+    if sys.version_info < (3, 3):
+        # timestamp added in python 3.3
+        import time
+
+        def timestamp_func(dt):
+            return time.mktime(dt.timetuple())
+    else:
+        timestamp_func = datetime.timestamp
+
+    total = [timestamp_func(dt) for dt in dt_list]
+    return datetime.fromtimestamp(sum(total) / len(total))
+
+
+def combine_metadata(*metadata_objects, **kwargs):
     """Combine the metadata of two or more Datasets.
 
+    If any keys are not equal or do not exist in all provided dictionaries
+    then they are not included in the returned dictionary.
+    By default any keys with the word 'time' in them and consisting
+    of datetime objects will be averaged. This is to handle cases where
+    data were observed at almost the same time but not exactly.
+
     Args:
         *metadata_objects: MetadataObject or dict objects to combine
+        average_times (bool): Average any keys with 'time' in the name
 
     Returns:
-        the combined metadata
+        dict: the combined metadata
 
     """
+    average_times = kwargs.get('average_times', True)  # python 2 compatibility (no kwarg after *args)
     shared_keys = None
     info_dicts = []
-    # grab all of the dictionary objects provided and make a set of the shared
-    # keys
+    # grab all of the dictionary objects provided and make a set of the shared keys
     for metadata_object in metadata_objects:
         if isinstance(metadata_object, dict):
             metadata_dict = metadata_object
@@ -82,6 +118,8 @@ def combine_metadata(*metadata_objects):
         if any_arrays:
             if all(np.all(val == values[0]) for val in values[1:]):
                 shared_info[k] = values[0]
+        elif 'time' in k and isinstance(values[0], datetime) and average_times:
+            shared_info[k] = average_datetimes(values)
         elif all(val == values[0] for val in values[1:]):
             shared_info[k] = values[0]
 

satpy/readers/ahi_hsd.py

@@ -1,31 +1,41 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
+#
 # Copyright (c) 2014-2018 PyTroll developers
-
+#
 # Author(s):
-
+#
 #   Adam.Dybbroe <adam.dybbroe@smhi.se>
 #   Cooke, Michael.C, UK Met Office
 #   Martin Raspaud <martin.raspaud@smhi.se>
-
+#
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
-
+#
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
-
+#
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-"""Advanced Himawari Imager (AHI) standard format data reader
+"""Advanced Himawari Imager (AHI) standard format data reader. The HSD format
+that this reader reads are described at the URL below:
 
 http://www.data.jma.go.jp/mscweb/en/himawari89/space_segment/spsg_ahi.html
 
+Time Information
+****************
+
+AHI observations use the idea of a "scheduled" time and an "observation time.
+The "scheduled" time is when the instrument was told to record the data,
+usually at a specific and consistent interval. The "observation" time is when
+the data was actually observed. Scheduled time can be accessed from the
+`scheduled_time` metadata key and observation time from the `start_time` key.
+
 """
 
 import logging
@@ -44,11 +54,6 @@
                      "6", "7", "8", "9", "10",
                      "11", "12", "13", "14", "15", "16")
 
-
-class CalibrationError(ValueError):
-    pass
-
-
 logger = logging.getLogger('ahi_hsd')
 
 # Basic information block:
@@ -215,9 +220,7 @@ class CalibrationError(ValueError):
 
 
 class AHIHSDFileHandler(BaseFileHandler):
-
-    """AHI standard format reader
-    """
+    """AHI standard format reader."""
 
     def __init__(self, filename, filename_info, filetype_info):
         """Initialize the reader."""
@@ -250,18 +253,19 @@ def __init__(self, filename, filename_info, filetype_info):
         self.platform_name = np2str(self.basic_info['satellite'])
         self.sensor = 'ahi'
 
-    def get_shape(self, dsid, ds_info):
-        return int(self.data_info['number_of_lines']), int(self.data_info['number_of_columns'])
-
     @property
     def start_time(self):
-        return (datetime(1858, 11, 17) +
-                timedelta(days=float(self.basic_info['observation_start_time'])))
+        return datetime(1858, 11, 17) + timedelta(days=float(self.basic_info['observation_start_time']))
 
     @property
     def end_time(self):
-        return (datetime(1858, 11, 17) +
-                timedelta(days=float(self.basic_info['observation_end_time'])))
+        return datetime(1858, 11, 17) + timedelta(days=float(self.basic_info['observation_end_time']))
+
+    @property
+    def scheduled_time(self):
+        """Time this band was scheduled to be recorded."""
+        timeline = "{:04d}".format(self.basic_info['observation_timeline'][0])
+        return self.start_time.replace(hour=int(timeline[:2]), minute=int(timeline[2:4]), second=0, microsecond=0)
 
     def get_dataset(self, key, info):
         return self.read_band(key, info)
@@ -311,10 +315,6 @@ def get_area_def(self, dsid):
         self.area = area
         return area
 
-    def get_lonlats(self, key, info, lon_out, lat_out):
-        logger.debug('Computing area for %s', str(key))
-        lon_out[:], lat_out[:] = self.area.get_lonlats()
-
     def geo_mask(self):
         """Masking the space pixels from geometry info."""
         cfac = np.uint32(self.proj_info['CFAC'])
@@ -346,7 +346,7 @@ def ellipse(line, col):
         return ellipse(lines_idx[:, None], cols_idx[None, :])
 
     def read_band(self, key, info):
-        """Read the data"""
+        """Read the data."""
         tic = datetime.now()
         header = {}
         with open(self.filename, "rb") as fp_:
@@ -428,19 +428,17 @@ def read_band(self, key, info):
                                           dtype='<u2',  shape=(nlines, ncols), mode='r'),
                                 chunks=CHUNK_SIZE)
         res = da.where(res == 65535, np.float32(np.nan), res)
-
         self._header = header
 
         logger.debug("Reading time " + str(datetime.now() - tic))
-
         res = self.calibrate(res, key.calibration)
-
         new_info = dict(units=info['units'],
                         standard_name=info['standard_name'],
                         wavelength=info['wavelength'],
                         resolution='resolution',
                         id=key,
                         name=key.name,
+                        scheduled_time=self.scheduled_time,
                         platform_name=self.platform_name,
                         sensor=self.sensor,
                         satellite_longitude=float(
@@ -473,23 +471,20 @@ def calibrate(self, data, calibration):
         return data
 
     def convert_to_radiance(self, data):
-        """Calibrate to radiance.
-        """
+        """Calibrate to radiance."""
 
         gain = self._header["block5"]["gain_count2rad_conversion"][0]
         offset = self._header["block5"]["offset_count2rad_conversion"][0]
 
         return data * gain + offset
 
     def _vis_calibrate(self, data):
-        """Visible channel calibration only.
-        """
+        """Visible channel calibration only."""
         coeff = self._header["calibration"]["coeff_rad2albedo_conversion"]
         return (data * coeff * 100).clip(0)
 
     def _ir_calibrate(self, data):
-        """IR calibration
-        """
+        """IR calibration."""
 
         cwl = self._header['block5']["central_wave_length"][0] * 1e-6
         c__ = self._header['calibration']["speed_of_light"][0]

satpy/tests/reader_tests/test_ahi_hsd.py

@@ -27,6 +27,8 @@
 except ImportError:
     import mock
 
+import numpy as np
+from datetime import datetime
 from satpy.readers.ahi_hsd import AHIHSDFileHandler
 
 
@@ -118,14 +120,61 @@ def test_segment(self, fromfile, np2str):
                                                     5500000.035542117, -2200000.0142168473))
 
 
-# pprint.pprint(dict(zip((item[0] for item in self.proj_info.dtype.descr), self.proj_info)))
+class TestAHIHSDFileHandler(unittest.TestCase):
+    @mock.patch('satpy.readers.ahi_hsd.np2str')
+    @mock.patch('satpy.readers.ahi_hsd.np.fromfile')
+    def setUp(self, fromfile, np2str):
+        """Create a test file handler."""
+        np2str.side_effect = lambda x: x
+        m = mock.mock_open()
+        with mock.patch('satpy.readers.ahi_hsd.open', m, create=True):
+            fh = AHIHSDFileHandler(None, {'segment_number': 8, 'total_segments': 10}, None)
+            fh.proj_info = {'CFAC': 40932549,
+                            'COFF': 5500.5,
+                            'LFAC': 40932549,
+                            'LOFF': 5500.5,
+                            'blocklength': 127,
+                            'coeff_for_sd': 1737122264.0,
+                            'distance_from_earth_center': 42164.0,
+                            'earth_equatorial_radius': 6378.137,
+                            'earth_polar_radius': 6356.7523,
+                            'hblock_number': 3,
+                            'req2_rpol2': 1.006739501,
+                            'req2_rpol2_req2': 0.0066943844,
+                            'resampling_size': 4,
+                            'resampling_types': 0,
+                            'rpol2_req2': 0.993305616,
+                            'spare': '',
+                            'sub_lon': 140.7}
+
+            fh.data_info = {'blocklength': 50,
+                            'compression_flag_for_data': 0,
+                            'hblock_number': 2,
+                            'number_of_bits_per_pixel': 16,
+                            'number_of_columns': 11000,
+                            'number_of_lines': 1100,
+                            'spare': ''}
+            fh.basic_info = {
+                'observation_start_time': np.array([58413.12523839]),
+                'observation_end_time': np.array([58413.12562439]),
+                'observation_timeline': np.array([300]),
+            }
+
+            self.fh = fh
+
+    def test_time_properties(self):
+        """Test start/end/scheduled time properties."""
+        self.assertEqual(self.fh.start_time, datetime(2018, 10, 22, 3, 0, 20, 596896))
+        self.assertEqual(self.fh.end_time, datetime(2018, 10, 22, 3, 0, 53, 947296))
+        self.assertEqual(self.fh.scheduled_time, datetime(2018, 10, 22, 3, 0, 0, 0))
 
 
 def suite():
     """The test suite for test_scene."""
     loader = unittest.TestLoader()
     mysuite = unittest.TestSuite()
     mysuite.addTest(loader.loadTestsFromTestCase(TestAHIHSDNavigation))
+    mysuite.addTest(loader.loadTestsFromTestCase(TestAHIHSDFileHandler))
     return mysuite