Reader for Hydrology SAF precipitation products (#814)

Reader for Hydrology SAF precipitation products

doc/source/index.rst

@@ -205,6 +205,10 @@ the base Satpy installation.
     * - TROPOMI L2 data in NetCDF4 format
       - `tropomi_l2`
       - Beta
+    * - Hydrology SAF products in GRIB format
+      - `hsaf_grib`
+      - | Beta
+        | Only the h03, h03b, h05 and h05B products are supported at-present
 
 
 Indices and tables

satpy/etc/composites/hsaf.yaml

@@ -0,0 +1,23 @@
+sensor_name: hsaf
+
+composites:
+  instantaneous_rainrate_3:
+    compositor: !!python/name:satpy.composites.GenericCompositor
+    prerequisites:
+    - name: h03
+    standard_name: instantaneous_rainrate_3
+  instantaneous_rainrate_3b:
+    compositor: !!python/name:satpy.composites.GenericCompositor
+    prerequisites:
+    - name: h03B
+    standard_name: instantaneous_rainrate_3b
+  accum_rainrate_5:
+    compositor: !!python/name:satpy.composites.GenericCompositor
+    prerequisites:
+    - name: h05
+    standard_name: accum_rainrate_5
+  accum_rainrate_5b:
+    compositor: !!python/name:satpy.composites.GenericCompositor
+    prerequisites:
+    - name: h05B
+    standard_name: accum_rainrate_5b

satpy/etc/readers/hsaf_grib.yaml

@@ -0,0 +1,47 @@
+reader:
+  description: Reader for Hydrology SAF products
+  name: hsaf_grib
+  reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader
+  sensors: [hsaf]
+
+file_types:
+  hsafgrib:
+    file_reader: !!python/name:satpy.readers.hsaf_grib.HSAFFileHandler
+    file_patterns: ['h03_{sensing_time:%Y%m%d_%H%M}_{region:3s}.grb',
+                    'h05_{sensing_time:%Y%m%d_%H%M}_{accum_time:2s}_{region:3s}.grb',
+                    'h03B_{sensing_time:%Y%m%d_%H%M}_{region:3s}.grb',
+                    'h05B_{sensing_time:%Y%m%d_%H%M}_{accum_time:2s}_{region:3s}.grb']
+
+datasets:
+  h03:
+    name: h03
+    msg_name: irrate
+    sensor: hsaf
+    resolution: 3000
+    standard_name: instantaneous_rainfall_rate
+    units: kg m-2 s-1
+    file_type: hsafgrib
+  h03B:
+    name: h03B
+    msg_name: irrate
+    sensor: hsaf
+    resolution: 3000
+    standard_name: instantaneous_rainfall_rate
+    units: kg m-2 s-1
+    file_type: hsafgrib
+  h05:
+    name: h05
+    msg_name: accumrain
+    sensor: hsaf
+    resolution: 3000
+    standard_name: accumulated_rainfall_rate
+    units: kg m-2
+    file_type: hsafgrib
+  h05B:
+    name: h05B
+    msg_name: accumrain
+    sensor: hsaf
+    resolution: 3000
+    standard_name: accumulated_rainfall_rate
+    units: kg m-2
+    file_type: hsafgrib

satpy/readers/hsaf_grib.py

@@ -0,0 +1,171 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2019.
+#
+# This file is part of satpy.
+#
+# satpy 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.
+#
+# satpy 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
+# satpy.  If not, see <http://www.gnu.org/licenses/>.
+"""A reader for files produced by the Hydrology SAF
+
+Currently this reader depends on the `pygrib` python package. The `eccodes`
+package from ECMWF is preferred, but does not support python 3 at the time
+of writing.
+
+"""
+import logging
+import numpy as np
+import xarray as xr
+import dask.array as da
+from pyresample import geometry
+from datetime import datetime, timedelta
+
+from satpy import CHUNK_SIZE
+from satpy.readers.file_handlers import BaseFileHandler
+import pygrib
+
+LOG = logging.getLogger(__name__)
+
+
+CF_UNITS = {
+    'none': '1',
+}
+
+
+class HSAFFileHandler(BaseFileHandler):
+
+    def __init__(self, filename, filename_info, filetype_info):
+        super(HSAFFileHandler, self).__init__(filename,
+                                              filename_info,
+                                              filetype_info)
+
+        self._msg_datasets = {}
+        self._start_time = None
+        self._end_time = None
+        try:
+            with pygrib.open(self.filename) as grib_file:
+                first_msg = grib_file.message(1)
+                analysis_time = self._get_datetime(first_msg)
+                self._analysis_time = analysis_time
+                self.metadata = self.get_metadata(first_msg)
+
+        except (RuntimeError, KeyError):
+            raise IOError("Unknown GRIB file format: {}".format(self.filename))
+
+    @staticmethod
+    def _get_datetime(msg):
+        dtstr = str(msg['dataDate']) + str(msg['dataTime']).zfill(4)
+        return datetime.strptime(dtstr, "%Y%m%d%H%M")
+
+    @property
+    def analysis_time(self):
+        """
+        Get validity time of this file
+        """
+        return self._analysis_time
+
+    def get_metadata(self, msg):
+        try:
+            center_description = msg['centreDescription']
+        except (RuntimeError, KeyError):
+            center_description = None
+        ds_info = {
+            'filename': self.filename,
+            'shortName': msg['shortName'],
+            'long_name': msg['name'],
+            'units': msg['units'],
+            'centreDescription': center_description,
+            'data_time': self._analysis_time,
+            'nx': msg['Nx'],
+            'ny': msg['Ny'],
+            'projparams': msg.projparams
+        }
+        return ds_info
+
+    def get_area_def(self, dsid):
+        """
+        Get area definition for message.
+        """
+        msg = self._get_message(1)
+        try:
+            return self._get_area_def(msg)
+        except (RuntimeError, KeyError):
+            raise RuntimeError("Unknown GRIB projection information")
+
+    def _get_area_def(self, msg):
+        """
+        Get the area definition of the datasets in the file.
+        """
+
+        proj_param = msg.projparams.copy()
+
+        Rx = 2 * np.arcsin(1. / msg['NrInRadiusOfEarth']) / msg['dx']
+        Ry = 2 * np.arcsin(1. / msg['NrInRadiusOfEarth']) / msg['dy']
+
+        x_0 = - msg['XpInGridLengths']
+        x_1 = msg['Nx'] - msg['XpInGridLengths']
+        y_0 = (msg['Ny'] - msg['YpInGridLengths']) * -1
+        y_1 = msg['YpInGridLengths']
+
+        min_x = (x_0 * Rx) * proj_param['h']
+        max_x = (x_1 * Rx) * proj_param['h']
+
+        min_y = (y_0 * Ry) * proj_param['h']
+        max_y = (y_1 * Ry) * proj_param['h']
+
+        area_extent = (min_x, min_y, max_x, max_y)
+
+        area = geometry.AreaDefinition('hsaf_region',
+                                       'A region from H-SAF',
+                                       'geos',
+                                       proj_param,
+                                       msg['Nx'],
+                                       msg['Ny'],
+                                       area_extent)
+
+        return area
+
+    def _get_message(self, idx):
+        with pygrib.open(self.filename) as grib_file:
+            msg = grib_file.message(idx)
+            return msg
+
+    def get_dataset(self, ds_id, ds_info):
+        """Read a GRIB message into an xarray DataArray."""
+        if (ds_id.name not in self.filename):
+            raise IOError("File does not contain {} data".format(ds_id.name))
+
+        msg = self._get_message(1)
+
+        ds_info = self.get_metadata(msg)
+        ds_info['end_time'] = ds_info['data_time']
+
+        if (ds_id.name == 'h05' or ds_id.name == 'h05B'):
+            flen = len(self.filename)
+            timedelt = self.filename[flen-10:flen-8]
+            ds_info['start_time'] = (ds_info['end_time'] -
+                                     timedelta(hours=int(timedelt)))
+        else:
+            ds_info['start_time'] = ds_info['end_time']
+        fill = msg['missingValue']
+        data = msg.values.astype(np.float32)
+        if msg.valid_key('jScansPositively') and msg['jScansPositively'] == 1:
+            data = data[::-1]
+
+        if isinstance(data, np.ma.MaskedArray):
+            data = data.filled(np.nan)
+            data = da.from_array(data, chunks=CHUNK_SIZE)
+        else:
+            data[data == fill] = np.nan
+            data = da.from_array(data, chunks=CHUNK_SIZE)
+
+        return xr.DataArray(data, attrs=ds_info, dims=('y', 'x'))

satpy/tests/reader_tests/__init__.py

@@ -38,7 +38,8 @@
                                       test_hdfeos_base, test_modis_l2,
                                       test_electrol_hrit, test_mersi2_l1b,
                                       test_avhrr_l1b_gaclac, test_vaisala_gld360,
-                                      test_fci_l1c_fdhsi, test_tropomi_l2)
+                                      test_fci_l1c_fdhsi, test_tropomi_l2,
+                                      test_hsaf_grib)
 
 
 if sys.version_info < (2, 7):
@@ -92,5 +93,6 @@ def suite():
     mysuite.addTests(test_vaisala_gld360.suite())
     mysuite.addTests(test_fci_l1c_fdhsi.suite())
     mysuite.addTests(test_tropomi_l2.suite())
+    mysuite.addTests(test_hsaf_grib.suite())
 
     return mysuite