FIX: Loading every optical band in reflectance (fixed for Maxar, `P…

…lanet` and `Vision-1` data) (sertit#30)

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

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(

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:

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 <https://earth.esa.int/eogateway/documents/20142/37627/DigitalGlobe-Standard-Imagery.pdf>`_
 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 <https://apollomapping.com/image_downloads/Maxar_AbsRadCalDataSheet2018v0.pdf>_` 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 <https://apollomapping.com/image_downloads/Maxar_AbsRadCalDataSheet2018v0.pdf>_` 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 <https://dg-cms-uploads-production.s3.amazonaws.com/uploads/document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf>`_
+        and
+        `Improvements in Calibration, and Validation of the Absolute Radiometric Response of MAXAR Earth-Observing Sensors
+        <https://dg-cms-uploads-production.s3.amazonaws.com/uploads/document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf>`_
+        (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 <https://dg-cms-uploads-production.s3.amazonaws.com/uploads/document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf>`_
+        and
+        `Improvements in Calibration, and Validation of the Absolute Radiometric Response of MAXAR Earth-Observing Sensors
+        <https://dg-cms-uploads-production.s3.amazonaws.com/uploads/document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf>`_
+        (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)