diff --git a/satpy/etc/composites/ec_msi.yaml b/satpy/etc/composites/ec_msi.yaml new file mode 100644 index 0000000000..13a00898d8 --- /dev/null +++ b/satpy/etc/composites/ec_msi.yaml @@ -0,0 +1,39 @@ +sensor_name: visir/ec_msi + + +modifiers: + sunz_corrected: + modifier: !!python/name:satpy.modifiers.SunZenithCorrector + + rayleigh_corrected: + modifier: !!python/name:satpy.modifiers.PSPRayleighReflectance + atmosphere: us-standard + aerosol_type: rayleigh_only + prerequisites: + - name: VIS + modifiers: [sunz_corrected] + optional_prerequisites: + - satellite_azimuth_angle + - satellite_zenith_angle + - solar_azimuth_angle + - solar_zenith_angle + +composites: + natural_color_nocorr: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - SWIR1 + - NIR + - VIS + standard_name: natural_color + + natural_color: + compositor: !!python/name:satpy.composites.GenericCompositor + prerequisites: + - name: SWIR1 + modifiers: [sunz_corrected] + - name: NIR + modifiers: [sunz_corrected] + - name: VIS + modifiers: [sunz_corrected] + standard_name: natural_color diff --git a/satpy/etc/composites/msi.yaml b/satpy/etc/composites/sen2_msi.yaml similarity index 99% rename from satpy/etc/composites/msi.yaml rename to satpy/etc/composites/sen2_msi.yaml index 010bd240b0..b600db5d95 100644 --- a/satpy/etc/composites/msi.yaml +++ b/satpy/etc/composites/sen2_msi.yaml @@ -1,4 +1,4 @@ -sensor_name: visir/msi +sensor_name: visir/sen2_msi modifiers: diff --git a/satpy/etc/readers/msi_l1c_earthcare.yaml b/satpy/etc/readers/msi_l1c_earthcare.yaml new file mode 100644 index 0000000000..c25d5abd64 --- /dev/null +++ b/satpy/etc/readers/msi_l1c_earthcare.yaml @@ -0,0 +1,262 @@ +reader: + name: msi_l1c_earthcare + short_name: MSI EarthCARE + long_name: Multispectral Imager for EarthCARE + description: Multispectral Imager for EarthCARE Level 1C (regridded) Reader + status: Nominal + supports_fsspec: true + sensors: [ec_msi] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + + +file_types: + msi_l1c_earthcare_rgr: + file_reader: !!python/name:satpy.readers.msi_ec_l1c_h5.MSIECL1CFileHandler + file_patterns: + - '{mission_id:s}_{processing_institute:s}_{sensor_id:s}_{file_id:s}_{proc_level:s}_{start_time:%Y%m%dT%H%M%S}Z_{end_time:%Y%m%dT%H%M%S}Z_{orbit_number:s}{frame_id:s}.h5' + + +datasets: +# Science measurement datasets + VIS: + name: VIS + sensor: ec_msi + wavelength: [0.66, 0.67, 0.68] + resolution: 500 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance + units: W m-2 sr-1 + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 0 + coordinates: [longitude, latitude] + NIR: + name: NIR + sensor: ec_msi + wavelength: [0.855, 0.865, 0.875] + resolution: 500 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance + units: W m-2 sr-1 + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 1 + coordinates: [longitude, latitude] + SWIR1: + name: SWIR1 + sensor: ec_msi + wavelength: [1.64, 1.67, 1.70] + resolution: 500 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance + units: W m-2 sr-1 + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 2 + coordinates: [longitude, latitude] + SWIR2: + name: SWIR2 + sensor: ec_msi + wavelength: [2.16, 2.21, 2.26] + resolution: 500 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance + units: W m-2 sr-1 + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 3 + coordinates: [longitude, latitude] + TIR1: + name: TIR1 + sensor: ec_msi + wavelength: [8.35, 8.80, 9.25] + resolution: 500 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 4 + coordinates: [longitude, latitude] + TIR2: + name: TIR2 + sensor: ec_msi + wavelength: [10.35, 10.80, 11.25] + resolution: 500 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 5 + coordinates: [longitude, latitude] + TIR3: + name: TIR3 + sensor: ec_msi + wavelength: [11.55,12.00,12.45] + resolution: 500 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values + band_index: 6 + coordinates: [longitude, latitude] + + # Relative error datasets + VIS_rel_error: + name: VIS_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_radiance + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 0 + NIR_rel_error: + name: NIR_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_radiance + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 1 + SWIR1_rel_error: + name: SWIR1_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_radiance + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 2 + SWIR2_rel_error: + name: SWIR2_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_radiance + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 3 + TIR1_rel_error: + name: TIR1_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_brightness_temperature + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 4 + TIR2_rel_error: + name: TIR2_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_brightness_temperature + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 5 + TIR3_rel_error: + name: TIR3_rel_error + sensor: ec_msi + resolution: 500 + standard_name: relative_error_in_toa_brightness_temperature + units: "%" + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/pixel_values_relative_error + band_index: 6 + +# Geolocation data + longitude: + name: longitude + units: degrees_east + standard_name: longitude + resolution: 500 + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/longitude + latitude: + name: latitude + units: degrees_north + standard_name: latitude + resolution: 500 + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/latitude + solar_azimuth_angle: + name: solar_azimuth_angle + units: degrees + standard_name: solar_azimuth_angle + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/solar_azimuth_angle + sensor_azimuth_angle: + name: sensor_azimuth_angle + units: degrees + standard_name: sensor_azimuth_angle + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/sensor_azimuth_angle + sensor_view_angle: + name: sensor_zenith_angle + units: degrees + standard_name: sensor_zenith_angle + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: NonStandard/sensor_view_angle + solar_zenith_angle: + name: solar_zenith_angle + units: degrees + standard_name: solar_zenith_angle + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: NonStandard/solar_zenith_angle + +# Ancillary data + land_flag: + name: land_water_mask + units: 1 + standard_name: land_water_mask + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/land_flag + surface_elevation: + name: surface_elevation + units: m + standard_name: surface_elevation + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: ScienceData/surface_elevation + surface_index: + name: surface_index + units: 1 + standard_name: surface_index + resolution: 500 + coordinates: [longitude, latitude] + file_type: msi_l1c_earthcare_rgr + file_key: NonStandard/surface_index diff --git a/satpy/etc/readers/msi_safe.yaml b/satpy/etc/readers/msi_safe.yaml index 20b324a036..b04fa8623b 100644 --- a/satpy/etc/readers/msi_safe.yaml +++ b/satpy/etc/readers/msi_safe.yaml @@ -5,7 +5,7 @@ reader: description: SAFE Reader for MSI data (Sentinel-2) status: Nominal supports_fsspec: false - sensors: [msi] + sensors: [sen2_msi] default_channels: [] reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader diff --git a/satpy/readers/msi_ec_l1c_h5.py b/satpy/readers/msi_ec_l1c_h5.py new file mode 100644 index 0000000000..c5d2a0c103 --- /dev/null +++ b/satpy/readers/msi_ec_l1c_h5.py @@ -0,0 +1,102 @@ +#!/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 . +"""A reader for Level 1C data produced by the MSI instrument aboard EarthCARE.""" +import logging + +import numpy as np + +from satpy.readers.hdf5_utils import HDF5FileHandler +from satpy.utils import get_legacy_chunk_size + +LOG = logging.getLogger(__name__) +CHUNK_SIZE = get_legacy_chunk_size() + + +class MSIECL1CFileHandler(HDF5FileHandler): + """File handler for MSI L1c H5 files.""" + + def __init__(self, filename, filename_info, filetype_info): + """Init the file handler.""" + super(MSIECL1CFileHandler, self).__init__(filename, + filename_info, + filetype_info) + + @property + def end_time(self): + """Get end time.""" + return self.filename_info["end_time"] + + @property + def start_time(self): + """Get start time.""" + return self.filename_info["start_time"] + + def get_dataset(self, dataset_id, ds_info): + """Load data variable and metadata and calibrate if needed.""" + file_key = ds_info.get("file_key", dataset_id["name"]) + data = self[file_key] + + # Band data is stored in a 3d array (Band x Along_Track x Across_Track). + # This means we have to select a single 2d array for a given band, + # and the correct index is given in the reader YAML. + band_index = ds_info.get("band_index") + if band_index is not None: + data = data[band_index] + + # The dataset has incorrect units attribute (due to storing multiple types). Fix it here. + data.attrs.update(ds_info) + data.attrs.update({"units": ds_info.get("units")}) + # VIS/SWIR data can have radiance or reflectance calibration. + if "calibration" in ds_info: + cal_type = ds_info["calibration"].name + data = self._calibrate(data, band_index, cal_type) + + # Rename dimensions, as some have incorrect names (notably the pixel value data). + if "dim_1" in data.dims: + data = data.rename({"dim_1": "y", "dim_2": "x"}) + + # The dimension list is usually a reference to an H5 variable, which is problematic + # when making a copy of the data. This sorts out the dimensions and sets them correctly + # following the process done in the OMPS reader. + if "DIMENSION_LIST" in data.attrs: + data = self._fix_dims(data, file_key) + + return data + + def _fix_dims(self, data, file_key): + """The pixel data has badly named coordinates, this fixes them.""" + data.attrs.pop("DIMENSION_LIST") + dimensions = self.get_reference(file_key, "DIMENSION_LIST") + dim_dict = {} + # We have to loop over dimensions to match dim sizes as the pixel data is 3d rather than 2d. + for i in range(0, len(data.dims)): + c_dim = data.dims[i] + for r_dim in dimensions: + if data.shape[i] == r_dim[0].shape[0]: + dim_dict[c_dim] = r_dim[0] + data.assign_coords(dim_dict) + return data + + def _calibrate(self, chan_data, band_index, cal_type): + """Calibrate the data.""" + if cal_type == "reflectance": + sol_irrad = self["NonStandard/solar_irradiance"] + chan_data.data = chan_data.data * 100. * np.pi / float(sol_irrad[band_index]) + return chan_data + + return chan_data diff --git a/satpy/tests/reader_tests/test_msi_ec_l1c.py b/satpy/tests/reader_tests/test_msi_ec_l1c.py new file mode 100644 index 0000000000..79e53f740b --- /dev/null +++ b/satpy/tests/reader_tests/test_msi_ec_l1c.py @@ -0,0 +1,162 @@ +"""Tests for the EarthCARE MSI L1c reader.""" +import os +from unittest import mock + +import dask.array as da +import numpy as np +import pytest +import xarray as xr + +from satpy.tests.reader_tests.test_hdf5_utils import FakeHDF5FileHandler + +N_BANDS = 7 +N_SCANS = 20 +N_COLS = 2048 +SHAPE_SC = (300, 6000) +SOL_IRRAD = [30.9, 19.59, 14.77, 8.25, 0., 0., 0.] +DIMLIST = np.ones((N_BANDS, N_SCANS, N_COLS)) + + +class FakeHDF5FileHandler2(FakeHDF5FileHandler): + """Swap-in HDF5 File Handler.""" + + + def _setup_test_data(self, N_BANDS, N_SCANS, N_COLS): + # Set some default attributes + data = { + "ScienceData/pixel_values": + xr.DataArray( + da.ones((N_BANDS, N_SCANS, N_COLS), chunks=1024, dtype=np.float32), + attrs={"units": "Wm-2 sr-1 or K", "DIMENSION_LIST": DIMLIST}, + dims=("band", "dim_2", "dim_1")), + "ScienceData/land_flag": + xr.DataArray( + da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.uint16), + attrs={"units": ""}, + dims=("along_track", "across_track")), + "ScienceData/solar_azimuth_angle": + xr.DataArray( + da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.float32), + attrs={"units": "degrees"}, + dims=("along_track", "across_track")), + "ScienceData/longitude": + xr.DataArray( + da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.float32), + attrs={"units": "degrees"}, + dims=("along_track", "across_track")), + "ScienceData/latitude": + xr.DataArray( + da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.float32), + attrs={"units": "degrees"}, + dims=("along_track", "across_track")), + "NonStandard/solar_irradiance": + xr.DataArray( + da.array(SOL_IRRAD), + attrs={"units": "W m-2"}, + dims=("band")), + } + + return data + + def get_test_content(self, filename, filename_info, filetype_info): + """Mimic reader input file content.""" + test_content = self._setup_test_data(N_BANDS, N_SCANS, N_COLS) + return test_content + + +class ECMSIL1CTester: + """Test MSI/EarthCARE L1C Reader.""" + + def setup_method(self): + """Wrap HDF5 file handler with our own fake handler.""" + from satpy._config import config_search_paths + from satpy.readers.msi_ec_l1c_h5 import MSIECL1CFileHandler + self.reader_configs = config_search_paths(os.path.join("readers", self.yaml_file)) + # http://stackoverflow.com/questions/12219967/how-to-mock-a-base-class-with-python-mock-library + self.p = mock.patch.object(MSIECL1CFileHandler, "__bases__", (FakeHDF5FileHandler2,)) + self.fake_handler = self.p.start() + self.p.is_local = True + + def teardown_method(self): + """Stop wrapping the HDF5 file handler.""" + self.p.stop() + + +class TestECMSIL1C(ECMSIL1CTester): + """Test the EarthCARE MSI L1C reader.""" + + yaml_file = "msi_l1c_earthcare.yaml" + filename = "ECA_EXAA_MSI_RGR_1C_20250625T005649Z_20250625T024013Z_42043E.h5" + + + @mock.patch("satpy.readers.hdf5_utils.HDF5FileHandler._get_reference") + @mock.patch("h5py.File") + def test_get_pixvalues(self, mock_h5py_file, mock_hdf5_utils_get_reference): + """Test loadingpixel values from file.""" + from satpy.readers import load_reader + mock_h5py_file.return_value = mock.MagicMock() + mock_hdf5_utils_get_reference.return_value = DIMLIST + reader = load_reader(self.reader_configs) + files = reader.select_files_from_pathnames([self.filename]) + assert 1 == len(files) + reader.create_filehandlers(files) + # Make sure we have some files + assert reader.file_handlers + available_datasets = list(reader.available_dataset_ids) + assert len(available_datasets) == 27 + + res = reader.load(["VIS", "NIR", "TIR1", "TIR3", "solar_azimuth_angle", "land_water_mask"]) + assert len(res) == 6 + with pytest.raises(KeyError): + res["TIR2"] + with pytest.raises(KeyError): + res["SWIR1"] + + assert res["VIS"].shape == (20, 2048) + assert res["VIS"].attrs["calibration"] == "reflectance" + assert res["VIS"].attrs["units"] == "%" + + assert res["TIR1"].shape == (20, 2048) + assert res["TIR1"].attrs["calibration"] == "brightness_temperature" + assert res["TIR1"].attrs["units"] == "K" + assert res["TIR1"].dtype == np.float32 + + assert res["solar_azimuth_angle"].shape == (20, 2048) + assert res["solar_azimuth_angle"].attrs["units"] == "degrees" + assert res["solar_azimuth_angle"].dtype == np.float32 + + assert res["land_water_mask"].shape == (20, 2048) + assert res["land_water_mask"].attrs["units"] == 1 + assert res["land_water_mask"].dtype == np.uint16 + + + @mock.patch("satpy.readers.hdf5_utils.HDF5FileHandler._get_reference") + @mock.patch("h5py.File") + def test_calibration(self, mock_h5py_file, mock_hdf5_utils_get_reference): + """Test loadingpixel values from file.""" + from satpy.readers import load_reader + from satpy.tests.utils import make_dataid + + mock_h5py_file.return_value = mock.MagicMock() + mock_hdf5_utils_get_reference.return_value = DIMLIST + + reader = load_reader(self.reader_configs) + files = reader.select_files_from_pathnames([self.filename]) + reader.create_filehandlers(files) + + with pytest.raises(KeyError): + reader.load([make_dataid(name="VIS", calibration="counts")]) + with pytest.raises(KeyError): + reader.load([make_dataid(name="TIR1", calibration="counts")]) + with pytest.raises(KeyError): + reader.load([make_dataid(name="TIR1", calibration="radiance")]) + + res = reader.load([make_dataid(name="VIS", calibration="radiance")]) + assert res["VIS"].attrs["calibration"] == "radiance" + assert res["VIS"].attrs["units"] == "W m-2 sr-1" + assert np.all(np.array(res["VIS"].data) == 1) + + res = reader.load([make_dataid(name="NIR", calibration="reflectance")]) + assert res["NIR"].attrs["calibration"] == "reflectance" + assert res["NIR"].attrs["units"] == "%" + assert np.all(np.array(res["NIR"].data) == 1 * np.pi * 100 / SOL_IRRAD[1])