From ec149369a5a9b018cff5aa6fab40307dfd4b9265 Mon Sep 17 00:00:00 2001 From: BRAUN REMI Date: Thu, 7 Apr 2022 16:36:30 +0200 Subject: [PATCH] FIX: Loading every optical band in reflectance (fixed for `Maxar`, `Planet` and `Vision-1` data) (#30) --- CHANGES.md | 1 + eoreader/products/optical/dimap_product.py | 21 +- eoreader/products/optical/landsat_product.py | 21 +- eoreader/products/optical/maxar_product.py | 252 +++++++++++++++++- eoreader/products/optical/optical_product.py | 45 +++- eoreader/products/optical/pla_product.py | 26 +- eoreader/products/optical/s2_product.py | 20 +- eoreader/products/optical/s2_theia_product.py | 20 +- eoreader/products/optical/s3_olci_product.py | 6 - eoreader/products/optical/s3_product.py | 28 +- eoreader/products/optical/s3_slstr_product.py | 6 - eoreader/products/optical/vis1_product.py | 67 ++++- 12 files changed, 476 insertions(+), 37 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 2db6bdc6..d6ccc656 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -16,6 +16,7 @@ ### Bug Fixes +- FIX: Loading every optical band in reflectance (fixed for `Maxar`, `Planet` and `Vision-1` data) ([#30](https://github.com/sertit/eoreader/issues/30)) - FIX: Fixing `ReferenceError: weakly-referenced object no longer exists` when deleting an object - FIX: Do not set sea values to nodata when orthorectifying SAR data with SNAP - FIX: Handle `Sentinel-2` data with PB < 02.07 as `L2Ap` products diff --git a/eoreader/products/optical/dimap_product.py b/eoreader/products/optical/dimap_product.py index 0b6e1142..2401edbe 100644 --- a/eoreader/products/optical/dimap_product.py +++ b/eoreader/products/optical/dimap_product.py @@ -31,6 +31,7 @@ import geopandas as gpd import numpy as np import rasterio +import xarray as xr from cloudpathlib import CloudPath from lxml import etree from rasterio import crs as riocrs @@ -508,8 +509,24 @@ def _manage_invalid_pixels( return self._set_nodata_mask(band_arr, nodata) def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ # Get MTD XML file root, _ = self.read_mtd() rad_proc = DimapRadiometricProcessing.from_value( diff --git a/eoreader/products/optical/landsat_product.py b/eoreader/products/optical/landsat_product.py index f9bc9bb9..91e0d28c 100644 --- a/eoreader/products/optical/landsat_product.py +++ b/eoreader/products/optical/landsat_product.py @@ -26,6 +26,7 @@ import geopandas as gpd import numpy as np import pandas as pd +import xarray as xr from cloudpathlib import CloudPath from lxml import etree from lxml.builder import E @@ -520,8 +521,24 @@ def _read_mtd(self, force_pd=False) -> (etree._Element, dict): return mtd_data def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ # Get band name: the last number of the filename: # ie: 'LC08_L1TP_200030_20191218_20191226_01_T1_B1' if self.is_archived: diff --git a/eoreader/products/optical/maxar_product.py b/eoreader/products/optical/maxar_product.py index f795fc84..bcf25559 100644 --- a/eoreader/products/optical/maxar_product.py +++ b/eoreader/products/optical/maxar_product.py @@ -15,18 +15,20 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -Maxar super class. +Maxar sensors (GeoEye, WorldViews...) class. See `here `_ for more information. """ import logging from abc import abstractmethod +from collections import namedtuple from datetime import datetime from enum import unique from pathlib import Path from typing import Union import numpy as np +import xarray as xr from cloudpathlib import CloudPath from lxml import etree from rasterio import crs as riocrs @@ -43,6 +45,152 @@ LOGGER = logging.getLogger(EOREADER_NAME) +MAXAR_BAND_MTD = { + obn.NIR: "N", + obn.NARROW_NIR: "N", + obn.RED: "R", + obn.GREEN: "G", + obn.BLUE: "B", + obn.WV: "N2", + obn.VRE_1: "RE", + obn.VRE_2: "RE", + obn.VRE_3: "RE", + obn.YELLOW: "Y", + obn.CA: "C", + obn.PAN: "P", +} + +GainOffset = namedtuple("GainOffset", ["gain", "offset"], defaults=[1.0, 0.0]) +MAXAR_GAIN_OFFSET = { + Platform.GE01: { + obn.PAN: GainOffset(gain=1.001, offset=0.0), + obn.BLUE: GainOffset(gain=1.041, offset=0.0), + obn.GREEN: GainOffset(gain=0.972, offset=0.0), + obn.RED: GainOffset(gain=0.979, offset=0.0), + obn.NIR: GainOffset(gain=0.951, offset=0.0), + obn.NARROW_NIR: GainOffset(gain=0.951, offset=0.0), + }, # 2018v0 + Platform.WV02: { + obn.PAN: GainOffset(gain=0.949, offset=-5.523), + obn.CA: GainOffset(gain=1.203, offset=-11.839), + obn.BLUE: GainOffset(gain=1.002, offset=-9.835), + obn.GREEN: GainOffset(gain=0.953, offset=-7.218), + obn.YELLOW: GainOffset(gain=0.946, offset=-5.675), + obn.RED: GainOffset(gain=0.955, offset=-5.046), + obn.VRE_1: GainOffset(gain=0.980, offset=-6.114), + obn.VRE_2: GainOffset(gain=0.980, offset=-6.114), + obn.VRE_3: GainOffset(gain=0.980, offset=-6.114), + obn.NIR: GainOffset(gain=0.966, offset=-5.096), + obn.NARROW_NIR: GainOffset(gain=0.966, offset=-5.096), + obn.WV: GainOffset(gain=1.01, offset=-4.059), + }, # 2018v0 + Platform.WV03: { + obn.PAN: GainOffset(gain=0.955, offset=-5.505), + obn.CA: GainOffset(gain=0.938, offset=-13.099), + obn.BLUE: GainOffset(gain=0.946, offset=-9.409), + obn.GREEN: GainOffset(gain=0.958, offset=-7.771), + obn.YELLOW: GainOffset(gain=0.979, offset=-5.489), + obn.RED: GainOffset(gain=0.969, offset=-4.579), + obn.VRE_1: GainOffset(gain=1.027, offset=-5.552), + obn.VRE_2: GainOffset(gain=1.027, offset=-5.552), + obn.VRE_3: GainOffset(gain=1.027, offset=-5.552), + obn.NIR: GainOffset(gain=0.977, offset=-6.508), + obn.NARROW_NIR: GainOffset(gain=0.977, offset=-6.508), + obn.WV: GainOffset(gain=1.007, offset=-3.699), + }, # 2018v0 + Platform.WV04: { + obn.PAN: GainOffset(gain=1.0, offset=0.0), + obn.BLUE: GainOffset(gain=1.0, offset=0.0), + obn.GREEN: GainOffset(gain=1.0, offset=0.0), + obn.RED: GainOffset(gain=1.0, offset=0.0), + obn.NIR: GainOffset(gain=1.0, offset=0.0), + obn.NARROW_NIR: GainOffset(gain=1.0, offset=0.0), + }, # 2017v0 + Platform.QB: { + obn.PAN: GainOffset(gain=0.870, offset=-1.491), + obn.BLUE: GainOffset(gain=1.105, offset=-2.820), + obn.GREEN: GainOffset(gain=1.071, offset=-3.338), + obn.RED: GainOffset(gain=1.060, offset=-2.954), + obn.NIR: GainOffset(gain=1.020, offset=-4.722), + obn.NARROW_NIR: GainOffset(gain=1.020, offset=-4.722), + }, # 2016v0.Int + Platform.WV01: { + obn.PAN: GainOffset(gain=1.016, offset=-1.824), + }, # 2016v0.Int +} +""" +The TDI specific abscalfactor and effectiveBandwidth are delivered with the imagery in the metadata file. The +digital number, DN, is the pixel value found in the imagery. The Gain and Offset are the absolute radiometric +calibration band dependent adjustment factors that are given in Table 1. Note that these are not necessarily +stagnant values and they are revisited annually. + +Only using last calibration as the Maxar sensors have been found to be very stable throughout their lifetime. + +See `here _` for the values. +""" + +MAXAR_E0 = { + Platform.GE01: { + obn.PAN: 1610.73, + obn.BLUE: 1993.18, + obn.GREEN: 1828.83, + obn.RED: 1491.49, + obn.NIR: 1022.58, + obn.NARROW_NIR: 1022.58, + }, + Platform.WV02: { + obn.PAN: 1571.36, + obn.CA: 1773.81, + obn.BLUE: 2007.27, + obn.GREEN: 1829.62, + obn.YELLOW: 1701.85, + obn.RED: 1538.85, + obn.VRE_1: 1346.09, + obn.VRE_2: 1346.09, + obn.VRE_3: 1346.09, + obn.NIR: 1053.21, + obn.NARROW_NIR: 1053.21, + obn.WV: 856.599, + }, + Platform.WV03: { + obn.PAN: 1574.41, + obn.CA: 1757.89, + obn.BLUE: 2004.61, + obn.GREEN: 1830.18, + obn.YELLOW: 1712.07, + obn.RED: 1535.33, + obn.VRE_1: 1348.08, + obn.VRE_2: 1348.08, + obn.VRE_3: 1348.08, + obn.NIR: 1055.94, + obn.NARROW_NIR: 1055.94, + obn.WV: 858.77, + }, + Platform.WV04: { + obn.PAN: 1608.01, + obn.BLUE: 2009.45, + obn.GREEN: 1831.88, + obn.RED: 1492.12, + obn.NIR: 937.8, + obn.NARROW_NIR: 937.8, + }, + Platform.QB: { + obn.PAN: 1370.92, + obn.BLUE: 1949.59, + obn.GREEN: 1823.64, + obn.RED: 1553.78, + obn.NIR: 1102.85, + obn.NARROW_NIR: 1102.85, + }, + Platform.WV01: { + obn.PAN: 1478.62, + }, +} +""" +Esun is the band-averaged Solar exoatmospheric irradiance @1AU (see Slide 7&8). DG calibration team uses the Thuillier 2003 solar curve for their calculations. +See `here _` for the values. +""" + @unique class MaxarProductType(ListEnum): @@ -531,11 +679,31 @@ def _read_mtd(self) -> (etree._Element, dict): return self._read_mtd_xml(mtd_from_path, mtd_archived) def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ # Delivered in uint16 - # TODO: Convert DN into radiance ! - # TODO: Convert radiance into reflectance ! + + # Convert DN into radiance + band_arr = self._dn_to_toa_rad(band_arr, band) + + # Convert radiance into reflectance + band_arr = self._toa_rad_to_toa_refl(band_arr, band) # To float32 if band_arr.dtype != np.float32: @@ -580,3 +748,77 @@ def _get_tile_path(self) -> Union[CloudPath, Path]: Union[CloudPath, Path]: DIMAP filepath """ return self._get_path(extension="TIL") + + def _dn_to_toa_rad(self, dn_arr: xr.DataArray, band: BandNames) -> xr.DataArray: + """ + Compute DN to TOA radiance + + See + `Absolute Radiometric Calibration: 2016v0 `_ + and + `Improvements in Calibration, and Validation of the Absolute Radiometric Response of MAXAR Earth-Observing Sensors + `_ + (3.2.2) for more information. + + Args: + rad_arr (xr.DataArray): TOA Radiance array + band (BandNames): Band + + Returns: + xr.DataArray: TOA Reflectance array + """ + band_mtd_str = f"BAND_{MAXAR_BAND_MTD[band]}" + + # Get MTD XML file + root, _ = self.read_mtd() + + # Open zenith and azimuth angle + try: + band_mtd = root.find(f".//{band_mtd_str}") + abs_factor = float(band_mtd.findtext(".//ABSCALFACTOR")) + effective_bandwidth = float(band_mtd.findtext(".//EFFECTIVEBANDWIDTH")) + except TypeError: + raise InvalidProductError( + "ABSCALFACTOR or EFFECTIVEBANDWIDTH not found in metadata!" + ) + + gain, offset = MAXAR_GAIN_OFFSET[self.platform][band] + + coeff = gain * (abs_factor / effective_bandwidth) + offset + + LOGGER.debug(f"DN to rad coeff = {coeff}") + + toa_rad = coeff * dn_arr.data + + return dn_arr.copy(data=toa_rad) + + def _toa_rad_to_toa_refl( + self, rad_arr: xr.DataArray, band: BandNames + ) -> xr.DataArray: + """ + Compute TOA reflectance from TOA radiance + + See + `Absolute Radiometric Calibration: 2016v0 `_ + and + `Improvements in Calibration, and Validation of the Absolute Radiometric Response of MAXAR Earth-Observing Sensors + `_ + (3.2.2) for more information. + + Args: + rad_arr (xr.DataArray): TOA Radiance array + band (BandNames): Band + + Returns: + xr.DataArray: TOA Reflectance array + """ + _, sun_zenith_angle = self.get_mean_sun_angles() + toa_refl_coeff = ( + np.pi + * self._sun_earth_distance_variation() ** 2 + / (MAXAR_E0[self.platform][band] * np.cos(np.deg2rad(sun_zenith_angle))) + ) + + LOGGER.debug(f"rad to refl coeff = {toa_refl_coeff}") + + return rad_arr.copy(data=toa_refl_coeff * rad_arr) diff --git a/eoreader/products/optical/optical_product.py b/eoreader/products/optical/optical_product.py index ac491c40..aa48122a 100644 --- a/eoreader/products/optical/optical_product.py +++ b/eoreader/products/optical/optical_product.py @@ -17,6 +17,7 @@ """ Super class for optical products """ import logging from abc import abstractmethod +from datetime import datetime from enum import unique from pathlib import Path from typing import Union @@ -24,6 +25,7 @@ import geopandas as gpd import numpy as np import rasterio +import xarray as xr from cloudpathlib import CloudPath from rasterio import crs as riocrs from rasterio.enums import Resampling @@ -227,6 +229,27 @@ def get_existing_band_paths(self) -> dict: existing_bands = self.get_existing_bands() return self.get_band_paths(band_list=existing_bands) + def _to_reflectance( + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ + raise NotImplementedError + def _open_bands( self, band_paths: dict, @@ -257,7 +280,8 @@ def _open_bands( ) band_arr.attrs["radiometry"] = "as_is" - if TO_REFLECTANCE: + if kwargs.get(TO_REFLECTANCE, True): + # TODO: check if needs reflectance (NO TIR !!!) Landsats, S3 band_arr = self._to_reflectance(band_arr, band_path, band) band_arr.attrs["radiometry"] = "reflectance" @@ -698,3 +722,22 @@ def _get_cloud_band_path( return self._get_band_folder(writable).joinpath( f"{self.condensed_name}_{band.name}_{res_str.replace('.', '-')}.tif", ) + + @cache + def _sun_earth_distance_variation(self): + """ + Compute Sun-Earth distance variation. + It utilises the inverse square law of irradiance, under which, + the intensity (or irradiance) of light radiating from a point source is inversely proportional to the square of the distance from the source. + + See `here `_ for more information. + + Returns: + + """ + # julian_date is the Julian Day corresponding to the acquisition date (reference day: 01/01/1950). + ref_julian_date = datetime(year=1950, month=1, day=1) + julian_date = (self.date - ref_julian_date).days + 1 + d = 1 / (1 - 0.01673 * np.cos(0.0172 * (julian_date - 2) ** 2)) + LOGGER.debug(f"d = {d}") + return d diff --git a/eoreader/products/optical/pla_product.py b/eoreader/products/optical/pla_product.py index 3ae42252..a9b91d84 100644 --- a/eoreader/products/optical/pla_product.py +++ b/eoreader/products/optical/pla_product.py @@ -29,7 +29,7 @@ import geopandas as gpd import numpy as np -import xarray +import xarray as xr from cloudpathlib import CloudPath from lxml import etree from rasterio.enums import Resampling @@ -455,8 +455,24 @@ def _read_band( return band_arr def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ # Get MTD XML file root, nsmap = self.read_mtd() @@ -741,7 +757,7 @@ def open_mask( mask_id: str, resolution: float = None, size: Union[list, tuple] = None, - ) -> Union[xarray.DataArray, None]: + ) -> Union[xr.DataArray, None]: """ Open a Planet UDM2 (Usable Data Mask) mask, band by band, as a xarray. Returns None if the mask is not available. @@ -810,7 +826,7 @@ def _load_nodata( self, resolution: float = None, size: Union[list, tuple] = None, - ) -> Union[xarray.DataArray, None]: + ) -> Union[xr.DataArray, None]: """ Load nodata (unimaged pixels) as a numpy array. diff --git a/eoreader/products/optical/s2_product.py b/eoreader/products/optical/s2_product.py index 9d146f82..88be16df 100644 --- a/eoreader/products/optical/s2_product.py +++ b/eoreader/products/optical/s2_product.py @@ -553,8 +553,24 @@ def _read_band( ) def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ if str(path).endswith(".jp2"): try: # Get MTD XML file diff --git a/eoreader/products/optical/s2_theia_product.py b/eoreader/products/optical/s2_theia_product.py index a9baea70..cf4a9a5c 100644 --- a/eoreader/products/optical/s2_theia_product.py +++ b/eoreader/products/optical/s2_theia_product.py @@ -267,8 +267,24 @@ def get_band_paths( return band_paths def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ # Compute the correct radiometry of the band (Theia product are stored into int16 bits) original_dtype = band_arr.encoding.get("dtype", band_arr.dtype) if original_dtype == "int16": diff --git a/eoreader/products/optical/s3_olci_product.py b/eoreader/products/optical/s3_olci_product.py index f1ab4f0c..92bbd068 100644 --- a/eoreader/products/optical/s3_olci_product.py +++ b/eoreader/products/optical/s3_olci_product.py @@ -315,12 +315,6 @@ def _geocode( **{"SRC_METHOD": "GCP_TPS"}, ) - def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): - # Do nothing, managed elsewhere - return band_arr - def _rad_2_refl(self, band_arr: xr.DataArray, band: obn = None) -> xr.DataArray: """ Convert radiance to reflectance diff --git a/eoreader/products/optical/s3_product.py b/eoreader/products/optical/s3_product.py index 6dbd855a..54a8fe74 100644 --- a/eoreader/products/optical/s3_product.py +++ b/eoreader/products/optical/s3_product.py @@ -48,6 +48,7 @@ from eoreader.bands import BandNames from eoreader.bands import OpticalBandNames as obn from eoreader.exceptions import InvalidProductError +from eoreader.keywords import TO_REFLECTANCE from eoreader.products import OpticalProduct from eoreader.utils import DATETIME_FMT, EOREADER_NAME @@ -391,7 +392,10 @@ def get_band_paths( else: # Pre-process the wanted band (does nothing if existing) band_paths[band] = self._preprocess( - band, resolution=resolution, **kwargs + band, + resolution=resolution, + to_reflectance=kwargs.get(TO_REFLECTANCE, True), + **kwargs, ) return band_paths @@ -432,6 +436,28 @@ def _read_band( # Read band return band.astype(np.float32) * band.scale_factor + def _to_reflectance( + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ + # Do nothing, managed elsewhere + return band_arr + @abstractmethod def _manage_invalid_pixels( self, band_arr: XDS_TYPE, band: obn, **kwargs diff --git a/eoreader/products/optical/s3_slstr_product.py b/eoreader/products/optical/s3_slstr_product.py index d02be68a..5906ae87 100644 --- a/eoreader/products/optical/s3_slstr_product.py +++ b/eoreader/products/optical/s3_slstr_product.py @@ -350,12 +350,6 @@ def _set_product_type(self) -> None: } ) - def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): - # Do nothing, managed elsewhere - return band_arr - def _preprocess( self, band: Union[obn, str], diff --git a/eoreader/products/optical/vis1_product.py b/eoreader/products/optical/vis1_product.py index 8e1c7f02..89f485db 100644 --- a/eoreader/products/optical/vis1_product.py +++ b/eoreader/products/optical/vis1_product.py @@ -28,6 +28,7 @@ from typing import Union import numpy as np +import xarray as xr from cloudpathlib import CloudPath from lxml import etree from rasterio import crs as riocrs @@ -43,6 +44,20 @@ LOGGER = logging.getLogger(EOREADER_NAME) +VIS1_E0 = { + obn.PAN: 1828, + obn.BLUE: 2003, + obn.GREEN: 1828, + obn.RED: 1618, + obn.NIR: 1042, + obn.NARROW_NIR: 1042, +} +""" +Solar spectral irradiance, E0b, (commonly known as ESUN) is a constant value specific to each band of the Vision-1 imager. +It is determined by using well know models of Solar Irradiance with the measured spectral transmission of the imager for each incident wavelength. +It has units of Wm-2μm-1. The applicable values for Vision-1 are provided in the table. +""" + @unique class Vis1BandCombination(ListEnum): @@ -315,21 +330,35 @@ def get_mean_sun_angles(self) -> (float, float): return azimuth_angle, zenith_angle def _to_reflectance( - self, band_arr, path: Union[Path, CloudPath], band: BandNames, **kwargs - ): + self, + band_arr: xr.DataArray, + path: Union[Path, CloudPath], + band: BandNames, + **kwargs, + ) -> xr.DataArray: + """ + Converts band to reflectance + + Args: + band_arr (xr.DataArray): + path (Union[Path, CloudPath]): + band (BandNames): + **kwargs: Other keywords + + Returns: + xr.DataArray: Band in reflectance + """ # Compute the correct radiometry of the band original_dtype = band_arr.encoding.get("dtype", band_arr.dtype) if original_dtype == "uint16": band_arr /= 100.0 - # TODO: Convert into reflectance ! - # To float32 if band_arr.dtype != np.float32: band_arr = band_arr.astype(np.float32) - return band_arr + return self._toa_rad_to_toa_refl(band_arr, band) @cache def _read_mtd(self) -> (etree._Element, dict): @@ -394,3 +423,31 @@ def _get_ortho_path(self, **kwargs) -> Union[CloudPath, Path]: rpcs = utils.open_rpc_file(rpcs_file) return super()._get_ortho_path(rpcs=rpcs, **kwargs) + + def _toa_rad_to_toa_refl( + self, rad_arr: xr.DataArray, band: BandNames + ) -> xr.DataArray: + """ + Compute TOA reflectance from TOA radiance + + See + `here `_ + (3.2.2) for more information. + + Args: + rad_arr (xr.DataArray): TOA Radiance array + band (BandNames): Band + + Returns: + xr.DataArray: TOA Reflectance array + """ + _, sun_zenith_angle = self.get_mean_sun_angles() + toa_refl_coeff = ( + np.pi + * self._sun_earth_distance_variation() ** 2 + / (VIS1_E0[band] * np.cos(np.deg2rad(sun_zenith_angle))) + ) + + LOGGER.debug(f"rad to refl coeff = {toa_refl_coeff}") + + return rad_arr.copy(data=toa_refl_coeff * rad_arr)