From fe9d5541309d9002a838c35de67a4a801a90c1d7 Mon Sep 17 00:00:00 2001 From: BRAUN REMI Date: Fri, 2 Sep 2022 14:30:21 +0200 Subject: [PATCH] - ENH: Adding the support of Landsat Level-2 products #49 - FIX: Computing Brightness Temperature of `Landsat` TIR bands instead of leaving them as is --- CHANGES.md | 2 + eoreader/products/optical/landsat_product.py | 157 +++++++++++++------ eoreader/reader.py | 10 +- 3 files changed, 115 insertions(+), 54 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index d3b65c5b..2d3b14c9 100644 --- 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 index 615ab8c9..c0fc6b3b 100644 --- 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 index 88068d6e..fb1f9d21 100644 --- 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",