- ENH: Adding the support of Landsat Level-2 products #49

- FIX: Computing Brightness Temperature of `Landsat` TIR bands instead of leaving them as is
sertit · Sep 2, 2022 · fe9d554 · fe9d554
1 parent 0f30b3e
commit fe9d554
Show file tree

Hide file tree

Showing 3 changed files with 115 additions and 54 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -5,6 +5,7 @@
 ### Enhancements
 
 - **ENH: Adding the support of RapidEye constellation**
+- **ENH: Adding the support of Landsat Level-2 products** ([#49](https://github.com/sertit/eoreader/issues/49))
 
 ### Bug Fixes
 
@@ -15,6 +16,7 @@
 - FIX: Fixing computation of invalid pixels for `Sentinel-2` and `DIMAP` products (do not remove straylight mask)
 - FIX: Fixing pandas FutureWarning `The frame.append method is deprecated and will be removed from pandas in a future version.`
 - FIX: Fixing DeprecationWarning `invalid escape sequence \.`
+- FIX: Computing Brightness Temperature of `Landsat` TIR bands instead of leaving them as is
 
 ### Optimizations
 

diff --git a/eoreader/products/optical/landsat_product.py b/eoreader/products/optical/landsat_product.py
@@ -66,26 +66,27 @@ class LandsatProductType(ListEnum):
     Each Level-1 data product includes individual spectral band files, a metadata file, and additional ancillary files.
     """
 
-    ARD = "ARD"
+    L2 = "L2"
     """
-    Uses Landsat Collections Level-1 data as input
-    to provide data that is processed to the highest scientific standards and placed in a tile-based structure to support time-series analysis.
+    Level-2 Science Products are time-series observational data of sufficient length, consistency, and continuity to record effects of climate change,
+    and serve as input into Landsat Level-3 Science Products.
 
-    Not handled by EOReader.
+    Only for Landsat 4 to 9.
+
+    Only Surface Reflectances and Temperatures are handled by EOReader
     """
 
-    L2 = "L2"
+    L3 = "L3"
     """
-    Level-2 and Level-3 products that are processed to include
-    atmospherically corrected data, surface reflectance, provisional surface temperature, and biophysical properties of the Earth’s surface.
+    Level-3 science products represent biophysical properties of the Earth's surface and are generated from Landsat U.S. Analysis Ready Data (ARD) inputs.
 
     Not handled by EOReader.
     """
 
-    L3 = "L3"
+    ARD = "ARD"
     """
-    Level-2 and Level-3 products that are processed to include
-    atmospherically corrected data, surface reflectance, provisional surface temperature, and biophysical properties of the Earth’s surface.
+    Uses Landsat Collections Level-1 data as input
+    to provide data that is processed to the highest scientific standards and placed in a tile-based structure to support time-series analysis.
 
     Not handled by EOReader.
     """
@@ -334,9 +335,9 @@ def _set_product_type(self) -> None:
 
             self.product_type = LandsatProductType.from_value(proc_lvl[:-2])
 
-            if self.product_type != LandsatProductType.L1:
+            if self.product_type not in [LandsatProductType.L1, LandsatProductType.L2]:
                 LOGGER.warning(
-                    "Only Landsat level 1 have been tested on EOReader, ise it at your own risk."
+                    "Only Landsat level 1 and 2 have been tested on EOReader, ise it at your own risk."
                 )
             else:
                 # Warning if GS (L1 only)
@@ -1099,53 +1100,111 @@ def _to_reflectance(
         if not (self._pixel_quality_id in filename or self._radsat_id in filename):
             # Convert raw bands from DN to correct reflectance
             if not filename.startswith(self.condensed_name):
-                # Original band name
-                band_name = filename[-1]
-
                 # Open mtd
                 mtd_data, _ = self._read_mtd()
 
                 # Get band nb and corresponding coeff
-                c_mul_str = "REFLECTANCE_MULT_BAND_" + band_name
-                c_add_str = "REFLECTANCE_ADD_BAND_" + band_name
-
-                # Get coeffs to convert DN to reflectance
+                band_name = filename[-1]
                 try:
-                    c_mul = mtd_data.findtext(f".//{c_mul_str}")
-                    c_add = mtd_data.findtext(f".//{c_add_str}")
+                    # Thermal (10/11)
+                    if band in [spb.TIR_1, spb.TIR_2]:
+                        band_name = filename[-2:]
+                        band_arr = self._to_tb(band_arr, band, mtd_data, band_name)
 
-                    # Manage some cases where the values are set to NULL
-                    if c_mul == "NULL":
-                        c_mul = 1.0
-                    else:
-                        c_mul = float(c_mul)
-                    if c_add == "NULL":
-                        c_add = 1.0
                     else:
-                        c_add = float(c_add)
+                        # Original band name
+                        band_arr = self._to_refl(band_arr, band, mtd_data, band_name)
                 except TypeError:
-                    if band in [spb.TIR_1, spb.TIR_2]:
-                        c_mul = 1.0
-                        c_add = 0.0
-                    else:
-                        raise InvalidProductError(
-                            f"Cannot find additive or multiplicative "
-                            f"rescaling factor for bands ({band.name}, "
-                            f"number {band_name}) in metadata"
-                        )
+                    raise InvalidProductError(
+                        f"Cannot find additive or multiplicative "
+                        f"rescaling factor for bands ({band.name}, "
+                        f"number {band_name}) in metadata"
+                    )
 
-                # Manage NULL values
-                try:
-                    c_mul = float(c_mul)
-                except ValueError:
-                    c_mul = 1.0
-                try:
-                    c_add = float(c_add)
-                except ValueError:
-                    c_add = 0.0
+        return band_arr
+
+    def _to_tb(
+        self,
+        band_arr: xr.DataArray,
+        band: BandNames,
+        mtd_data,
+        band_name,
+        **kwargs,
+    ) -> xr.DataArray:
+        """
+        Converts band to brightness temperatue
+
+        https://www.usna.edu/Users/oceano/pguth/md_help/remote_sensing_course/landsat_thermal.htm
+
+
+        Args:
+
+        Returns:
+            xr.DataArray: Band in brightness temperature
+        """
+        # Get coeffs to convert DN to radiance
+        c_mul_str = "RADIANCE_MULT_BAND_" + band_name
+        c_add_str = "RADIANCE_ADD_BAND_" + band_name
+        c_mul = mtd_data.findtext(f".//{c_mul_str}")
+        c_add = mtd_data.findtext(f".//{c_add_str}")
+
+        # Manage NULL values
+        try:
+            c_mul = float(c_mul)
+        except ValueError:
+            c_mul = 1.0
+        try:
+            c_add = float(c_add)
+        except ValueError:
+            c_add = 0.0
+
+        band_arr = c_mul * band_arr + c_add
+
+        # Get coeffs to convert radiance to TB
+        k1_str = "K1_CONSTANT_BAND_" + band_name
+        k2_str = "K2_CONSTANT_BAND_" + band_name
+        k1 = float(mtd_data.findtext(f".//{k1_str}"))
+        k2 = float(mtd_data.findtext(f".//{k2_str}"))
+
+        band_arr = k2 / np.log(k1 / band_arr + 1)
+
+        return band_arr
+
+    def _to_refl(
+        self,
+        band_arr: xr.DataArray,
+        band: BandNames,
+        mtd_data,
+        band_name,
+        **kwargs,
+    ) -> xr.DataArray:
+        """
+        Converts band to reflectance
+
+        Args:
+
+        Returns:
+            xr.DataArray: Band in reflectance
+        """
+        c_mul_str = "REFLECTANCE_MULT_BAND_" + band_name
+        c_add_str = "REFLECTANCE_ADD_BAND_" + band_name
+
+        # Get coeffs to convert DN to reflectance
+        c_mul = mtd_data.findtext(f".//{c_mul_str}")
+        c_add = mtd_data.findtext(f".//{c_add_str}")
+
+        # Manage NULL values
+        try:
+            c_mul = float(c_mul)
+        except ValueError:
+            c_mul = 1.0
+        try:
+            c_add = float(c_add)
+        except ValueError:
+            c_add = 0.0
 
-                # Compute the correct reflectance of the band and set no data to 0
-                band_arr = c_mul * band_arr + c_add  # Already in float
+        # Compute the correct reflectance of the band and set no data to 0
+        band_arr = c_mul * band_arr + c_add  # Already in float
 
         return band_arr
 

diff --git a/eoreader/reader.py b/eoreader/reader.py
@@ -192,11 +192,11 @@ class Constellation(ListEnum):
     Constellation.S2_THEIA: r"SENTINEL2[AB]_\d{8}-\d{6}-\d{3}_L(2A|1C)_T\d{2}\w{3}_[CDH](_V\d-\d|)",
     Constellation.S3_OLCI: r"S3[AB]_OL_[012]_\w{6}_\d{8}T\d{6}_\d{8}T\d{6}_\d{8}T\d{6}_\w{17}_\w{3}_[OFDR]_(NR|ST|NT)_\d{3}",
     Constellation.S3_SLSTR: r"S3[AB]_SL_[012]_\w{6}_\d{8}T\d{6}_\d{8}T\d{6}_\d{8}T\d{6}_\w{17}_\w{3}_[OFDR]_(NR|ST|NT)_\d{3}",
-    Constellation.L9: r"L[OTC]09_L1(GT|TP)_\d{6}_\d{8}_\d{8}_\d{2}_(RT|T1|T2)",
-    Constellation.L8: r"L[OTC]08_L1(GT|TP)_\d{6}_\d{8}_\d{8}_\d{2}_(RT|T1|T2)",
-    Constellation.L7: r"LE07_L1(GT|TP|GS)_\d{6}_\d{8}_\d{8}_\d{2}_(RT|T1|T2)",
-    Constellation.L5: r"L[TM]05_L1(TP|GS)_\d{6}_\d{8}_\d{8}_\d{2}_(T1|T2)",
-    Constellation.L4: r"L[TM]04_L1(TP|GS)_\d{6}_\d{8}_\d{8}_\d{2}_(T1|T2)",
+    Constellation.L9: r"L[OTC]09_(L1(GT|TP)|L2(SP|SR))_\d{6}_\d{8}_\d{8}_\d{2}_(RT|T1|T2)",
+    Constellation.L8: r"L[OTC]08_(L1(GT|TP)|L2(SP|SR))_\d{6}_\d{8}_\d{8}_\d{2}_(RT|T1|T2)",
+    Constellation.L7: r"LE07_(L1(GT|TP|GS)|L2(SP|SR))_\d{6}_\d{8}_\d{8}_\d{2}_(RT|T1|T2)",
+    Constellation.L5: r"L[TM]05_(L1(TP|GS)|L2(SP|SR))_\d{6}_\d{8}_\d{8}_\d{2}_(T1|T2)",
+    Constellation.L4: r"L[TM]04_(L1(TP|GS)|L2(SP|SR))_\d{6}_\d{8}_\d{8}_\d{2}_(T1|T2)",
     Constellation.L3: r"LM03_L1(TP|GS)_\d{6}_\d{8}_\d{8}_\d{2}_T2",
     Constellation.L2: r"LM02_L1(TP|GS)_\d{6}_\d{8}_\d{8}_\d{2}_T2",
     Constellation.L1: r"LM01_L1(TP|GS)_\d{6}_\d{8}_\d{8}_\d{2}_T2",