From 969ac2345d6e9908b36b5e9fcf6872c0c297cade Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 27 May 2026 09:24:28 +0200 Subject: [PATCH 1/5] Update `true_color_reproduction` composite to reflect MTG-I1 becoming Meteosat-12 --- satpy/enhancements/ahi.py | 2 +- satpy/etc/readers/msi_eopf.yaml | 270 ++++++++++++++++++++ satpy/readers/msi_eopf.py | 423 ++++++++++++++++++++++++++++++++ 3 files changed, 694 insertions(+), 1 deletion(-) create mode 100644 satpy/etc/readers/msi_eopf.yaml create mode 100644 satpy/readers/msi_eopf.py diff --git a/satpy/enhancements/ahi.py b/satpy/enhancements/ahi.py index 1ae5e3deff..e239adc536 100644 --- a/satpy/enhancements/ahi.py +++ b/satpy/enhancements/ahi.py @@ -75,7 +75,7 @@ def _jma_true_color_reproduction(img_data, platform=None): [-0.0150, 0.8605, 0.1317], [-0.0174, -0.1009, 1.0512]]), - "mtg-i1": np.array([[0.9007, 0.2086, -0.0100], + "meteosat-12": np.array([[0.9007, 0.2086, -0.0100], [-0.0475, 1.0662, -0.0414], [-0.0123, -0.1342, 1.0794]]), diff --git a/satpy/etc/readers/msi_eopf.yaml b/satpy/etc/readers/msi_eopf.yaml new file mode 100644 index 0000000000..f64d07842c --- /dev/null +++ b/satpy/etc/readers/msi_eopf.yaml @@ -0,0 +1,270 @@ +reader: + name: msi_safe + short_name: MSI SAFE L1C + long_name: Sentinel-2 A and B MSI L1C data in SAFE format + description: SAFE Reader for MSI L1C data (Sentinel-2) + status: Nominal + supports_fsspec: false + sensors: [sen2_msi] + default_channels: [] + reader: !!python/name:satpy.readers.core.yaml_reader.FileYAMLReader + +file_types: + l1c_safe_granule: + file_reader: !!python/name:satpy.readers.msi_safe.SAFEMSIL1C + file_patterns: ['{fmission_id:3s}_MSI{process_level:3s}_{observation_time:%Y%m%dT%H%M%S}_N{fprocessing_baseline_number:4d}_R{relative_orbit_number:3d}_T{dtile_number:5s}_{dproduct_discriminator:%Y%m%dT%H%M%S}.SAFE/GRANULE/L1C_T{gtile_number:5s}_A{absolute_orbit_number:6d}_{gfile_discriminator:%Y%m%dT%H%M%S}/IMG_DATA/T{tile_number:5s}_{file_discriminator:%Y%m%dT%H%M%S}_{band_name:3s}.jp2'] + requires: [l1c_safe_metadata, l1c_safe_tile_metadata] + l1c_safe_tile_metadata: + file_reader: !!python/name:satpy.readers.msi_safe.SAFEMSITileMDXML + file_patterns: ['{fmission_id:3s}_MSI{process_level:3s}_{observation_time:%Y%m%dT%H%M%S}_N{fprocessing_baseline_number:4d}_R{relative_orbit_number:3d}_T{dtile_number:5s}_{dproduct_discriminator:%Y%m%dT%H%M%S}.SAFE/GRANULE/L1C_T{gtile_number:5s}_A{absolute_orbit_number:6d}_{gfile_discriminator:%Y%m%dT%H%M%S}/MTD_TL.xml'] + l1c_safe_metadata: + file_reader: !!python/name:satpy.readers.msi_safe.SAFEMSIMDXML + file_patterns: ['{fmission_id:3s}_MSI{process_level:3s}_{observation_time:%Y%m%dT%H%M%S}_N{fprocessing_baseline_number:4d}_R{relative_orbit_number:3d}_T{dtile_number:5s}_{dproduct_discriminator:%Y%m%dT%H%M%S}.SAFE/MTD_MSIL1C.xml'] + +datasets: + B01: + name: B01 + sensor: msi + wavelength: [0.415, 0.443, 0.470] + resolution: 60 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B02: + name: B02 + sensor: msi + wavelength: [0.440, 0.490, 0.540] + resolution: 10 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B03: + name: B03 + sensor: msi + wavelength: [0.540, 0.560, 0.580] + resolution: 10 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B04: + name: B04 + sensor: msi + wavelength: [0.645, 0.665, 0.685] + resolution: 10 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B05: + name: B05 + sensor: msi + wavelength: [0.695, 0.705, 0.715] + resolution: 20 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B06: + name: B06 + sensor: msi + wavelength: [0.731, 0.740, 0.749] + resolution: 20 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B07: + name: B07 + sensor: msi + wavelength: [0.764, 0.783, 0.802] + resolution: 20 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B08: + name: B08 + sensor: msi + wavelength: [0.780, 0.842, 0.905] + resolution: 10 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B8A: + name: B8A + sensor: msi + wavelength: [0.855, 0.865, 0.875] + resolution: 20 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B09: + name: B09 + sensor: msi + wavelength: [0.935, 0.945, 0.955] + resolution: 60 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B10: + name: B10 + sensor: msi + wavelength: [1.365, 1.375, 1.385] + resolution: 60 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B11: + name: B11 + sensor: msi + wavelength: [1.565, 1.610, 1.655] + resolution: 20 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + B12: + name: B12 + sensor: msi + wavelength: [2.100, 2.190, 2.280] + resolution: 20 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + counts: + standard_name: counts + units: "1" + file_type: l1c_safe_granule + + solar_zenith_angle: + name: solar_zenith_angle + resolution: [10, 20, 60] + file_type: l1c_safe_tile_metadata + xml_tag: Sun_Angles_Grid/Zenith + + solar_azimuth_angle: + name: solar_azimuth_angle + resolution: [10, 20, 60] + file_type: l1c_safe_tile_metadata + xml_tag: Sun_Angles_Grid/Azimuth + + satellite_azimuth_angle: + name: satellite_azimuth_angle + resolution: [10, 20, 60] + file_type: l1c_safe_tile_metadata + xml_tag: Viewing_Incidence_Angles_Grids + xml_item: Azimuth + + satellite_zenith_angle: + name: satellite_zenith_angle + resolution: [10, 20, 60] + file_type: l1c_safe_tile_metadata + xml_tag: Viewing_Incidence_Angles_Grids + xml_item: Zenith diff --git a/satpy/readers/msi_eopf.py b/satpy/readers/msi_eopf.py new file mode 100644 index 0000000000..d5e6cf0a76 --- /dev/null +++ b/satpy/readers/msi_eopf.py @@ -0,0 +1,423 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2016-2020 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 . +"""SAFE MSI L1C/L2A reader. + +The MSI data has a special value for saturated pixels. By default, these +pixels are set to np.inf, but for some applications it might be desirable +to have these pixels left untouched. +For this case, the `mask_saturated` flag is available in the reader, and can be +toggled with ``reader_kwargs`` upon Scene creation:: + + scene = satpy.Scene(filenames, + reader='msi_safe', + reader_kwargs={'mask_saturated': False}) + scene.load(['B01']) + +L1C/L2A format description for the files read here: + + https://sentinels.copernicus.eu/documents/247904/685211/S2-PDGS-TAS-DI-PSD-V14.9.pdf/3d3b6c9c-4334-dcc4-3aa7-f7c0deffbaf7?t=1643013091529 + +NOTE: At present, L1B data is not supported. If the user needs radiance data instead of counts or reflectances, these +are retrieved by first calculating the reflectance and then working back to the radiance. L1B radiance data support +will be added once the data is published onto the Copernicus data ecosystem. + +""" + +import logging +from datetime import datetime + +import dask.array as da +import defusedxml.ElementTree as ET +import numpy as np +import xarray as xr +from pyresample import geometry + +from satpy._compat import cached_property +from satpy.readers.core.file_handlers import BaseFileHandler +from satpy.utils import get_legacy_chunk_size + +logger = logging.getLogger(__name__) +CHUNK_SIZE = get_legacy_chunk_size() + +PLATFORMS = {"S2A": "Sentinel-2A", + "S2B": "Sentinel-2B", + "S2C": "Sentinel-2C", + "S2D": "Sentinel-2D"} + + +class SAFEMSIL1C(BaseFileHandler): + """File handler for SAFE MSI files (jp2).""" + + def __init__(self, filename, filename_info, filetype_info, mda, tile_mda, + mask_saturated=True): + """Initialize the reader.""" + super(SAFEMSIL1C, self).__init__(filename, filename_info, + filetype_info) + del mask_saturated + self._channel = filename_info["band_name"] + self.process_level = filename_info["process_level"] + if self.process_level not in ["L1C", "L2A"]: + raise ValueError(f"Unsupported process level: {self.process_level}") + self._tile_mda = tile_mda + self._mda = mda + self.platform_name = PLATFORMS[filename_info["fmission_id"]] + self._start_time = self._tile_mda.start_time() + self._end_time = filename_info["observation_time"] + + def get_dataset(self, key, info): + """Load a dataset.""" + if self._channel != key["name"]: + return + + logger.debug("Reading %s.", key["name"]) + + proj = self._read_from_file(key) + if proj is None: + return + proj.attrs = info.copy() + proj.attrs["platform_name"] = self.platform_name + return proj + + def _read_from_file(self, key): + proj = xr.open_dataset(self.filename, engine="rasterio", chunks=CHUNK_SIZE)["band_data"] + proj = proj.squeeze("band") + if key["calibration"] == "reflectance": + return self._mda.calibrate_to_reflectances(proj, self._channel) + if key["calibration"] == "radiance": + # The calibration procedure differs for L1B and L1C/L2A data! + if self.process_level in ["L1C", "L2A"]: + # For higher level data, radiances must be computed from the reflectance. + # By default, we use the mean solar angles so that the user does not need to resample, + # but the user can also choose to use the solar angles from the tile metadata. + # This is on a coarse grid so for most bands must be resampled before use. + dq = dict(name="solar_zenith_angle", resolution=key["resolution"]) + zen = self._tile_mda.get_dataset(dq, dict(xml_tag="Sun_Angles_Grid/Zenith")) + tmp_refl = self._mda.calibrate_to_reflectances(proj, self._channel) + return self._mda.calibrate_to_radiances(tmp_refl, zen, self._channel) + else: + # For L1B the radiances can be directly computed from the digital counts. + return self._mda.calibrate_to_radiances_l1b(proj, self._channel) + + + if key["calibration"] == "counts": + return self._mda._sanitize_data(proj) + if key["calibration"] in ["aerosol_thickness", "water_vapor"]: + return self._mda.calibrate_to_atmospheric(proj, self._channel) + + @property + def start_time(self): + """Get the start time.""" + return self._start_time + + @property + def end_time(self): + """Get the end time.""" + return self._start_time + + def get_area_def(self, dsid): + """Get the area def.""" + if self._channel != dsid["name"]: + return + + return self._tile_mda.get_area_def(dsid) + + +class SAFEMSIXMLMetadata(BaseFileHandler): + """Base class for SAFE MSI XML metadata filehandlers.""" + + def __init__(self, filename, filename_info, filetype_info, mask_saturated=True): + """Init the reader.""" + super().__init__(filename, filename_info, filetype_info) + self._start_time = filename_info["observation_time"] + self._end_time = filename_info["observation_time"] + self.root = ET.parse(self.filename) + self.tile = filename_info["dtile_number"] + self.process_level = filename_info["process_level"] + self.platform_name = PLATFORMS[filename_info["fmission_id"]] + self.mask_saturated = mask_saturated + import bottleneck # noqa + import geotiepoints # noqa + + @property + def end_time(self): + """Get end time.""" + return self._start_time + + @property + def start_time(self): + """Get start time.""" + return self._start_time + + +class SAFEMSIMDXML(SAFEMSIXMLMetadata): + """File handle for sentinel 2 safe XML generic metadata.""" + + def calibrate_to_reflectances(self, data, band_name): + """Calibrate *data* using the radiometric information for the metadata.""" + quantification = int(self.root.find(".//QUANTIFICATION_VALUE").text) if self.process_level[:2] == "L1" else \ + int(self.root.find(".//BOA_QUANTIFICATION_VALUE").text) + data = self._sanitize_data(data) + return (data + self.band_offset(band_name)) / quantification * 100 + + def calibrate_to_atmospheric(self, data, band_name): + """Calibrate L2A AOT/WVP product.""" + atmospheric_bands = ["AOT", "WVP"] + if self.process_level == "L1C" or self.process_level == "L1B": + return + elif self.process_level == "L2A" and band_name not in atmospheric_bands: + return + + quantification = float(self.root.find(f".//{band_name}_QUANTIFICATION_VALUE").text) + data = self._sanitize_data(data) + return data / quantification + + def _sanitize_data(self, data): + data = data.where(data != self.no_data) + if self.mask_saturated: + data = data.where(data != self.saturated, np.inf) + return data + + def band_offset(self, band): + """Get the band offset for *band*.""" + band_index = self._band_index(band) + return self.band_offsets.get(band_index, 0) + + def _band_index(self, band): + band_indices = self.band_indices + band_conversions = {"B01": "B1", "B02": "B2", "B03": "B3", "B04": "B4", "B05": "B5", "B06": "B6", "B07": "B7", + "B08": "B8", "B8A": "B8A", "B09": "B9", "B10": "B10", "B11": "B11", "B12": "B12"} + band_index = band_indices[band_conversions[band]] + return band_index + + @cached_property + def band_indices(self): + """Get the band indices from the metadata.""" + spectral_info = self.root.findall(".//Spectral_Information") + band_indices = {spec.attrib["physicalBand"]: int(spec.attrib["bandId"]) for spec in spectral_info} + return band_indices + + @cached_property + def band_offsets(self): + """Get the band offsets from the metadata.""" + offsets = self.root.find(".//Radiometric_Offset_List") if self.process_level[:2] == "L1" else \ + self.root.find(".//BOA_ADD_OFFSET_VALUES_LIST") + if offsets is not None: + band_offsets = {int(off.attrib["band_id"]): float(off.text) for off in offsets} + else: + band_offsets = {} + return band_offsets + + def solar_irradiance(self, band_name): + """Get the solar irradiance for a given *band_name*.""" + band_index = self._band_index(band_name) + return self.solar_irradiances[band_index] + + @cached_property + def solar_irradiances(self): + """Get the TOA solar irradiance values from the metadata.""" + irrads = self.root.find(".//Solar_Irradiance_List") + + if irrads is not None: + solar_irrad = {int(irr.attrib["bandId"]): float(irr.text) for irr in irrads} + if len(solar_irrad) > 0: + return solar_irrad + raise ValueError("No solar irradiance values were found in the metadata.") + + @cached_property + def solar_correction_factor(self): + """Get the solar correction factor from the metadata.""" + sed = self.root.find(".//U") + if sed.text is not None: + return float(sed.text) + raise ValueError("Solar correction factor, U, in metadata is missing.") + + @cached_property + def special_values(self): + """Get the special values from the metadata.""" + special_values = self.root.findall(".//Special_Values") + special_values_dict = {value[0].text: float(value[1].text) for value in special_values} + return special_values_dict + + @property + def no_data(self): + """Get the nodata value from the metadata.""" + return self.special_values["NODATA"] + + @property + def saturated(self): + """Get the saturated value from the metadata.""" + return self.special_values["SATURATED"] + + def calibrate_to_radiances_l1b(self, data, band_name): + """Calibrate *data* to radiance using the radiometric information for the metadata.""" + physical_gain = self.physical_gain(band_name) + data = self._sanitize_data(data) + return (data + self.band_offset(band_name)) / physical_gain + + def calibrate_to_radiances(self, data, solar_zenith, band_name): + """Calibrate *data* to radiance using the radiometric information for the metadata. + + This follows the principle set out in the user guide: + https://sentiwiki.copernicus.eu/web/s2-products#S2Products-Level-1CProductsS2-Products-L1Ctrue + """ + ucor = self.solar_correction_factor + solar_irrad_band = self.solar_irradiance(band_name) + + solar_zenith = np.deg2rad(solar_zenith) + + return (data / 100.) * solar_irrad_band * np.cos(solar_zenith) * ucor / np.pi + + def physical_gain(self, band_name): + """Get the physical gain for a given *band_name*.""" + band_index = self._band_index(band_name) + return self.physical_gains[band_index] + + @cached_property + def physical_gains(self): + """Get the physical gains dictionary.""" + physical_gains = {int(elt.attrib["bandId"]): float(elt.text) for elt in self.root.findall(".//PHYSICAL_GAINS")} + return physical_gains + + +def _fill_swath_edges(angles): + """Fill gaps at edges of swath.""" + darr = xr.DataArray(angles, dims=["y", "x"]) + darr = darr.bfill("x") + darr = darr.ffill("x") + darr = darr.bfill("y") + darr = darr.ffill("y") + angles = darr.data + return angles + + +class SAFEMSITileMDXML(SAFEMSIXMLMetadata): + """File handle for sentinel 2 safe XML tile metadata.""" + + def __init__(self, filename, filename_info, filetype_info, mask_saturated=True): + """Init the reader.""" + super().__init__(filename, filename_info, filetype_info, mask_saturated) + self.geocoding = self.root.find(".//Tile_Geocoding") + + def get_area_def(self, dsid): + """Get the area definition of the dataset.""" + area_extent = self._area_extent(dsid["resolution"]) + cols, rows = self._shape(dsid["resolution"]) + area = geometry.AreaDefinition( + self.tile, + "On-the-fly area", + self.tile, + self.projection, + cols, + rows, + area_extent) + return area + + @cached_property + def projection(self): + """Get the geographic projection.""" + from pyproj import CRS + epsg = self.geocoding.find("HORIZONTAL_CS_CODE").text + return CRS(epsg) + + def _area_extent(self, resolution): + cols, rows = self._shape(resolution) + geoposition = self.geocoding.find('Geoposition[@resolution="' + str(resolution) + '"]') + ulx = float(geoposition.find("ULX").text) + uly = float(geoposition.find("ULY").text) + xdim = float(geoposition.find("XDIM").text) + ydim = float(geoposition.find("YDIM").text) + area_extent = (ulx, uly + rows * ydim, ulx + cols * xdim, uly) + return area_extent + + def _shape(self, resolution): + rows = int(self.geocoding.find('Size[@resolution="' + str(resolution) + '"]/NROWS').text) + cols = int(self.geocoding.find('Size[@resolution="' + str(resolution) + '"]/NCOLS').text) + return cols, rows + + def start_time(self): + """Get the observation time from the tile metadata.""" + timestr = self.root.find(".//SENSING_TIME").text + return datetime.strptime(timestr, "%Y-%m-%dT%H:%M:%S.%fZ") + + @staticmethod + def _do_interp(minterp, xcoord, ycoord): + interp_points2 = np.vstack((ycoord.ravel(), xcoord.ravel())) + res = minterp(interp_points2) + return res.reshape(xcoord.shape) + + def interpolate_angles(self, angles, resolution): + """Interpolate the angles.""" + from geotiepoints.multilinear import MultilinearInterpolator + + cols, rows = self._shape(resolution) + + smin = [0, 0] + smax = np.array(angles.shape) - 1 + orders = angles.shape + minterp = MultilinearInterpolator(smin, smax, orders) + minterp.set_values(da.atleast_2d(angles.ravel())) + + y = da.arange(rows, dtype=angles.dtype, chunks=CHUNK_SIZE) / (rows-1) * (angles.shape[0] - 1) + x = da.arange(cols, dtype=angles.dtype, chunks=CHUNK_SIZE) / (cols-1) * (angles.shape[1] - 1) + xcoord, ycoord = da.meshgrid(x, y) + return da.map_blocks(self._do_interp, minterp, xcoord, ycoord, dtype=angles.dtype, chunks=xcoord.chunks) + + def _get_coarse_dataset(self, key, info): + """Get the coarse dataset refered to by `key` from the XML data.""" + angles = self.root.find(".//Tile_Angles") + if key["name"] in ["solar_zenith_angle", "solar_azimuth_angle"]: + angles = self._get_solar_angles(angles, info) + elif key["name"] in ["satellite_zenith_angle", "satellite_azimuth_angle"]: + angles = self._get_satellite_angles(angles, info) + else: + angles = None + return angles + + def _get_solar_angles(self, angles, info): + angles = self._get_values_from_tag(angles, info["xml_tag"]) + return angles + + @staticmethod + def _get_values_from_tag(xml_tree, xml_tag): + elts = xml_tree.findall(xml_tag + "/Values_List/VALUES") + return np.array([[val for val in elt.text.split()] for elt in elts], + dtype=np.float64) + + def _get_satellite_angles(self, angles, info): + arrays = [] + elts = angles.findall(info["xml_tag"] + '[@bandId="1"]') + for elt in elts: + arrays.append(self._get_values_from_tag(elt, info["xml_item"])) + angles = np.nanmean(np.dstack(arrays), axis=-1) + return angles + + def get_dataset(self, key, info): + """Get the dataset referred to by `key`.""" + angles = self._get_coarse_dataset(key, info) + if angles is None: + return None + + angles = _fill_swath_edges(angles) + + res = self.interpolate_angles(angles, key["resolution"]) + + proj = xr.DataArray(res, dims=["y", "x"]) + proj.attrs = info.copy() + proj.attrs["units"] = "degrees" + proj.attrs["platform_name"] = self.platform_name + return proj From 6312c201824e3fea294db562bb5f6f48cdc46d9f Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 27 May 2026 09:32:33 +0200 Subject: [PATCH 2/5] Remove accidentally committed files --- satpy/etc/readers/msi_eopf.yaml | 270 -------------------- satpy/readers/msi_eopf.py | 423 -------------------------------- 2 files changed, 693 deletions(-) delete mode 100644 satpy/etc/readers/msi_eopf.yaml delete mode 100644 satpy/readers/msi_eopf.py diff --git a/satpy/etc/readers/msi_eopf.yaml b/satpy/etc/readers/msi_eopf.yaml deleted file mode 100644 index f64d07842c..0000000000 --- a/satpy/etc/readers/msi_eopf.yaml +++ /dev/null @@ -1,270 +0,0 @@ -reader: - name: msi_safe - short_name: MSI SAFE L1C - long_name: Sentinel-2 A and B MSI L1C data in SAFE format - description: SAFE Reader for MSI L1C data (Sentinel-2) - status: Nominal - supports_fsspec: false - sensors: [sen2_msi] - default_channels: [] - reader: !!python/name:satpy.readers.core.yaml_reader.FileYAMLReader - -file_types: - l1c_safe_granule: - file_reader: !!python/name:satpy.readers.msi_safe.SAFEMSIL1C - file_patterns: ['{fmission_id:3s}_MSI{process_level:3s}_{observation_time:%Y%m%dT%H%M%S}_N{fprocessing_baseline_number:4d}_R{relative_orbit_number:3d}_T{dtile_number:5s}_{dproduct_discriminator:%Y%m%dT%H%M%S}.SAFE/GRANULE/L1C_T{gtile_number:5s}_A{absolute_orbit_number:6d}_{gfile_discriminator:%Y%m%dT%H%M%S}/IMG_DATA/T{tile_number:5s}_{file_discriminator:%Y%m%dT%H%M%S}_{band_name:3s}.jp2'] - requires: [l1c_safe_metadata, l1c_safe_tile_metadata] - l1c_safe_tile_metadata: - file_reader: !!python/name:satpy.readers.msi_safe.SAFEMSITileMDXML - file_patterns: ['{fmission_id:3s}_MSI{process_level:3s}_{observation_time:%Y%m%dT%H%M%S}_N{fprocessing_baseline_number:4d}_R{relative_orbit_number:3d}_T{dtile_number:5s}_{dproduct_discriminator:%Y%m%dT%H%M%S}.SAFE/GRANULE/L1C_T{gtile_number:5s}_A{absolute_orbit_number:6d}_{gfile_discriminator:%Y%m%dT%H%M%S}/MTD_TL.xml'] - l1c_safe_metadata: - file_reader: !!python/name:satpy.readers.msi_safe.SAFEMSIMDXML - file_patterns: ['{fmission_id:3s}_MSI{process_level:3s}_{observation_time:%Y%m%dT%H%M%S}_N{fprocessing_baseline_number:4d}_R{relative_orbit_number:3d}_T{dtile_number:5s}_{dproduct_discriminator:%Y%m%dT%H%M%S}.SAFE/MTD_MSIL1C.xml'] - -datasets: - B01: - name: B01 - sensor: msi - wavelength: [0.415, 0.443, 0.470] - resolution: 60 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B02: - name: B02 - sensor: msi - wavelength: [0.440, 0.490, 0.540] - resolution: 10 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B03: - name: B03 - sensor: msi - wavelength: [0.540, 0.560, 0.580] - resolution: 10 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B04: - name: B04 - sensor: msi - wavelength: [0.645, 0.665, 0.685] - resolution: 10 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B05: - name: B05 - sensor: msi - wavelength: [0.695, 0.705, 0.715] - resolution: 20 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B06: - name: B06 - sensor: msi - wavelength: [0.731, 0.740, 0.749] - resolution: 20 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B07: - name: B07 - sensor: msi - wavelength: [0.764, 0.783, 0.802] - resolution: 20 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B08: - name: B08 - sensor: msi - wavelength: [0.780, 0.842, 0.905] - resolution: 10 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B8A: - name: B8A - sensor: msi - wavelength: [0.855, 0.865, 0.875] - resolution: 20 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B09: - name: B09 - sensor: msi - wavelength: [0.935, 0.945, 0.955] - resolution: 60 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B10: - name: B10 - sensor: msi - wavelength: [1.365, 1.375, 1.385] - resolution: 60 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B11: - name: B11 - sensor: msi - wavelength: [1.565, 1.610, 1.655] - resolution: 20 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - B12: - name: B12 - sensor: msi - wavelength: [2.100, 2.190, 2.280] - resolution: 20 - calibration: - reflectance: - standard_name: toa_bidirectional_reflectance - units: "%" - radiance: - standard_name: toa_outgoing_radiance_per_unit_wavelength - units: W m-2 um-1 sr-1 - counts: - standard_name: counts - units: "1" - file_type: l1c_safe_granule - - solar_zenith_angle: - name: solar_zenith_angle - resolution: [10, 20, 60] - file_type: l1c_safe_tile_metadata - xml_tag: Sun_Angles_Grid/Zenith - - solar_azimuth_angle: - name: solar_azimuth_angle - resolution: [10, 20, 60] - file_type: l1c_safe_tile_metadata - xml_tag: Sun_Angles_Grid/Azimuth - - satellite_azimuth_angle: - name: satellite_azimuth_angle - resolution: [10, 20, 60] - file_type: l1c_safe_tile_metadata - xml_tag: Viewing_Incidence_Angles_Grids - xml_item: Azimuth - - satellite_zenith_angle: - name: satellite_zenith_angle - resolution: [10, 20, 60] - file_type: l1c_safe_tile_metadata - xml_tag: Viewing_Incidence_Angles_Grids - xml_item: Zenith diff --git a/satpy/readers/msi_eopf.py b/satpy/readers/msi_eopf.py deleted file mode 100644 index d5e6cf0a76..0000000000 --- a/satpy/readers/msi_eopf.py +++ /dev/null @@ -1,423 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright (c) 2016-2020 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 . -"""SAFE MSI L1C/L2A reader. - -The MSI data has a special value for saturated pixels. By default, these -pixels are set to np.inf, but for some applications it might be desirable -to have these pixels left untouched. -For this case, the `mask_saturated` flag is available in the reader, and can be -toggled with ``reader_kwargs`` upon Scene creation:: - - scene = satpy.Scene(filenames, - reader='msi_safe', - reader_kwargs={'mask_saturated': False}) - scene.load(['B01']) - -L1C/L2A format description for the files read here: - - https://sentinels.copernicus.eu/documents/247904/685211/S2-PDGS-TAS-DI-PSD-V14.9.pdf/3d3b6c9c-4334-dcc4-3aa7-f7c0deffbaf7?t=1643013091529 - -NOTE: At present, L1B data is not supported. If the user needs radiance data instead of counts or reflectances, these -are retrieved by first calculating the reflectance and then working back to the radiance. L1B radiance data support -will be added once the data is published onto the Copernicus data ecosystem. - -""" - -import logging -from datetime import datetime - -import dask.array as da -import defusedxml.ElementTree as ET -import numpy as np -import xarray as xr -from pyresample import geometry - -from satpy._compat import cached_property -from satpy.readers.core.file_handlers import BaseFileHandler -from satpy.utils import get_legacy_chunk_size - -logger = logging.getLogger(__name__) -CHUNK_SIZE = get_legacy_chunk_size() - -PLATFORMS = {"S2A": "Sentinel-2A", - "S2B": "Sentinel-2B", - "S2C": "Sentinel-2C", - "S2D": "Sentinel-2D"} - - -class SAFEMSIL1C(BaseFileHandler): - """File handler for SAFE MSI files (jp2).""" - - def __init__(self, filename, filename_info, filetype_info, mda, tile_mda, - mask_saturated=True): - """Initialize the reader.""" - super(SAFEMSIL1C, self).__init__(filename, filename_info, - filetype_info) - del mask_saturated - self._channel = filename_info["band_name"] - self.process_level = filename_info["process_level"] - if self.process_level not in ["L1C", "L2A"]: - raise ValueError(f"Unsupported process level: {self.process_level}") - self._tile_mda = tile_mda - self._mda = mda - self.platform_name = PLATFORMS[filename_info["fmission_id"]] - self._start_time = self._tile_mda.start_time() - self._end_time = filename_info["observation_time"] - - def get_dataset(self, key, info): - """Load a dataset.""" - if self._channel != key["name"]: - return - - logger.debug("Reading %s.", key["name"]) - - proj = self._read_from_file(key) - if proj is None: - return - proj.attrs = info.copy() - proj.attrs["platform_name"] = self.platform_name - return proj - - def _read_from_file(self, key): - proj = xr.open_dataset(self.filename, engine="rasterio", chunks=CHUNK_SIZE)["band_data"] - proj = proj.squeeze("band") - if key["calibration"] == "reflectance": - return self._mda.calibrate_to_reflectances(proj, self._channel) - if key["calibration"] == "radiance": - # The calibration procedure differs for L1B and L1C/L2A data! - if self.process_level in ["L1C", "L2A"]: - # For higher level data, radiances must be computed from the reflectance. - # By default, we use the mean solar angles so that the user does not need to resample, - # but the user can also choose to use the solar angles from the tile metadata. - # This is on a coarse grid so for most bands must be resampled before use. - dq = dict(name="solar_zenith_angle", resolution=key["resolution"]) - zen = self._tile_mda.get_dataset(dq, dict(xml_tag="Sun_Angles_Grid/Zenith")) - tmp_refl = self._mda.calibrate_to_reflectances(proj, self._channel) - return self._mda.calibrate_to_radiances(tmp_refl, zen, self._channel) - else: - # For L1B the radiances can be directly computed from the digital counts. - return self._mda.calibrate_to_radiances_l1b(proj, self._channel) - - - if key["calibration"] == "counts": - return self._mda._sanitize_data(proj) - if key["calibration"] in ["aerosol_thickness", "water_vapor"]: - return self._mda.calibrate_to_atmospheric(proj, self._channel) - - @property - def start_time(self): - """Get the start time.""" - return self._start_time - - @property - def end_time(self): - """Get the end time.""" - return self._start_time - - def get_area_def(self, dsid): - """Get the area def.""" - if self._channel != dsid["name"]: - return - - return self._tile_mda.get_area_def(dsid) - - -class SAFEMSIXMLMetadata(BaseFileHandler): - """Base class for SAFE MSI XML metadata filehandlers.""" - - def __init__(self, filename, filename_info, filetype_info, mask_saturated=True): - """Init the reader.""" - super().__init__(filename, filename_info, filetype_info) - self._start_time = filename_info["observation_time"] - self._end_time = filename_info["observation_time"] - self.root = ET.parse(self.filename) - self.tile = filename_info["dtile_number"] - self.process_level = filename_info["process_level"] - self.platform_name = PLATFORMS[filename_info["fmission_id"]] - self.mask_saturated = mask_saturated - import bottleneck # noqa - import geotiepoints # noqa - - @property - def end_time(self): - """Get end time.""" - return self._start_time - - @property - def start_time(self): - """Get start time.""" - return self._start_time - - -class SAFEMSIMDXML(SAFEMSIXMLMetadata): - """File handle for sentinel 2 safe XML generic metadata.""" - - def calibrate_to_reflectances(self, data, band_name): - """Calibrate *data* using the radiometric information for the metadata.""" - quantification = int(self.root.find(".//QUANTIFICATION_VALUE").text) if self.process_level[:2] == "L1" else \ - int(self.root.find(".//BOA_QUANTIFICATION_VALUE").text) - data = self._sanitize_data(data) - return (data + self.band_offset(band_name)) / quantification * 100 - - def calibrate_to_atmospheric(self, data, band_name): - """Calibrate L2A AOT/WVP product.""" - atmospheric_bands = ["AOT", "WVP"] - if self.process_level == "L1C" or self.process_level == "L1B": - return - elif self.process_level == "L2A" and band_name not in atmospheric_bands: - return - - quantification = float(self.root.find(f".//{band_name}_QUANTIFICATION_VALUE").text) - data = self._sanitize_data(data) - return data / quantification - - def _sanitize_data(self, data): - data = data.where(data != self.no_data) - if self.mask_saturated: - data = data.where(data != self.saturated, np.inf) - return data - - def band_offset(self, band): - """Get the band offset for *band*.""" - band_index = self._band_index(band) - return self.band_offsets.get(band_index, 0) - - def _band_index(self, band): - band_indices = self.band_indices - band_conversions = {"B01": "B1", "B02": "B2", "B03": "B3", "B04": "B4", "B05": "B5", "B06": "B6", "B07": "B7", - "B08": "B8", "B8A": "B8A", "B09": "B9", "B10": "B10", "B11": "B11", "B12": "B12"} - band_index = band_indices[band_conversions[band]] - return band_index - - @cached_property - def band_indices(self): - """Get the band indices from the metadata.""" - spectral_info = self.root.findall(".//Spectral_Information") - band_indices = {spec.attrib["physicalBand"]: int(spec.attrib["bandId"]) for spec in spectral_info} - return band_indices - - @cached_property - def band_offsets(self): - """Get the band offsets from the metadata.""" - offsets = self.root.find(".//Radiometric_Offset_List") if self.process_level[:2] == "L1" else \ - self.root.find(".//BOA_ADD_OFFSET_VALUES_LIST") - if offsets is not None: - band_offsets = {int(off.attrib["band_id"]): float(off.text) for off in offsets} - else: - band_offsets = {} - return band_offsets - - def solar_irradiance(self, band_name): - """Get the solar irradiance for a given *band_name*.""" - band_index = self._band_index(band_name) - return self.solar_irradiances[band_index] - - @cached_property - def solar_irradiances(self): - """Get the TOA solar irradiance values from the metadata.""" - irrads = self.root.find(".//Solar_Irradiance_List") - - if irrads is not None: - solar_irrad = {int(irr.attrib["bandId"]): float(irr.text) for irr in irrads} - if len(solar_irrad) > 0: - return solar_irrad - raise ValueError("No solar irradiance values were found in the metadata.") - - @cached_property - def solar_correction_factor(self): - """Get the solar correction factor from the metadata.""" - sed = self.root.find(".//U") - if sed.text is not None: - return float(sed.text) - raise ValueError("Solar correction factor, U, in metadata is missing.") - - @cached_property - def special_values(self): - """Get the special values from the metadata.""" - special_values = self.root.findall(".//Special_Values") - special_values_dict = {value[0].text: float(value[1].text) for value in special_values} - return special_values_dict - - @property - def no_data(self): - """Get the nodata value from the metadata.""" - return self.special_values["NODATA"] - - @property - def saturated(self): - """Get the saturated value from the metadata.""" - return self.special_values["SATURATED"] - - def calibrate_to_radiances_l1b(self, data, band_name): - """Calibrate *data* to radiance using the radiometric information for the metadata.""" - physical_gain = self.physical_gain(band_name) - data = self._sanitize_data(data) - return (data + self.band_offset(band_name)) / physical_gain - - def calibrate_to_radiances(self, data, solar_zenith, band_name): - """Calibrate *data* to radiance using the radiometric information for the metadata. - - This follows the principle set out in the user guide: - https://sentiwiki.copernicus.eu/web/s2-products#S2Products-Level-1CProductsS2-Products-L1Ctrue - """ - ucor = self.solar_correction_factor - solar_irrad_band = self.solar_irradiance(band_name) - - solar_zenith = np.deg2rad(solar_zenith) - - return (data / 100.) * solar_irrad_band * np.cos(solar_zenith) * ucor / np.pi - - def physical_gain(self, band_name): - """Get the physical gain for a given *band_name*.""" - band_index = self._band_index(band_name) - return self.physical_gains[band_index] - - @cached_property - def physical_gains(self): - """Get the physical gains dictionary.""" - physical_gains = {int(elt.attrib["bandId"]): float(elt.text) for elt in self.root.findall(".//PHYSICAL_GAINS")} - return physical_gains - - -def _fill_swath_edges(angles): - """Fill gaps at edges of swath.""" - darr = xr.DataArray(angles, dims=["y", "x"]) - darr = darr.bfill("x") - darr = darr.ffill("x") - darr = darr.bfill("y") - darr = darr.ffill("y") - angles = darr.data - return angles - - -class SAFEMSITileMDXML(SAFEMSIXMLMetadata): - """File handle for sentinel 2 safe XML tile metadata.""" - - def __init__(self, filename, filename_info, filetype_info, mask_saturated=True): - """Init the reader.""" - super().__init__(filename, filename_info, filetype_info, mask_saturated) - self.geocoding = self.root.find(".//Tile_Geocoding") - - def get_area_def(self, dsid): - """Get the area definition of the dataset.""" - area_extent = self._area_extent(dsid["resolution"]) - cols, rows = self._shape(dsid["resolution"]) - area = geometry.AreaDefinition( - self.tile, - "On-the-fly area", - self.tile, - self.projection, - cols, - rows, - area_extent) - return area - - @cached_property - def projection(self): - """Get the geographic projection.""" - from pyproj import CRS - epsg = self.geocoding.find("HORIZONTAL_CS_CODE").text - return CRS(epsg) - - def _area_extent(self, resolution): - cols, rows = self._shape(resolution) - geoposition = self.geocoding.find('Geoposition[@resolution="' + str(resolution) + '"]') - ulx = float(geoposition.find("ULX").text) - uly = float(geoposition.find("ULY").text) - xdim = float(geoposition.find("XDIM").text) - ydim = float(geoposition.find("YDIM").text) - area_extent = (ulx, uly + rows * ydim, ulx + cols * xdim, uly) - return area_extent - - def _shape(self, resolution): - rows = int(self.geocoding.find('Size[@resolution="' + str(resolution) + '"]/NROWS').text) - cols = int(self.geocoding.find('Size[@resolution="' + str(resolution) + '"]/NCOLS').text) - return cols, rows - - def start_time(self): - """Get the observation time from the tile metadata.""" - timestr = self.root.find(".//SENSING_TIME").text - return datetime.strptime(timestr, "%Y-%m-%dT%H:%M:%S.%fZ") - - @staticmethod - def _do_interp(minterp, xcoord, ycoord): - interp_points2 = np.vstack((ycoord.ravel(), xcoord.ravel())) - res = minterp(interp_points2) - return res.reshape(xcoord.shape) - - def interpolate_angles(self, angles, resolution): - """Interpolate the angles.""" - from geotiepoints.multilinear import MultilinearInterpolator - - cols, rows = self._shape(resolution) - - smin = [0, 0] - smax = np.array(angles.shape) - 1 - orders = angles.shape - minterp = MultilinearInterpolator(smin, smax, orders) - minterp.set_values(da.atleast_2d(angles.ravel())) - - y = da.arange(rows, dtype=angles.dtype, chunks=CHUNK_SIZE) / (rows-1) * (angles.shape[0] - 1) - x = da.arange(cols, dtype=angles.dtype, chunks=CHUNK_SIZE) / (cols-1) * (angles.shape[1] - 1) - xcoord, ycoord = da.meshgrid(x, y) - return da.map_blocks(self._do_interp, minterp, xcoord, ycoord, dtype=angles.dtype, chunks=xcoord.chunks) - - def _get_coarse_dataset(self, key, info): - """Get the coarse dataset refered to by `key` from the XML data.""" - angles = self.root.find(".//Tile_Angles") - if key["name"] in ["solar_zenith_angle", "solar_azimuth_angle"]: - angles = self._get_solar_angles(angles, info) - elif key["name"] in ["satellite_zenith_angle", "satellite_azimuth_angle"]: - angles = self._get_satellite_angles(angles, info) - else: - angles = None - return angles - - def _get_solar_angles(self, angles, info): - angles = self._get_values_from_tag(angles, info["xml_tag"]) - return angles - - @staticmethod - def _get_values_from_tag(xml_tree, xml_tag): - elts = xml_tree.findall(xml_tag + "/Values_List/VALUES") - return np.array([[val for val in elt.text.split()] for elt in elts], - dtype=np.float64) - - def _get_satellite_angles(self, angles, info): - arrays = [] - elts = angles.findall(info["xml_tag"] + '[@bandId="1"]') - for elt in elts: - arrays.append(self._get_values_from_tag(elt, info["xml_item"])) - angles = np.nanmean(np.dstack(arrays), axis=-1) - return angles - - def get_dataset(self, key, info): - """Get the dataset referred to by `key`.""" - angles = self._get_coarse_dataset(key, info) - if angles is None: - return None - - angles = _fill_swath_edges(angles) - - res = self.interpolate_angles(angles, key["resolution"]) - - proj = xr.DataArray(res, dims=["y", "x"]) - proj.attrs = info.copy() - proj.attrs["units"] = "degrees" - proj.attrs["platform_name"] = self.platform_name - return proj From 801c00232d11473b53f6c9e0bae0d1731774cbac Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 27 May 2026 21:12:38 +0200 Subject: [PATCH 3/5] Update true color reproduction code to use TSIS-1 v2 solar spectrum. Add Sentinels. --- satpy/enhancements/ahi.py | 73 ++++++++++++++++++++++++--------------- 1 file changed, 45 insertions(+), 28 deletions(-) diff --git a/satpy/enhancements/ahi.py b/satpy/enhancements/ahi.py index e239adc536..b55fb6671b 100644 --- a/satpy/enhancements/ahi.py +++ b/satpy/enhancements/ahi.py @@ -54,34 +54,51 @@ def _jma_true_color_reproduction(img_data, platform=None): """ # Conversion matrix dictionaries specifying sensor and platform. - ccm_dict = {"himawari-8": np.array([[1.1629, 0.1539, -0.2175], - [-0.0252, 0.8725, 0.1300], - [-0.0204, -0.1100, 1.0633]]), - - "himawari-9": np.array([[1.1619, 0.1542, -0.2168], - [-0.0271, 0.8749, 0.1295], - [-0.0202, -0.1103, 1.0634]]), - - "goes-16": np.array([[1.1425, 0.1819, -0.2250], - [-0.0951, 0.9363, 0.1360], - [-0.0113, -0.1179, 1.0621]]), - "goes-17": np.array([[1.1437, 0.1818, -0.2262], - [-0.0952, 0.9354, 0.1371], - [-0.0113, -0.1178, 1.0620]]), - "goes-18": np.array([[1.1629, 0.1539, -0.2175], - [-0.0252, 0.8725, 0.1300], - [-0.0204, -0.1100, 1.0633]]), - "goes-19": np.array([[0.9481, 0.3706, -0.2194], - [-0.0150, 0.8605, 0.1317], - [-0.0174, -0.1009, 1.0512]]), - - "meteosat-12": np.array([[0.9007, 0.2086, -0.0100], - [-0.0475, 1.0662, -0.0414], - [-0.0123, -0.1342, 1.0794]]), - - "geo-kompsat-2a": np.array([[1.1661, 0.1489, -0.2157], - [-0.0255, 0.8745, 0.1282], - [-0.0205, -0.1103, 1.0637]]), + ccm_dict = {"himawari-8": np.array([[1.1495, 0.1534, -0.2203], + [-0.0245, 0.8741, 0.1321], + [-0.0203, -0.1102, 1.0685]]), + + "himawari-9": np.array([[1.1485, 0.1537, -0.2196], + [-0.0264, 0.8765, 0.1316], + [-0.0200, -0.1105, 1.0686]]), + + "goes-16": np.array([[1.1293, 0.1810, -0.2277], + [-0.0940, 0.9377, 0.1380], + [-0.0111, -0.1180, 1.0673]]), + "goes-17": np.array([[1.1306, 0.1809, -0.2289], + [-0.0941, 0.9367, 0.1391], + [-0.0112, -0.1179, 1.0672]]), + "goes-18": np.array([[1.1328, 0.1815, -0.2317], + [-0.0943, 0.9342, 0.1417], + [-0.0112, -0.1176, 1.0669]]), + "goes-19": np.array([[1.1253, 0.1818, -0.2244], + [-0.0937, 0.9403, 0.1351], + [-0.0111, -0.1184, 1.0676]]), + + "meteosat-12": np.array([[0.8860, 0.2081, -0.0115], + [-0.0468, 1.0691, -0.0406], + [-0.0121, -0.1345, 1.0847]]), + + "geo-kompsat-2a": np.array([[1.1527, 0.1484, -0.2185], + [-0.0248, 0.8762, 0.1303], + [-0.0203, -0.1105, 1.0689]]), + + "sentinel-3a": np.array([[2.3828, 0.0046, -1.3048], + [-0.2172, 0.0342, 1.1647], + [-0.0211, -0.0043, 0.9636]]), + "sentinel-3b": np.array([[2.3718, 0.0058, -1.2950], + [-0.2161, 0.0434, 1.1544], + [-0.0211, -0.0055, 0.9646]]), + + "sentinel-2a": np.array([[2.0301, 0.0231, -0.9706], + [-0.1801, 0.1302, 1.0316], + [-0.0186, -0.0164, 0.9732]]), + "sentinel-2b": np.array([[1.9942, 0.0289, -0.9405], + [-0.1774, 0.1617, 0.9974], + [-0.0182, -0.0204, 0.9767]]), + "sentinel-2c": np.array([[1.7759, 0.0666, -0.7598], + [-0.1589, 0.3731, 0.7675], + [-0.0161, -0.0470, 1.0013]]), } # A conversion matrix, sensor name and platform name is required From f79ec2266f04efbf5bd919e491b5f176b889631c Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Thu, 28 May 2026 10:08:43 +0200 Subject: [PATCH 4/5] Update tests for true color reproduction --- satpy/tests/enhancement_tests/test_ahi.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/satpy/tests/enhancement_tests/test_ahi.py b/satpy/tests/enhancement_tests/test_ahi.py index ec6e0df144..37fd6e6ec4 100644 --- a/satpy/tests/enhancement_tests/test_ahi.py +++ b/satpy/tests/enhancement_tests/test_ahi.py @@ -39,14 +39,14 @@ def test_jma_true_color_reproduction(self): from satpy.enhancements.ahi import jma_true_color_reproduction - expected = [[[-109.93, 10.993, 131.916, 252.839, 373.762], - [494.685, 615.608, 736.531, 857.454, 978.377]], + expected = [[[-108.260, 10.826, 129.912, 248.998, 368.084], + [487.170, 606.256, 725.342, 844.428, 963.514]], - [[-97.73, 9.773, 117.276, 224.779, 332.282], - [439.785, 547.288, 654.791, 762.294, 869.797]], + [[ -98.170, 9.817, 117.804, 225.791, 333.778], + [441.765, 549.752, 657.739, 765.726, 873.713]], - [[-93.29, 9.329, 111.948, 214.567, 317.186], - [419.805, 522.424, 625.043, 727.662, 830.281]]] + [[-93.800, 9.380, 112.560, 215.740, 318.920], + [422.100, 525.280, 628.460, 731.640, 834.820]]] img = XRImage(self.rgb) jma_true_color_reproduction(img) From 78937529a26edf2bf83e5dbf94f953b4f3450a41 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Thu, 28 May 2026 10:17:19 +0200 Subject: [PATCH 5/5] Update true color reproduction spaces --- satpy/enhancements/ahi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/enhancements/ahi.py b/satpy/enhancements/ahi.py index b55fb6671b..3f957b6b88 100644 --- a/satpy/enhancements/ahi.py +++ b/satpy/enhancements/ahi.py @@ -58,7 +58,7 @@ def _jma_true_color_reproduction(img_data, platform=None): [-0.0245, 0.8741, 0.1321], [-0.0203, -0.1102, 1.0685]]), - "himawari-9": np.array([[1.1485, 0.1537, -0.2196], + "himawari-9": np.array([[1.1485, 0.1537, -0.2196], [-0.0264, 0.8765, 0.1316], [-0.0200, -0.1105, 1.0686]]),