Harmonise AreaDefinition namings in EUM geos readers and sort geos areas in areas.yaml

satpy/readers/_geos_area.py

@@ -127,6 +127,10 @@ def get_area_definition(pdict, a_ext):
     Returns:
         a_def: An area definition for the scene
 
+    .. note::
+
+        The AreaDefinition `proj_id` attribute is being deprecated.
+
     """
     proj_dict = {'a': float(pdict['a']),
                  'b': float(pdict['b']),
@@ -154,9 +158,8 @@ def sampling_to_lfac_cfac(sampling):
     Appendix E.2.
 
 
-    .. _MSG Ground Segment LRIT HRIT Mission Specific Implementation: http://www.eumetsat.int/\
-website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05057_SPE_MSG_LRIT_HRI\
-&RevisionSelectionMethod=LatestReleased&Rendition=Web
+    .. _MSG Ground Segment LRIT HRIT Mission Specific Implementation:
+    https://www-cdn.eumetsat.int/files/2020-04/pdf_ten_05057_spe_msg_lrit_hri.pdf
 
     Args:
         sampling: float
@@ -166,3 +169,61 @@ def sampling_to_lfac_cfac(sampling):
         Line/column scaling factor (deg-1)
     """
     return 2.0 ** 16 / np.rad2deg(sampling)
+
+
+def get_geos_area_naming(input_dict):
+    """Get a dictionary containing formatted AreaDefinition naming.
+
+    Args:
+        input_dict: dict
+            Dictionary with keys `platform_name`, `instrument_name`, `service_name`, `service_desc`, `resolution` .
+            The resolution is expected in meters.
+    Returns:
+        area_naming_dict with `area_id`, `description`  keys, values are strings.
+
+    .. note::
+
+        The AreaDefinition `proj_id` attribute is being deprecated and is therefore not formatted here.
+        An empty string is to be used until the attribute is fully removed.
+
+    """
+    area_naming_dict = {}
+
+    resolution_strings = get_resolution_and_unit_strings(input_dict['resolution'])
+
+    area_naming_dict['area_id'] = '{}_{}_{}_{}{}'.format(input_dict['platform_name'].lower(),
+                                                         input_dict['instrument_name'].lower(),
+                                                         input_dict['service_name'].lower(),
+                                                         resolution_strings['value'],
+                                                         resolution_strings['unit']
+                                                         )
+
+    area_naming_dict['description'] = '{} {} {} area definition ' \
+                                      'with {} {} resolution'.format(input_dict['platform_name'].upper(),
+                                                                     input_dict['instrument_name'].upper(),
+                                                                     input_dict['service_desc'],
+                                                                     resolution_strings['value'],
+                                                                     resolution_strings['unit']
+                                                                     )
+
+    return area_naming_dict
+
+
+def get_resolution_and_unit_strings(resolution):
+    """Get the resolution value and unit as strings.
+
+    If the resolution is larger than 1000 m, use kilometer as unit. If lower, use meter.
+
+    Args:
+        resolution: scalar
+            Resolution in meters.
+
+    Returns:
+        Dictionary with `value` and `unit` keys, values are strings.
+    """
+    if resolution >= 1000:
+        return {'value': '{:.0f}'.format(resolution*1e-3),
+                'unit': 'km'}
+
+    return {'value': '{:.0f}'.format(resolution),
+            'unit': 'm'}

satpy/readers/eum_base.py

@@ -82,3 +82,18 @@ def recarray2dict(arr):
                 res[key] = data.squeeze()
 
     return res
+
+
+def get_service_mode(instrument_name, ssp_lon):
+    """Get information about service mode for a given instrument and subsatellite longitude."""
+    service_modes = {'seviri': {'0.0':  {'service_name': 'fes', 'service_desc': 'Full Earth Scanning service'},
+                                '9.5':  {'service_name': 'rss', 'service_desc': 'Rapid Scanning Service'},
+                                '41.5': {'service_name': 'iodc', 'service_desc': 'Indian Ocean Data Coverage service'}
+                                },
+                     'fci':    {'0.0':  {'service_name': 'fdss', 'service_desc': 'Full Disk Scanning Service'},
+                                '9.5':  {'service_name': 'rss', 'service_desc': 'Rapid Scanning Service'},
+                                },
+                     }
+    unknown_modes = {'service_name': 'unknown', 'service_desc': 'unknown'}
+
+    return service_modes.get(instrument_name, unknown_modes).get('{:.1f}'.format(ssp_lon), unknown_modes)

satpy/readers/fci_l1c_fdhsi.py

@@ -84,6 +84,8 @@
 
 from pyresample import geometry
 from netCDF4 import default_fillvals
+from satpy.readers._geos_area import get_geos_area_naming
+from satpy.readers.eum_base import get_service_mode
 
 from .netcdf_utils import NetCDF4FileHandler
 
@@ -337,15 +339,17 @@ def get_area_def(self, key, info=None):
                      'units': 'm',
                      "sweep": sweep}
 
-        # follow the format <platform>_<instrument>_<service>_<resolution>
-        area_id_name = 'mtg_fci_fdss_{:.0f}km'.format(key['resolution']/1000)
-        area_id_description = 'FCI Full-Disk Scanning Service area definition with {:.0f} km resolution' \
-                              ''.format(key['resolution']/1000)
+        area_naming_input_dict = {'platform_name': 'mtg',
+                                  'instrument_name': 'fci',
+                                  'resolution': int(key['resolution'])
+                                  }
+        area_naming = get_geos_area_naming({**area_naming_input_dict,
+                                            **get_service_mode('fci', lon_0)})
 
         area = geometry.AreaDefinition(
-            area_id_name,
-            area_id_description,
-            "FCI Full-Disk Scanning Service geostationary projection",
+            area_naming['area_id'],
+            area_naming['description'],
+            "",
             proj_dict,
             ncols,
             nlines,

satpy/readers/seviri_base.py

@@ -635,15 +635,6 @@ def create_coef_dict(coefs_nominal, coefs_gsics, radiance_type, ext_coefs):
     }
 
 
-def get_service_mode(ssp_lon):
-    """Get information about whether we're dealing with data from the FES, RSS or IODC service mode."""
-    seviri_service_modes = {'0.0': {'name': 'FES', 'desc': 'Full Earth scanning service'},
-                            '9.5': {'name': 'RSS', 'desc': 'Rapid scanning service'},
-                            '41.5': {'name': 'IODC', 'desc': 'Indian Ocean data coverage service'}}
-
-    return seviri_service_modes['%.1f' % ssp_lon]
-
-
 def get_padding_area(shape, dtype):
     """Create a padding area filled with no data."""
     if np.issubdtype(dtype, np.floating):

satpy/readers/seviri_l1b_hrit.py

@@ -135,7 +135,7 @@
 from pyresample import geometry
 from satpy import CHUNK_SIZE
 import satpy.readers.utils as utils
-from satpy.readers.eum_base import recarray2dict, time_cds_short
+from satpy.readers.eum_base import recarray2dict, time_cds_short, get_service_mode
 from satpy.readers.hrit_base import (HRITFileHandler, ancillary_text,
                                      annotation_header, base_hdr_map,
                                      image_data_function)
@@ -146,8 +146,7 @@
                                        pad_data_horizontally, create_coef_dict)
 from satpy.readers.seviri_l1b_native_hdr import (hrit_epilogue, hrit_prologue,
                                                  impf_configuration)
-from satpy.readers._geos_area import get_area_extent, get_area_definition
-
+from satpy.readers._geos_area import get_area_extent, get_area_definition, get_geos_area_naming
 
 logger = logging.getLogger('hrit_msg')
 
@@ -556,13 +555,20 @@ def get_area_def(self, dsid):
         else:
             pdict['scandir'] = 'S2N'
 
+        area_naming_input_dict = {'platform_name': 'msg',
+                                  'instrument_name': 'seviri',
+                                  'resolution': int(dsid['resolution'])
+                                  }
+        area_naming = get_geos_area_naming({**area_naming_input_dict,
+                                            **get_service_mode('seviri', pdict['ssp_lon'])})
+
         # Compute area definition for non-HRV channels:
         if dsid['name'] != 'HRV':
             pdict['loff'] = loff - nlines
             aex = self._get_area_extent(pdict)
-            pdict['a_name'] = 'geosmsg'
-            pdict['a_desc'] = 'MSG/SEVIRI low resolution channel area'
-            pdict['p_id'] = 'msg_lowres'
+            pdict['a_name'] = area_naming['area_id']
+            pdict['a_desc'] = area_naming['description']
+            pdict['p_id'] = ""
             area = get_area_definition(pdict, aex)
             self.area = area
             return self.area
@@ -574,8 +580,8 @@ def get_area_def(self, dsid):
                               * pdict['nlines'])
 
         # Or, if we are processing HRV:
-        pdict['a_name'] = 'geosmsg_hrv'
-        pdict['p_id'] = 'msg_hires'
+        pdict['a_name'] = area_naming['area_id']
+        pdict['p_id'] = ""
         bounds = self.epilogue['ImageProductionStats']['ActualL15CoverageHRV'].copy()
         if self.fill_hrv:
             bounds['UpperEastColumnActual'] = 1
@@ -594,15 +600,15 @@ def get_area_def(self, dsid):
         pdict['nlines'] = upper_south_line
         pdict['loff'] = loff - upper_south_line
         pdict['coff'] = lower_coff
-        pdict['a_desc'] = 'MSG/SEVIRI high resolution channel, lower window'
+        pdict['a_desc'] = area_naming['description']
         lower_area_extent = self._get_area_extent(pdict)
         lower_area = get_area_definition(pdict, lower_area_extent)
 
         # Now the upper window
         pdict['nlines'] = nlines - upper_south_line
         pdict['loff'] = loff - pdict['nlines'] - upper_south_line
         pdict['coff'] = upper_coff
-        pdict['a_desc'] = 'MSG/SEVIRI high resolution channel, upper window'
+        pdict['a_desc'] = area_naming['description']
         upper_area_extent = self._get_area_extent(pdict)
         upper_area = get_area_definition(pdict, upper_area_extent)
 

satpy/readers/seviri_l1b_native.py

@@ -37,16 +37,16 @@
 from pyresample import geometry
 
 from satpy.readers.file_handlers import BaseFileHandler
-from satpy.readers.eum_base import recarray2dict
+from satpy.readers.eum_base import recarray2dict, get_service_mode
 from satpy.readers.seviri_base import (
     SEVIRICalibrationHandler, CHANNEL_NAMES, SATNUM, dec10216,
     VISIR_NUM_COLUMNS, VISIR_NUM_LINES, HRV_NUM_COLUMNS, HRV_NUM_LINES,
-    create_coef_dict, get_service_mode, pad_data_horizontally,
+    create_coef_dict, pad_data_horizontally,
     pad_data_vertically
 )
 from satpy.readers.seviri_l1b_native_hdr import (GSDTRecords, native_header,
                                                  native_trailer)
-from satpy.readers._geos_area import get_area_definition
+from satpy.readers._geos_area import get_area_definition, get_geos_area_naming
 
 logger = logging.getLogger('native_msg')
 
@@ -267,23 +267,38 @@ def get_area_def(self, dataset_id):
         and corresponding number of image lines/columns. In case of FES HRV data, two area definitions are
         computed, stacked and squeezed. For other cases, the lists will only have one entry each, from which
         a single area definition is computed.
+
+        Note that the AreaDefinition area extents returned by this function for Native data will be slightly
+        different compared to the area extents returned by the SEVIRI HRIT reader.
+        This is due to slightly different pixel size values when calculated using the data available in the files. E.g.
+        for the 3 km grid:
+
+        ``Native: data15hd['ImageDescription']['ReferenceGridVIS_IR']['ColumnDirGridStep'] == 3000.4031658172607``
+        ``HRIT:                            np.deg2rad(2.**16 / pdict['lfac']) * pdict['h'] == 3000.4032785810186``
+
+        This results in the Native 3 km full-disk area extents being approx. 20 cm shorter in each direction.
+
+        The method for calculating the area extents used by the HRIT reader (CFAC/LFAC mechanism) keeps the
+        highest level of numeric precision and is used as reference by EUM. For this reason, the standard area
+        definitions defined in the `areas.yaml` file correspond to the HRIT ones.
+
         """
         pdict = {}
         pdict['a'] = self.mda['projection_parameters']['a']
         pdict['b'] = self.mda['projection_parameters']['b']
         pdict['h'] = self.mda['projection_parameters']['h']
         pdict['ssp_lon'] = self.mda['projection_parameters']['ssp_longitude']
 
-        if dataset_id['name'] == 'HRV':
-            res = 1.0
-            pdict['p_id'] = 'seviri_hrv'
-        else:
-            res = 3.0
-            pdict['p_id'] = 'seviri_visir'
+        area_naming_input_dict = {'platform_name': 'msg',
+                                  'instrument_name': 'seviri',
+                                  'resolution': int(dataset_id['resolution'])
+                                  }
+        area_naming = get_geos_area_naming({**area_naming_input_dict,
+                                            **get_service_mode('seviri', pdict['ssp_lon'])})
 
-        service_mode = get_service_mode(pdict['ssp_lon'])
-        pdict['a_name'] = 'msg_seviri_%s_%.0fkm' % (service_mode['name'], res)
-        pdict['a_desc'] = 'SEVIRI %s area definition with %.0f km resolution' % (service_mode['desc'], res)
+        pdict['a_name'] = area_naming['area_id']
+        pdict['a_desc'] = area_naming['description']
+        pdict['p_id'] = ""
 
         area_extent = self.get_area_extent(dataset_id)
         areas = list()

satpy/readers/seviri_l1b_nc.py

@@ -20,9 +20,10 @@
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers.seviri_base import (SEVIRICalibrationHandler,
                                        CHANNEL_NAMES, SATNUM)
+from satpy.readers.eum_base import get_service_mode
 import xarray as xr
 
-from satpy.readers._geos_area import get_area_definition
+from satpy.readers._geos_area import get_area_definition, get_geos_area_naming
 from satpy import CHUNK_SIZE
 
 import datetime
@@ -174,25 +175,48 @@ def _get_calib_coefs(self, dataset, channel):
         }
 
     def get_area_def(self, dataset_id):
-        """Get the area def."""
+        """Get the area def.
+
+        Note that the AreaDefinition area extents returned by this function for NetCDF data will be slightly
+        different compared to the area extents returned by the SEVIRI HRIT reader.
+        This is due to slightly different pixel size values when calculated using the data available in the files. E.g.
+        for the 3 km grid:
+
+        ``NetCDF:  self.nc.attrs['vis_ir_column_dir_grid_step'] == 3000.4031658172607``
+        ``HRIT: np.deg2rad(2.**16 / pdict['lfac']) * pdict['h'] == 3000.4032785810186``
+
+        This results in the Native 3 km full-disk area extents being approx. 20 cm shorter in each direction.
+
+        The method for calculating the area extents used by the HRIT reader (CFAC/LFAC mechanism) keeps the
+        highest level of numeric precision and is used as reference by EUM. For this reason, the standard area
+        definitions defined in the `areas.yaml` file correspond to the HRIT ones.
+
+        """
         pdict = {}
         pdict['a'] = self.mda['projection_parameters']['a']
         pdict['b'] = self.mda['projection_parameters']['b']
         pdict['h'] = self.mda['projection_parameters']['h']
         pdict['ssp_lon'] = self.mda['projection_parameters']['ssp_longitude']
 
+        area_naming_input_dict = {'platform_name': 'msg',
+                                  'instrument_name': 'seviri',
+                                  'resolution': int(dataset_id['resolution'])
+                                  }
+        area_naming = get_geos_area_naming({**area_naming_input_dict,
+                                            **get_service_mode('seviri', pdict['ssp_lon'])})
+
         if dataset_id['name'] == 'HRV':
             pdict['nlines'] = self.mda['hrv_number_of_lines']
             pdict['ncols'] = self.mda['hrv_number_of_columns']
-            pdict['a_name'] = 'geosmsg_hrv'
-            pdict['a_desc'] = 'MSG/SEVIRI high resolution channel area'
-            pdict['p_id'] = 'msg_hires'
+            pdict['a_name'] = area_naming['area_id']
+            pdict['a_desc'] = area_naming['description']
+            pdict['p_id'] = ""
         else:
             pdict['nlines'] = self.mda['number_of_lines']
             pdict['ncols'] = self.mda['number_of_columns']
-            pdict['a_name'] = 'geosmsg'
-            pdict['a_desc'] = 'MSG/SEVIRI low resolution channel area'
-            pdict['p_id'] = 'msg_lowres'
+            pdict['a_name'] = area_naming['area_id']
+            pdict['a_desc'] = area_naming['description']
+            pdict['p_id'] = ""
 
         area = get_area_definition(pdict, self.get_area_extent(dataset_id))