# Copyright (c) 2019-2023 Satpy developers
#
# 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/>.
"""Reader for the SEVIRI L2 products in GRIB2 format.
References:
FM 92 GRIB Edition 2
https://www.wmo.int/pages/prog/www/WMOCodes/Guides/GRIB/GRIB2_062006.pdf
EUMETSAT Product Navigator
https://navigator.eumetsat.int/
"""
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.file_handlers import BaseFileHandler
from satpy.readers.seviri_base import PLATFORM_DICT, REPEAT_CYCLE_DURATION, calculate_area_extent
from satpy.utils import 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")
CHUNK_SIZE = get_legacy_chunk_size()
logger = logging.getLogger(__name__)
[docs]
class SeviriL2GribFileHandler(BaseFileHandler):
"""Reader class for SEVIRI L2 products in GRIB format."""
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()
@property
def start_time(self):
"""Return the sensing start time."""
return self.filename_info["start_time"]
@property
def end_time(self):
"""Return the sensing end time."""
return self.start_time + timedelta(minutes=REPEAT_CYCLE_DURATION)
[docs]
def get_area_def(self, dataset_id):
"""Return the area definition for a dataset."""
self._area_dict["column_step"] = dataset_id["resolution"]
self._area_dict["line_step"] = dataset_id["resolution"]
area_extent = 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)
return area_def
[docs]
def get_dataset(self, dataset_id, dataset_info):
"""Get dataset using the parameter_number key in dataset_info.
In a previous version of the reader, the attributes (nrows, ncols, ssp_lon) and projection information
(pdict and area_dict) were computed while initializing the file handler. Also the code would break out from
the While-loop below as soon as the correct parameter_number was found. This has now been revised becasue the
reader would sometimes give corrupt information about the number of messages in the file and the dataset
dimensions within a given message if the file was only partly read (not looping over all messages) in an earlier
instance.
"""
logger.debug("Reading in file to get dataset with parameter number %d.",
dataset_info["parameter_number"])
xarr = None
message_found = False
with open(self.filename, "rb") as fh:
# Iterate over all messages and fetch data when the correct parameter number is found
while True:
gid = ec.codes_grib_new_from_file(fh)
if gid is None:
if not message_found:
# Could not obtain a valid message ID from the grib file
logger.warning("Could not find parameter_number %d in GRIB file, no valid Dataset created",
dataset_info["parameter_number"])
break
# Check if the parameter number in the GRIB message corresponds to the required key
parameter_number = self._get_from_msg(gid, "parameterNumber")
if parameter_number == dataset_info["parameter_number"]:
self._res = dataset_id["resolution"]
self._read_attributes(gid)
# Read the missing value
missing_value = self._get_from_msg(gid, "missingValue")
# Retrieve values and metadata from the GRIB message, masking the values equal to missing_value
xarr = self._get_xarray_from_msg(gid)
xarr.data = da.where(xarr.data == missing_value, np.nan, xarr.data)
ec.codes_release(gid)
# Combine all metadata into the dataset attributes and break out of the loop
xarr.attrs.update(dataset_info)
xarr.attrs.update(self._get_attributes())
message_found = True
else:
# The parameter number is not the correct one, release gid and skip to next message
ec.codes_release(gid)
return xarr
[docs]
def _read_attributes(self, gid):
"""Read the parameter attributes from the message and create the projection and area dictionaries."""
# Read SSP and date/time
self._ssp_lon = self._get_from_msg(gid, "longitudeOfSubSatellitePointInDegrees")
# Read number of points on the x and y axes
self._nrows = self._get_from_msg(gid, "Ny")
self._ncols = self._get_from_msg(gid, "Nx")
# Creates the projection and area dictionaries
self._pdict, self._area_dict = self._get_proj_area(gid)
[docs]
def _get_proj_area(self, gid):
"""Compute the dictionary with the projection and area definition from a GRIB message.
Args:
gid: The ID of the GRIB message.
Returns:
tuple: A tuple of two dictionaries for the projection and the area definition.
pdict:
a: Earth major axis [m]
b: Earth minor axis [m]
h: Height over surface [m]
ssp_lon: longitude of subsatellite point [deg]
nlines: number of lines
ncols: number of columns
a_name: name of the area
a_desc: description of the area
p_id: id of the projection
area_dict:
center_point: coordinate of the center point
north: coodinate of the north limit
east: coodinate of the east limit
west: coodinate of the west limit
south: coodinate of the south limit
"""
# Get name of area definition
area_naming_input_dict = {"platform_name": "msg",
"instrument_name": "seviri",
"resolution": self._res,
}
area_naming = get_geos_area_naming({**area_naming_input_dict,
**get_service_mode("seviri", 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)
nr_in_radius_of_earth = self._get_from_msg(gid, "NrInRadiusOfEarth")
xp_in_grid_lengths = self._get_from_msg(gid, "XpInGridLengths")
h_in_meters = earth_major_axis_in_meters * (nr_in_radius_of_earth - 1.0) # [m]
# Create the dictionary with the projection data
pdict = {
"a": earth_major_axis_in_meters,
"b": earth_minor_axis_in_meters,
"h": h_in_meters,
"ssp_lon": self._ssp_lon,
"nlines": self._ncols,
"ncols": self._nrows,
"a_name": area_naming["area_id"],
"a_desc": area_naming["description"],
"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,
}
return pdict, area_dict
[docs]
@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.
"""
scale_factor = 10 ** np.ceil(np.log10(1e6/data))
return data * scale_factor
[docs]
def _get_xarray_from_msg(self, gid):
"""Read the values from the GRIB message and return a DataArray object.
Args:
gid: The ID of the GRIB message.
Returns:
DataArray: The array containing the retrieved values.
"""
# Data from GRIB message are read into an Xarray...
xarr = xr.DataArray(da.from_array(ec.codes_get_values(
gid).reshape(self._nrows, self._ncols), CHUNK_SIZE), dims=("y", "x"))
return xarr
[docs]
def _get_attributes(self):
"""Create a dictionary of attributes to be added to the dataset.
Returns:
dict: A dictionary of parameter attributes.
ssp_lon: longitude of subsatellite point
sensor: name of sensor
platform_name: name of the platform
"""
orbital_parameters = {
"projection_longitude": self._ssp_lon
}
attributes = {
"orbital_parameters": orbital_parameters,
"sensor": "seviri",
"platform_name": PLATFORM_DICT[self.filename_info["spacecraft"]]
}
return attributes
[docs]
@staticmethod
def _get_from_msg(gid, key):
"""Get a value from the GRIB message based on the key, return None if missing.
Args:
gid: The ID of the GRIB message.
key: The key of the required attribute.
Returns:
The retrieved attribute or None if the key is missing.
"""
try:
attr = ec.codes_get(gid, key)
except ec.KeyValueNotFoundError:
logger.warning("Key %s not found in GRIB message", key)
attr = None
return attr