diff --git a/AUTHORS.md b/AUTHORS.md index 9078e441b4..1d98572541 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -87,3 +87,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) diff --git a/satpy/etc/readers/eum_l2_grib.yaml b/satpy/etc/readers/eum_l2_grib.yaml new file mode 100644 index 0000000000..80edd3b2e5 --- /dev/null +++ b/satpy/etc/readers/eum_l2_grib.yaml @@ -0,0 +1,387 @@ +reader: + name: eum_l2_grib + short_name: EUM L2 GRIB + long_name: MSG (Meteosat 8 to 11) SEVIRI Level products 2 and FCI L2 products in GRIB2 format + description: Reader for EUMETSAT MSG SEVIRI L2 files and FCI L2 files in GRIB format. + status: Alpha + supports_fsspec: false + sensors: [seviri,fci] + reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader + + +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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGAESE-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + # 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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGCLMK-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + # 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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGCLTH-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + # 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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGCRMN-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + # 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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGFIRG-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + # 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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGMPEG-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + # 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.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}' + - '{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}.grb' + - '{spacecraft:4s}-SEVI-MSGOCAE-{id1:4s}-{id2:4s}-{start_time:%Y%m%d%H%M%S}.000000000Z-NA.grb' + + 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: + + # EUMETSAT MSG SEVIRI L2 Aerosol Properties over Sea product + aerosol_optical_thickness_vis06: + name: aerosol_optical_thickness_vis06 + long_name: Aerosol optical Thickness at 0.6um + standard_name: atmosphere_absorption_optical_thickness_due_to_ambient_aerosol + resolution: 9001.209497451 + file_type: grib_seviri_aes + parameter_number: 20 + units: "1" + + aerosol_optical_thickness_vis08: + name: aerosol_optical_thickness_vis08 + long_name: Aerosol optical Thickness at 0.8um + standard_name: atmosphere_absorption_optical_thickness_due_to_ambient_aerosol + resolution: 9001.209497451 + file_type: grib_seviri_aes + parameter_number: 21 + units: "1" + + aerosol_optical_thickness_vis16: + name: aerosol_optical_thickness_vis16 + long_name: Aerosol optical Thickness at 1.6um + standard_name: atmosphere_absorption_optical_thickness_due_to_ambient_aerosol + resolution: 9001.209497451 + file_type: grib_seviri_aes + parameter_number: 22 + units: "1" + + angstroem_coefficient: + name: angstroem_coefficient + long_name: Angstroem Coefficient + standard_name: aerosol_angstrom_exponent + resolution: 9001.209497451 + file_type: grib_seviri_aes + parameter_number: 23 + units: "1" + + aes_quality: + name: aes_quality + long_name: AES Product Quality Flag + standard_name: quality_flag + resolution: 9001.209497451 + file_type: grib_seviri_aes + parameter_number: 192 + units: "1" + flag_values: [0, 1, 2, 3] + flag_meanings: ['clear sky over water','clear sky over land', 'cloudy', 'no data' ] + + + # EUMETSAT MSG SEVIRI L2 Cloud Mask product + cloud_mask: + name: cloud_mask + long_name: Cloud Classification + standard_name: cloud_classification + resolution: + 3000.403165817: {file_type: grib_seviri_clm} + 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', 'no data' ] + + + # EUMETSAT MSG SEVIRI L2 Cloud Top Height product + cloud_top_height: + name: cloud_top_height + long_name: Cloud Top Height + standard_name: height_at_cloud_top + resolution: 9001.209497451 + file_type: grib_seviri_cth + parameter_number: 2 + units: m + + cloud_top_quality: + name: cloud_top_quality + long_name: CTH Product Quality Flag + standard_name: height_at_cloud_top quality_flag + resolution: 9001.209497451 + file_type: grib_seviri_cth + parameter_number: 3 + units: "1" + flag_values: [0, 1] + flag_meanings: ['good quality retrieval','poor quality retrieval' ] + + + # EUMETSAT MSG SEVIRI L2 Clear-Sky Reflectance Map product + vis_refl_06: + name: vis_refl_06 + long_name: TOA Bidirectional Reflectance at 0.6um (7 days average) + standard_name: toa_bidirectional_reflectance + resolution: 3000.403165817 + wavelength: [0.56, 0.635, 0.71] + file_type: grib_seviri_crm + parameter_number: 9 + units: "%" + + vis_refl_08: + name: vis_refl_08 + long_name: TOA Bidirectional Reflectance at 0.8um (7 days average) + standard_name: toa_bidirectional_reflectance + resolution: 3000.403165817 + wavelength: [0.74, 0.81, 0.88] + file_type: grib_seviri_crm + parameter_number: 10 + units: "%" + + vis_refl_16: + name: vis_refl_16 + long_name: TOA Bidirectional Reflectance at 1.6um (7 days average) + standard_name: toa_bidirectional_reflectance + resolution: 3000.403165817 + wavelength: [1.5, 1.64, 1.78] + file_type: grib_seviri_crm + parameter_number: 11 + units: "%" + + nir_refl_39: + name: nir_refl_39 + long_name: TOA Bidirectional Reflectance at 3.9um (7 days average) + standard_name: toa_bidirectional_reflectance + resolution: 3000.403165817 + wavelength: [3.48, 3.92, 4.36] + file_type: grib_seviri_crm + parameter_number: 12 + units: "%" + + num_accumulations: + name: num_accumulations + long_name: CRM Product Number of Accumulations + standard_name: number_of_accumulations + resolution: 3000.403165817 + file_type: grib_seviri_crm + parameter_number: 6 + units: "1" + + solar_zenith_angle: + name: solar_zenith_angle + long_name: Solar Zenith Angle (7 days average) + standard_name: solar_zenith_angle + resolution: 3000.403165817 + file_type: grib_seviri_crm + parameter_number: 7 + units: degrees + + relative_azimuth_angle: + name: relative_azimuth_angle + long_name: Relative Azimuth Angle (7 days average) + standard_name: relative_sensor_azimuth_angle + resolution: 3000.403165817 + file_type: grib_seviri_crm + parameter_number: 8 + units: degrees + + + # EUMETSAT MSG SEVIRI L2 Active Fire Monitoring product + active_fires: + name: active_fires + long_name: Active Fire Classification + standard_name: active_fire_classification + resolution: 3000.403165817 + file_type: grib_seviri_fir + parameter_number: 9 + units: "1" + flag_values: [0, 1, 2, 3] + flag_meanings: ['no fire','possible fire', 'probable fire', 'missing' ] + + fire_probability: + name: fire_probability + long_name: Fire Probability + standard_name: fire_probability + resolution: 3000.403165817 + file_type: grib_seviri_fir + parameter_number: 192 + units: "%" + + + # EUMETSAT MSG SEVIRI L2 Multi-Sensor Precipitation Estimate product + instantaneous_rain_rate: + name: instantaneous_rain_rate + long_name: MPE Product Instantaneous Rain Rate + standard_name: rainfall_rate + resolution: 3000.403165817 + file_type: grib_seviri_mpe + parameter_number: 1 + units: "kg m-2 s-1" + + + # EUMETSAT MSG SEVIRI L2 Optimal Cloud Analysis product + pixel_scene_type: + name: pixel_scene_type + long_name: Cloud Type + standard_name: scene_classification + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 8 + units: "1" + flag_values: [24,111,112] + flag_meanings: ['multi-layered cloud','water cloud','ice cloud'] + + measurement_cost: + name: measurement_cost + long_name: OCA Cost Function - Measurement part + standard_name: cost_function + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 30 + units: "1" + + upper_layer_cloud_optical_depth: + name: upper_layer_cloud_optical_depth + long_name: Upper Cloud Layer Optical Depth + standard_name: atmosphere_optical_thickness_due_to_cloud + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 31 + units: "1" + + upper_layer_cloud_top_pressure: + name: upper_layer_cloud_top_pressure + long_name: Upper Cloud Top Pressure + standard_name: air_pressure_at_cloud_top + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 32 + units: Pa + + upper_layer_cloud_effective_radius: + name: upper_layer_cloud_effective_radius + long_name: Upper Cloud Particle Effective Radius + standard_name: effective_radius_of_cloud_condensed_water_particles_at_cloud_top + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 33 + units: m + + error_in_upper_layer_cloud_optical_depth: + name: error_in_upper_layer_cloud_optical_depth + long_name: Upper Cloud Optical Depth Error Estimate + standard_name: atmosphere_optical_thickness_due_to_cloud standard_error + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 34 + units: "1" + + error_in_upper_layer_cloud_top_pressure: + name: error_in_upper_layer_cloud_top_pressure + long_name: Upper Cloud Top Pressure Error Estimate + standard_name: air_pressure_at_cloud_top standard_error + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 35 + units: Pa + + error_in_upper_layer_cloud_effective_radius: + name: error_in_upper_layer_cloud_effective_radius + long_name: Upper Cloud Particle Effective Radius Error Estimate + standard_name: effective_radius_of_cloud_condensed_water_particles_at_cloud_top standard_error + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 36 + units: m + + lower_layer_cloud_optical_depth: + name: lower_layer_cloud_optical_depth + long_name: Lower Cloud Optical Depth + standard_name: atmosphere_optical_thickness_due_to_cloud_in_lower_atmosphere_layer + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 37 + units: "1" + + lower_layer_cloud_top_pressure: + name: lower_layer_cloud_top_pressure + long_name: Lower Cloud Top Pressure + standard_name: air_pressure_at_cloud_top_in_lower_atmosphere_layer + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 38 + units: Pa + + error_in_lower_layer_cloud_optical_depth: + name: error_in_lower_layer_cloud_optical_depth + long_name: Lower Cloud Optical Depth Error Estimate + standard_name: atmosphere_optical_thickness_due_to_cloud_in_lower_atmosphere_layer standard_error + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 39 + units: "1" + + error_in_lower_layer_cloud_top_pressure: + name: error_in_lower_layer_cloud_top_pressure + long_name: Lower Cloud Top Pressure Error Estimate + standard_name: air_pressure_at_cloud_top_in_lower_atmosphere_layer standard_error + resolution: 3000.403165817 + file_type: grib_seviri_oca + parameter_number: 40 + units: Pa diff --git a/satpy/readers/eum_l2_grib.py b/satpy/readers/eum_l2_grib.py new file mode 100644 index 0000000000..47cf9a0ba9 --- /dev/null +++ b/satpy/readers/eum_l2_grib.py @@ -0,0 +1,315 @@ +# 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 . + +"""Reader for both SEVIRI and FCI 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.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 +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") + +logger = logging.getLogger(__name__) + + +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.""" + return self.filename_info["start_time"] + + @property + def end_time(self): + """Return the sensing end time.""" + if self.sensor == "seviri": + return self.start_time + timedelta(minutes=REPEAT_CYCLE_DURATION) + elif self.sensor == "fci": + return self.filename_info["end_time"] + + 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"] + + 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) + + return area_def + + 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 + + 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) + + 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": self.sensor, + "resolution": self._res, + } + + area_naming = get_geos_area_naming({**area_naming_input_dict, + **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] + + 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") + 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": "", + } + + 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. + + """ + scale_factor = 10 ** np.ceil(np.log10(1e6/data)) + return data * scale_factor + + 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 + + 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": self.sensor + } + + + attributes["platform_name"] = self.PLATFORM_NAME + + return attributes + + @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 diff --git a/satpy/readers/fci_base.py b/satpy/readers/fci_base.py new file mode 100644 index 0000000000..c4a3714291 --- /dev/null +++ b/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 . +"""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 L2 grid 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) diff --git a/satpy/tests/reader_tests/test_eum_l2_grib.py b/satpy/tests/reader_tests/test_eum_l2_grib.py new file mode 100644 index 0000000000..3e4dee87a8 --- /dev/null +++ b/satpy/tests/reader_tests/test_eum_l2_grib.py @@ -0,0 +1,319 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (c) 2019 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 . + +"""EUM L2 GRIB-reader test package.""" + +import datetime +import sys +import unittest +from unittest import mock + +import numpy as np + +from satpy.tests.utils import make_dataid + +# Dictionary to be used as fake GRIB message +FAKE_SEVIRI_MESSAGE = { + "longitudeOfSubSatellitePointInDegrees": 9.5, + "dataDate": 20191020, + "dataTime": 1745, + "Nx": 1000, + "Ny": 1200, + "earthMajorAxis": 6400., + "earthMinorAxis": 6300., + "NrInRadiusOfEarth": 6., + "XpInGridLengths": 500, + "parameterNumber": 30, + "missingValue": 9999, +} + +FAKE_FCI_MESSAGE = { + "longitudeOfSubSatellitePointInDegrees": 0.0, + "dataDate": 20191020, + "dataTime": 1745, + "Nx": 5568, + "Ny": 5568, + "earthMajorAxis": 6378140., + "earthMinorAxis": 6356755., + "NrInRadiusOfEarth": 6.6107, + "XpInGridLengths": 2784.0, + "parameterNumber": 30, + "missingValue": 9999, +} + +# List to be used as fake GID source +FAKE_GID = [0, 1, 2, 3, None] + + +class Test_EUML2GribFileHandler(unittest.TestCase): + """Test the EUML2GribFileHandler reader.""" + + @mock.patch("satpy.readers.eum_l2_grib.ec") + def setUp(self, ec_): + """Set up the test by creating a mocked eccodes library.""" + fake_gid_generator = (i for i in FAKE_GID) + ec_.codes_grib_new_from_file.side_effect = lambda fh: next(fake_gid_generator) + self.ec_ = ec_ + + @unittest.skipIf(sys.platform.startswith("win"), "'eccodes' not supported on Windows") + @mock.patch("satpy.readers.eum_l2_grib.xr") + @mock.patch("satpy.readers.eum_l2_grib.da") + def test_seviri_data_reading(self, da_, xr_): + """Test the reading of data from the product.""" + from satpy.readers.eum_l2_grib import REPEAT_CYCLE_DURATION, EUML2GribFileHandler + from satpy.utils import get_legacy_chunk_size + CHUNK_SIZE = get_legacy_chunk_size() + + with mock.patch("builtins.open", mock.mock_open()) as mock_file: + with mock.patch("satpy.readers.eum_l2_grib.ec", self.ec_): + self.ec_.codes_get_values.return_value = np.ones(1000*1200) + self.ec_.codes_get.side_effect = lambda gid, key: FAKE_SEVIRI_MESSAGE[key] + self.reader = EUML2GribFileHandler( + filename="test.grib", + filename_info={ + "spacecraft": "MET11", + "start_time": datetime.datetime(year=2020, month=10, day=20, + hour=19, minute=45, second=0) + }, + filetype_info={ + "file_type" : "seviri" + } + ) + + dataset_id = make_dataid(name="dummmy", resolution=3000) + + # Checks that the codes_grib_multi_support_on function has been called + self.ec_.codes_grib_multi_support_on.assert_called() + + # Restarts the id generator and clears the call history + fake_gid_generator = (i for i in FAKE_GID) + self.ec_.codes_grib_new_from_file.side_effect = lambda fh: next(fake_gid_generator) + self.ec_.codes_grib_new_from_file.reset_mock() + self.ec_.codes_release.reset_mock() + + # Checks the correct execution of the get_dataset function with a valid parameter_number + valid_dataset = self.reader.get_dataset(dataset_id, {"parameter_number": 30}) + # Checks the correct file open call + mock_file.assert_called_with("test.grib", "rb") + # Checks that the dataset has been created as a DataArray object + assert valid_dataset._extract_mock_name() == "xr.DataArray()" + # Checks that codes_release has been called after each codes_grib_new_from_file call + # (except after the last one which has returned a None) + assert self.ec_.codes_grib_new_from_file.call_count == self.ec_.codes_release.call_count + 1 + + # Restarts the id generator and clears the call history + fake_gid_generator = (i for i in FAKE_GID) + self.ec_.codes_grib_new_from_file.side_effect = lambda fh: next(fake_gid_generator) + self.ec_.codes_grib_new_from_file.reset_mock() + self.ec_.codes_release.reset_mock() + + # Checks the correct execution of the get_dataset function with an invalid parameter_number + invalid_dataset = self.reader.get_dataset(dataset_id, {"parameter_number": 50}) + # Checks that the function returns None + assert invalid_dataset is None + # Checks that codes_release has been called after each codes_grib_new_from_file call + # (except after the last one which has returned a None) + assert self.ec_.codes_grib_new_from_file.call_count == self.ec_.codes_release.call_count + 1 + + # Checks the basic data reading + assert REPEAT_CYCLE_DURATION == 15 + + # Checks the correct execution of the _get_global_attributes and _get_metadata_from_msg functions + attributes = self.reader._get_attributes() + expected_attributes = { + "orbital_parameters": { + "projection_longitude": 9.5 + }, + "sensor": "seviri", + "platform_name": "Meteosat-11" + } + assert attributes == expected_attributes + + # Checks the reading of an array from the message + self.reader._get_xarray_from_msg(0) + + # Checks that dask.array has been called with the correct arguments + name, args, kwargs = da_.mock_calls[0] + assert np.all(args[0] == np.ones((1200, 1000))) + assert args[1] == CHUNK_SIZE + + # Checks that xarray.DataArray has been called with the correct arguments + name, args, kwargs = xr_.mock_calls[0] + assert kwargs["dims"] == ("y", "x") + + # Checks the correct execution of the _get_proj_area function + pdict, area_dict = self.reader._get_proj_area(0) + + expected_pdict = { + "a": 6400000., + "b": 6300000., + "h": 32000000., + "ssp_lon": 9.5, + "nlines": 1000, + "ncols": 1200, + "a_name": "msg_seviri_rss_3km", + "a_desc": "MSG SEVIRI Rapid Scanning Service area definition with 3 km resolution", + "p_id": "", + } + assert pdict == expected_pdict + expected_area_dict = { + "center_point": 500, + "north": 1200, + "east": 1, + "west": 1000, + "south": 1, + } + assert area_dict == expected_area_dict + + # Checks the correct execution of the get_area_def function + with mock.patch("satpy.readers.eum_l2_grib.seviri_calculate_area_extent", + mock.Mock(name="seviri_calculate_area_extent")) as cae: + with mock.patch("satpy.readers.eum_l2_grib.get_area_definition", mock.Mock()) as gad: + dataset_id = make_dataid(name="dummmy", resolution=400.) + self.reader.get_area_def(dataset_id) + # Asserts that seviri_calculate_area_extent has been called with the correct arguments + expected_args = ({"center_point": 500, "east": 1, "west": 1000, "south": 1, "north": 1200, + "column_step": 400., "line_step": 400.},) + name, args, kwargs = cae.mock_calls[0] + assert args == expected_args + # Asserts that get_area_definition has been called with the correct arguments + name, args, kwargs = gad.mock_calls[0] + assert args[0] == expected_pdict + # The second argument must be the return result of seviri_calculate_area_extent + assert args[1]._extract_mock_name() == "seviri_calculate_area_extent()" + + @unittest.skipIf(sys.platform.startswith("win"), "'eccodes' not supported on Windows") + @mock.patch("satpy.readers.eum_l2_grib.xr") + @mock.patch("satpy.readers.eum_l2_grib.da") + def test_fci_data_reading(self, da_, xr_): + """Test the reading of fci data from the product.""" + from satpy.readers.eum_l2_grib import EUML2GribFileHandler + from satpy.utils import get_legacy_chunk_size + CHUNK_SIZE = get_legacy_chunk_size() + + with mock.patch("builtins.open", mock.mock_open()) as mock_file: + with mock.patch("satpy.readers.eum_l2_grib.ec", self.ec_): + self.ec_.codes_get_values.return_value = np.ones(5568*5568) + self.ec_.codes_get.side_effect = lambda gid, key: FAKE_FCI_MESSAGE[key] + self.reader = EUML2GribFileHandler( + filename="test.grib", + filename_info={ + "spacecraft_id": "1", + "start_time": datetime.datetime(year=2020, month=10, day=20, + hour=19, minute=45, second=0) + }, + filetype_info={ + "file_type" : "fci" + } + ) + + dataset_id = make_dataid(name="dummmy", resolution=2000) + + # Checks that the codes_grib_multi_support_on function has been called + self.ec_.codes_grib_multi_support_on.assert_called() + + # Restarts the id generator and clears the call history + fake_gid_generator = (i for i in FAKE_GID) + self.ec_.codes_grib_new_from_file.side_effect = lambda fh: next(fake_gid_generator) + self.ec_.codes_grib_new_from_file.reset_mock() + self.ec_.codes_release.reset_mock() + + # Checks the correct execution of the get_dataset function with a valid parameter_number + valid_dataset = self.reader.get_dataset(dataset_id, {"parameter_number": 30}) + # Checks the correct file open call + mock_file.assert_called_with("test.grib", "rb") + # Checks that the dataset has been created as a DataArray object + assert valid_dataset._extract_mock_name() == "xr.DataArray()" + # Checks that codes_release has been called after each codes_grib_new_from_file call + # (except after the last one which has returned a None) + assert self.ec_.codes_grib_new_from_file.call_count == self.ec_.codes_release.call_count + 1 + + # Restarts the id generator and clears the call history + fake_gid_generator = (i for i in FAKE_GID) + self.ec_.codes_grib_new_from_file.side_effect = lambda fh: next(fake_gid_generator) + self.ec_.codes_grib_new_from_file.reset_mock() + self.ec_.codes_release.reset_mock() + + # Checks the correct execution of the get_dataset function with an invalid parameter_number + invalid_dataset = self.reader.get_dataset(dataset_id, {"parameter_number": 50}) + # Checks that the function returns None + assert invalid_dataset is None + # Checks that codes_release has been called after each codes_grib_new_from_file call + # (except after the last one which has returned a None) + assert self.ec_.codes_grib_new_from_file.call_count == self.ec_.codes_release.call_count + 1 + + # Checks the correct execution of the _get_global_attributes and _get_metadata_from_msg functions + attributes = self.reader._get_attributes() + expected_attributes = { + "orbital_parameters": { + "projection_longitude": 0.0 + }, + "sensor": "fci", + "platform_name": "MTG-i1" + } + assert attributes == expected_attributes + + # Checks the reading of an array from the message + self.reader._get_xarray_from_msg(0) + + # Checks that dask.array has been called with the correct arguments + name, args, kwargs = da_.mock_calls[0] + assert np.all(args[0] == np.ones((5568, 5568))) + assert args[1] == CHUNK_SIZE + + # Checks that xarray.DataArray has been called with the correct arguments + name, args, kwargs = xr_.mock_calls[0] + assert kwargs["dims"] == ("y", "x") + + # Checks the correct execution of the _get_proj_area function + pdict, area_dict = self.reader._get_proj_area(0) + + expected_pdict = { + "a": 6378140000.0, + "b": 6356755000.0, + "h": 35785830098.0, + "ssp_lon": 0.0, + "nlines": 5568, + "ncols": 5568, + "a_name": "msg_fci_fdss_2km", + "a_desc": "MSG FCI Full Disk Scanning Service area definition with 2 km resolution", + "p_id": "" + } + assert pdict == expected_pdict + expected_area_dict = { + "nlines": 5568, + "ncols": 5568 + } + assert area_dict == expected_area_dict + + # Checks the correct execution of the get_area_def function + with mock.patch("satpy.readers.eum_l2_grib.fci_calculate_area_extent", + mock.Mock(name="fci_calculate_area_extent")) as cae: + with mock.patch("satpy.readers.eum_l2_grib.get_area_definition", mock.Mock()) as gad: + dataset_id = make_dataid(name="dummmy", resolution=2000.) + self.reader.get_area_def(dataset_id) + # Asserts that seviri_calculate_area_extent has been called with the correct arguments + expected_args = ({"nlines": 5568, "ncols": 5568, + "column_step": 2000., "line_step": 2000.},) + name, args, kwargs = cae.mock_calls[0] + assert args == expected_args + # Asserts that get_area_definition has been called with the correct arguments + name, args, kwargs = gad.mock_calls[0] + assert args[0] == expected_pdict + # The second argument must be the return result of seviri_calculate_area_extent + assert args[1]._extract_mock_name() == "fci_calculate_area_extent()"