Add combined GRIB reader for both SEVIRI and FCI L2 products

AUTHORS.md

@@ -88,3 +88,4 @@ The following people have made contributions to this project:
 - [Xin Zhang (zxdawn)](https://github.com/zxdawn)
 - [Yufei Zhu (yufeizhu600)](https://github.com/yufeizhu600)
 - [Youva Aoun (YouvaEUMex)](https://github.com/YouvaEUMex)
+- [David Navia (dnaviap)](https://github.com/dnaviap)

satpy/etc/readers/fci_l2_grib.yaml

@@ -0,0 +1,28 @@
+reader:
+  name: fci_l2_grib
+  short_name: FCI L2 GRIB2
+  long_name: MTG FCI L2 data in GRIB2 format
+  description: Reader for EUMETSAT MTG FCI L2 files in GRIB2 format.
+  status: Nominal
+  supports_fsspec: false
+  sensors: [fci]
+  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
+
+file_types:
+  grib_fci_clm:
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
+    file_patterns:
+      - '{pflag}_{location_indicator},{data_designator},MTI{spacecraft_id:1d}+FCI-2-CLM-{subtype}-{coverage}-{subsetting}-{component1}-{component2}-{component3}-{purpose}-GRIB2_{oflag}_{originator}_{processing_time:%Y%m%d%H%M%S}_{facility_or_tool}_{environment}_{start_time:%Y%m%d%H%M%S}_{end_time:%Y%m%d%H%M%S}_{processing_mode}_{special_compression}_{disposition_mode}_{repeat_cycle_in_day:>04d}_{count_in_repeat_cycle:>04d}.bin'
+
+
+datasets:
+  cloud_mask:
+    name: cloud_mask
+    long_name: Cloud Classification
+    standard_name: cloud_classification
+    resolution: 2000
+    file_type: grib_fci_clm
+    parameter_number: 7
+    units: "1"
+    flag_values: [0, 1, 2, 3]
+    flag_meanings: ['clear sky over water','clear sky over land', 'cloudy', 'undefined' ]

satpy/etc/readers/seviri_l2_grib.yaml

@@ -14,7 +14,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Aerosol Properties over Sea product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:AES
   grib_seviri_aes:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'AESGRIBProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGAESE-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'
@@ -24,7 +24,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Cloud Mask product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:CLM
   grib_seviri_clm:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'CLMEncProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGCLMK-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'
@@ -34,7 +34,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Cloud Top Height product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:CTH
   grib_seviri_cth:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'CTHEncProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGCLTH-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'
@@ -44,7 +44,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Clear-Sky Reflectance Map product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:CRM
   grib_seviri_crm:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'CRMEncProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGCRMN-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'
@@ -54,7 +54,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Active Fire Monitoring product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:FIR
   grib_seviri_fir:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'FIREncProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGFIRG-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'
@@ -64,7 +64,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Multi-Sensor Precipitation Estimate product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:MPE-GRIB
   grib_seviri_mpe:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'MPEGRIBProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGMPEG-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'
@@ -74,7 +74,7 @@ file_types:
   # EUMETSAT MSG SEVIRI L2 Optimal Cloud Analysis product
   # https://navigator.eumetsat.int/product/EO:EUM:DAT:MSG:OCA
   grib_seviri_oca:
-    file_reader: !!python/name:satpy.readers.seviri_l2_grib.SeviriL2GribFileHandler
+    file_reader: !!python/name:satpy.readers.eum_l2_grib.EUML2GribFileHandler
     file_patterns:
       - 'OCAEncProd_{start_time:%Y%m%d%H%M%S}Z_00_{server:8s}_{spacecraft:5s}_{scan_mode:3s}_{sub_sat:5s}'
       - '{spacecraft:4s}-SEVI-MSGOCAE-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-{product_creation_time:%Y%m%d%H%M%S}-{ord_num:7s}'

satpy/readers/seviri_l2_grib.py → satpy/readers/eum_l2_grib.py

@@ -13,7 +13,7 @@
 # You should have received a copy of the GNU General Public License
 # along with satpy.  If not, see <http://www.gnu.org/licenses/>.
 
-"""Reader for the SEVIRI L2 products in GRIB2 format.
+"""Reader for both SEVIRI and FCI L2 products in GRIB2 format.
 
 References:
     FM 92 GRIB Edition 2
@@ -23,37 +23,49 @@
 """
 
 import logging
-from datetime import timedelta
 
 import dask.array as da
 import numpy as np
 import xarray as xr
 
 from satpy.readers._geos_area import get_area_definition, get_geos_area_naming
 from satpy.readers.eum_base import get_service_mode
+from satpy.readers.fci_base import calculate_area_extent as fci_calculate_area_extent
 from satpy.readers.file_handlers import BaseFileHandler
-from satpy.readers.seviri_base import PLATFORM_DICT, REPEAT_CYCLE_DURATION, calculate_area_extent
+from satpy.readers.seviri_base import PLATFORM_DICT, REPEAT_CYCLE_DURATION, REPEAT_CYCLE_DURATION_RSS
+from satpy.readers.seviri_base import calculate_area_extent as seviri_calculate_area_extent
 from satpy.utils import get_legacy_chunk_size
 
+CHUNK_SIZE = get_legacy_chunk_size()
+
 try:
     import eccodes as ec
 except ImportError:
     raise ImportError(
-        "Missing eccodes-python and/or eccodes C-library installation. Use conda to install eccodes")
+                "Missing eccodes-python and/or eccodes C-library installation. Use conda to install eccodes")
 
-CHUNK_SIZE = get_legacy_chunk_size()
 logger = logging.getLogger(__name__)
 
 
-class SeviriL2GribFileHandler(BaseFileHandler):
-    """Reader class for SEVIRI L2 products in GRIB format."""
+class EUML2GribFileHandler(BaseFileHandler):
+    """Reader class for EUM L2 products in GRIB format."""
+
+    calculate_area_extent = None
 
     def __init__(self, filename, filename_info, filetype_info):
         """Read the global attributes and prepare for dataset reading."""
         super().__init__(filename, filename_info, filetype_info)
         # Turn on support for multiple fields in single GRIB messages (required for SEVIRI L2 files)
         ec.codes_grib_multi_support_on()
 
+        if "seviri" in self.filetype_info["file_type"]:
+            self.sensor = "seviri"
+            self.PLATFORM_NAME = PLATFORM_DICT[self.filename_info["spacecraft"]]
+        elif "fci" in self.filetype_info["file_type"]:
+            self.sensor = "fci"
+            self.PLATFORM_NAME = f"MTG-i{self.filename_info['spacecraft_id']}"
+        pass
+
     @property
     def start_time(self):
         """Return the sensing start time."""
@@ -62,14 +74,21 @@ def start_time(self):
     @property
     def end_time(self):
         """Return the sensing end time."""
-        return self.start_time + timedelta(minutes=REPEAT_CYCLE_DURATION)
+        delta = REPEAT_CYCLE_DURATION_RSS if self._ssp_lon == 9.5 else REPEAT_CYCLE_DURATION
+        return self.start_time + delta
 
     def get_area_def(self, dataset_id):
         """Return the area definition for a dataset."""
+        # Compute the dictionary with the area extension
+
         self._area_dict["column_step"] = dataset_id["resolution"]
         self._area_dict["line_step"] = dataset_id["resolution"]
 
-        area_extent = calculate_area_extent(self._area_dict)
+        if self.sensor == "seviri":
+            area_extent = seviri_calculate_area_extent(self._area_dict)
+
+        elif self.sensor == "fci":
+            area_extent = fci_calculate_area_extent(self._area_dict)
 
         # Call the get_area_definition function to obtain the area
         area_def = get_area_definition(self._pdict, area_extent)
@@ -173,19 +192,20 @@ def _get_proj_area(self, gid):
         """
         # Get name of area definition
         area_naming_input_dict = {"platform_name": "msg",
-                                  "instrument_name": "seviri",
+                                  "instrument_name": self.sensor,
                                   "resolution": self._res,
                                   }
 
         area_naming = get_geos_area_naming({**area_naming_input_dict,
-                                            **get_service_mode("seviri", self._ssp_lon)})
+                                            **get_service_mode(self.sensor, self._ssp_lon)})
 
         # Read all projection and area parameters from the message
         earth_major_axis_in_meters = self._get_from_msg(gid, "earthMajorAxis") * 1000.0  # [m]
         earth_minor_axis_in_meters = self._get_from_msg(gid, "earthMinorAxis") * 1000.0  # [m]
 
-        earth_major_axis_in_meters = self._scale_earth_axis(earth_major_axis_in_meters)
-        earth_minor_axis_in_meters = self._scale_earth_axis(earth_minor_axis_in_meters)
+        if self.sensor == "seviri":
+            earth_major_axis_in_meters = self._scale_earth_axis(earth_major_axis_in_meters)
+            earth_minor_axis_in_meters = self._scale_earth_axis(earth_minor_axis_in_meters)
 
         nr_in_radius_of_earth = self._get_from_msg(gid, "NrInRadiusOfEarth")
         xp_in_grid_lengths = self._get_from_msg(gid, "XpInGridLengths")
@@ -204,25 +224,31 @@ def _get_proj_area(self, gid):
             "p_id": "",
         }
 
-        # Compute the dictionary with the area extension
-        area_dict = {
-            "center_point": xp_in_grid_lengths,
-            "north": self._nrows,
-            "east": 1,
-            "west": self._ncols,
-            "south": 1,
-        }
+        if self.sensor == "seviri":
+            # Compute the dictionary with the area extension
+            area_dict = {
+                "center_point": xp_in_grid_lengths,
+                "north": self._nrows,
+                "east": 1,
+                "west": self._ncols,
+                "south": 1,
+            }
+
+        elif self.sensor == "fci":
+            area_dict = {
+                "nlines": self._ncols,
+                "ncols": self._nrows,
+            }
 
         return pdict, area_dict
 
     @staticmethod
     def _scale_earth_axis(data):
         """Scale Earth axis data to make sure the value matched the expected unit [m].
 
-        The earthMinorAxis value stored in the aerosol over sea product is scaled incorrectly by a factor of 1e8. This
-        method provides a flexible temporarily workaraound by making sure that all earth axis values are scaled such
-        that they are on the order of millions of meters as expected by the reader. As soon as the scaling issue has
-        been resolved by EUMETSAT this workaround can be removed.
+        The earthMinorAxis value stored in the MPEF aerosol over sea product prior to December 12, 2022 has the wrong
+        unit and this method provides a flexible work-around by making sure that all earth axis values are scaled such
+        that they are on the order of millions of meters as expected by the reader.
 
         """
         scale_factor = 10 ** np.ceil(np.log10(1e6/data))
@@ -258,9 +284,12 @@ def _get_attributes(self):
 
         attributes = {
             "orbital_parameters": orbital_parameters,
-            "sensor": "seviri",
-            "platform_name": PLATFORM_DICT[self.filename_info["spacecraft"]]
+            "sensor": self.sensor
         }
+
+
+        attributes["platform_name"] = self.PLATFORM_NAME
+
         return attributes
 
     @staticmethod

satpy/readers/fci_base.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2017-2018 Satpy developers
+#
+# 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/>.
+"""Common functionality for FCI data readers."""
+from __future__ import annotations
+
+
+def calculate_area_extent(area_dict):
+        """Calculate the area extent seen by MTG FCI instrument.
+
+        Since the center of the FCI grids is located at the interface between the pixels, there are equally many
+        pixels (e.g. 5568/2 = 2784 for 2km grid) in each direction from the center points. Hence, the area extent
+        can be easily computed by simply adding and subtracting half the width and height from teh centre point (=0).
+
+        Args:
+            area_dict: A dictionary containing the required parameters
+                ncols: number of pixels in east-west direction
+                nlines: number of pixels in south-north direction
+                column_step: Pixel resulution in meters in east-west direction
+                line_step: Pixel resulution in meters in south-north direction
+        Returns:
+            tuple: An area extent for the scene defined by the lower left and
+                   upper right corners
+
+        """
+        ncols = area_dict["ncols"]
+        nlines = area_dict["nlines"]
+        column_step = area_dict["column_step"]
+        line_step = area_dict["line_step"]
+
+        ll_c = (0 - ncols / 2.) * column_step
+        ll_l = (0 + nlines / 2.) * line_step
+        ur_c = (0 + ncols / 2.) * column_step
+        ur_l = (0 - nlines / 2.) * line_step
+
+        return (ll_c, ll_l, ur_c, ur_l)

satpy/readers/seviri_base.py

@@ -212,6 +212,8 @@
 
 REPEAT_CYCLE_DURATION = 15
 
+REPEAT_CYCLE_DURATION_RSS = 5
+
 C1 = 1.19104273e-5
 C2 = 1.43877523