Initial commit for H-SAF products in NC format

satpy/etc/composites/hsaf.yaml

@@ -1,4 +1,4 @@
-sensor_name: hsaf
+sensor_name: visir/hsaf
 
 composites:
   instantaneous_rainrate_3:
@@ -21,3 +21,9 @@ composites:
     prerequisites:
     - name: h05B
     standard_name: accum_rainrate_5b
+
+  instantaneous_rainfall_rate:
+    compositor: !!python/name:satpy.composites.GenericCompositor
+    prerequisites:
+    - name: h63
+    standard_name: instantaneous_rainfall_rate

satpy/etc/enhancements/generic.yaml

@@ -1253,3 +1253,16 @@ enhancements:
          stretch: crude
          min_stretch: [0,0,0]
          max_stretch: [1,1,1]
+
+  instantaneous_rainfall_rate:
+    standard_name: instantaneous_rainfall_rate
+    operations:
+      - name: colorize
+        method: !!python/name:satpy.enhancements.colorize
+        kwargs:
+          palettes:
+             - {colors: 'spectral',
+                reverse: true,
+                min_value: 0.0,
+                max_value: 65.0,
+               }

satpy/etc/enhancements/hsaf.yaml

@@ -0,0 +1,13 @@
+enhancements:
+  instantaneous_rainfall_rate:
+    standard_name: instantaneous_rainfall_rate
+    operations:
+      - name: colorize
+        method: !!python/name:satpy.enhancements.colorize
+        kwargs:
+          palettes:
+             - {colors: 'spectral',
+                reverse: true,
+                min_value: 0.0,
+                max_value: 65.0,
+               }

satpy/etc/readers/hsaf_nc.yaml

@@ -0,0 +1,34 @@
+reader:
+  name: hsaf_nc
+  short_name: Hydrology SAF
+  long_name: Hydrology SAF products in NetCDF format
+  description: Reader for Hydrology SAF products
+  status: Beta, only h63 currently supported
+  supports_fsspec: false
+  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
+  sensors: [hsaf]
+
+file_types:
+  hsafnc_63:
+    file_reader: !!python/name:satpy.readers.hsaf_nc.HSAFNCFileHandler
+    file_patterns: ['h63_{start_time:%Y%m%d_%H%M}_{region:s}.nc']
+  hsafnc_90:
+    file_reader: !!python/name:satpy.readers.hsaf_nc.HSAFNCFileHandler
+    file_patterns: ['h90_{start_time:%Y%m%d_%H%M}_{acc_time:2s}_{region:s}.nc']
+datasets:
+  h63:
+    name: h63
+    file_key: rr
+    sensor: hsaf
+    resolution: 3000
+    standard_name: instantaneous_rainfall_rate
+    units: mm/h
+    file_type: hsafnc_63
+  h90:
+    name: h90
+    file_key: acc_rr
+    sensor: hsaf
+    resolution: 3000
+    standard_name: accumulated_precipitation_estimate
+    units: mm
+    file_type: hsafnc_90

satpy/readers/hsaf_nc.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2024.
+#
+# 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 NetCDF Hydrology SAF products. In this beta version, only H63 is supported."""
+import logging
+import os
+from contextlib import suppress
+from datetime import timedelta
+from pyresample import geometry
+
+import dask.array as da
+import numpy as np
+import xarray as xr
+from satpy.readers._geos_area import get_area_definition, get_area_extent
+
+from satpy.readers.file_handlers import BaseFileHandler
+from satpy.readers.utils import unzip_file
+from satpy.utils import get_chunk_size_limit
+
+logger = logging.getLogger(__name__)
+
+CHUNK_SIZE = get_chunk_size_limit()
+
+PLATFORM_NAMES = {"MSG1": "Meteosat-8",
+                  "MSG2": "Meteosat-9",
+                  "MSG3": "Meteosat-10",
+                  "MSG4": "Meteosat-11",
+                  "GOES16": "GOES-16",
+                  "GOES17": "GOES-17",
+                  }
+
+class HSAFNCFileHandler(BaseFileHandler):
+    """NWCSAF PPS&MSG NetCDF reader."""
+
+    def __init__(self, filename, filename_info, filetype_info):
+        """Init method."""
+        super(HSAFNCFileHandler, self).__init__(filename, filename_info,
+                                                filetype_info)
+
+        self._unzipped = unzip_file(self.filename)
+        if self._unzipped:
+            self.filename = self._unzipped
+
+        self.cache = {}
+        self.nc = xr.open_dataset(self.filename,
+                                  decode_cf=True,
+                                  mask_and_scale=False,
+                                  chunks=CHUNK_SIZE)
+
+        if 'xc' in self.nc.dims:
+            self.nc = self.nc.rename({"xc": "x", "yc": "y"})
+        elif 'nx' in self.nc.dims:
+            self.nc = self.nc.rename({"nx": "x", "ny": "y"})
+
+        try:
+            kwrgs = {"sat_id": self.nc.attrs["satellite_identifier"]}
+        except KeyError:
+            kwrgs = {"sat_id": None}
+
+        self.set_platform_and_sensor(**kwrgs)
+
+    def __del__(self):
+        """Delete the instance."""
+        if self._unzipped:
+            try:
+                os.remove(self._unzipped)
+            except OSError:
+                pass
+
+    def set_platform_and_sensor(self, **kwargs):
+        """Set some metadata: platform_name, sensors, and pps (identifying PPS or Geo)."""
+        self.platform_name = PLATFORM_NAMES.get(kwargs["sat_id"], 'N/A')
+
+        self.sensor = "seviri"
+
+    def get_dataset(self, dsid, info):
+            """Load a dataset."""
+            dsid_name = info["file_key"]
+            if dsid_name in self.cache:
+                logger.debug("Get the data set from cache: %s.", dsid_name)
+                return self.cache[dsid_name]
+
+            logger.debug("Reading %s.", dsid_name)
+            variable = self.nc[dsid_name]
+
+            # Data is transposed in file, fix it here
+            variable.data = variable.data.T
+
+            variable.attrs["start_time"] = self.start_time
+            variable.attrs["end_time"] = self.end_time
+
+            # Fill value is not defined as an attribute, manually specify
+            variable.data = da.where(variable.data > 0, variable.data, np.nan)
+
+            return variable
+
+    def get_area_def(self, dsid):
+        """Get the area definition of the band."""
+
+        val_dict = {}
+        for keyval in self.nc.attrs['cgms_projection'].split():
+            try:
+                key, val = keyval.split('=')
+                val_dict[key] = val
+            except:
+                pass
+
+        pdict = {}
+        pdict['scandir'] = 'N2S'
+        pdict["ssp_lon"] = np.float32(self.nc.attrs["sub-satellite_longitude"][:-1])
+        pdict['a'] = float(val_dict['+r_eq'])*1000
+        pdict['b'] = float(val_dict['+r_pol'])*1000
+        pdict['h'] = float(val_dict['+h'])*1000 - pdict['a']
+        pdict['loff'] = float(val_dict['+loff'])
+        pdict['coff'] = float(val_dict['+coff'])
+        pdict['lfac'] = -float(val_dict['+lfac'])
+        pdict['cfac'] = -float(val_dict['+cfac'])
+        pdict['ncols'] = self.nc.x.size
+        pdict['nlines'] = self.nc.y.size
+        pdict["a_name"] = "seviri_geos_fds"
+        pdict["a_desc"] = "SEVIRI full disk area at native resolution"
+        pdict["p_id"] = "seviri_fixed_grid"
+
+        area_extent = get_area_extent(pdict)
+        fg_area_def = get_area_definition(pdict, area_extent)
+        return fg_area_def
+
+    @property
+    def start_time(self):
+        """Return the start time of the object."""
+        return self.filename_info["start_time"]
+
+    @property
+    def end_time(self):
+        """Return the end time of the object.
+        
+        The product does not provide the end time, so the start time is used."""
+        return self.filename_info["start_time"]