Theia glint corrected

Author: cordmaur

waterdetect/Common.py

@@ -2,9 +2,8 @@
 from shutil import copy
 import configparser
 import ast
-
 from sklearn.model_selection import train_test_split
-from sklearn.preprocessing import MinMaxScaler
+from sklearn.preprocessing import MinMaxScaler, RobustScaler
 from sklearn.preprocessing import QuantileTransformer
 from waterdetect import DWProducts, gdal
 
@@ -449,6 +448,24 @@ def calc_normalized_difference(img1, img2, mask=None, compress_cte=0.02):
 
         return nd.filled(), nd.mask
 
+    @staticmethod
+    def calc_mbwi(bands, factor, mask):
+        # changement for negative SRE values scene
+        min_cte = np.min([np.min(bands['Green'][~mask]), np.min(bands['Red'][~mask]),
+                          np.min(bands['Nir'][~mask]), np.min(bands['Mir'][~mask]), np.min(bands['Mir2'][~mask])])
+        if min_cte <= 0:
+            min_cte = -min_cte + 0.001
+        else:
+            min_cte = 0
+        mbwi = factor * (bands['Green'] + min_cte) - (bands['Red'] + min_cte) - (bands['Nir'] + min_cte) \
+               - (bands['Mir'] + min_cte) - (bands['Mir2'] + min_cte)
+        mbwi[~mask] = RobustScaler(copy=False).fit_transform(mbwi[~mask].reshape(-1, 1)).reshape(-1)
+        mbwi[~mask] = MinMaxScaler(feature_range=(-1, 1), copy=False).fit_transform(mbwi[~mask].reshape(-1, 1)) \
+            .reshape(-1)
+        mask = np.isinf(mbwi) | np.isnan(mbwi) | mask
+        mbwi = np.ma.array(mbwi, mask=mask, fill_value=-9999)
+        return mbwi, mask
+
     @staticmethod
     def rgb_burn_in(red, green, blue, burn_in_array, color=None, min_value=None, max_value=None, colormap='viridis',
                     fade=1, uniform_distribution=False, no_data_value=-9999, valid_value=1, transp=0.0):

waterdetect/Glint.py

@@ -5,19 +5,32 @@
 
 
 class DWGlintProcessor:
+    supported_products = ['S2_S2COR', 'S2_THEIA']
+
     def __init__(self, image, limit_angle=30):
+
         self.image = image
         self.limit_angle = limit_angle
         self.adjustable_bands = ['Mir', 'Mir2', 'Nir', 'Nir2']
 
         try:
-            self.glint_array = self.create_glint_array(self.image.metadata)
+            self.glint_array = self.create_glint_array(self.image.metadata, image.product)
 
         except BaseException as err:
             self.glint_array = None
             print(f'### GLINT PROCESSOR ERROR #####')
             print(err)
 
+    @classmethod
+    def create(cls, image, limit_angle=30):
+        if image.product not in cls.supported_products:
+            print(f'Product {image.product} not supported by GlintProcessor')
+            print(f'Supported products: {cls.supported_products}')
+            return None
+        else:
+            return cls(image, limit_angle)
+
+
     @staticmethod
     def get_grid_values_from_xml(tree_node, xpath_str):
         """Receives a XML tree node and a XPath parsing string and search for children matching the string.
@@ -35,14 +48,16 @@ def get_grid_values_from_xml(tree_node, xpath_str):
         return np.nanmean(arrays_lst, axis=0)
 
     @staticmethod
-    def create_glint_array(xml_file):
+    def create_glint_array(xml_file, product):
         xml_file = Path(xml_file)
         parser = etree.XMLParser()
         root = etree.parse(xml_file.as_posix(), parser).getroot()
 
-        sun_zenith = np.deg2rad(DWGlintProcessor.get_grid_values_from_xml(root, './/Sun_Angles_Grid/Zenith'))[:-1, :-1]
-        sun_azimuth = np.deg2rad(DWGlintProcessor.get_grid_values_from_xml(root, './/Sun_Angles_Grid/Azimuth'))[:-1,
-                      :-1]
+        sun_angles = 'Sun_Angles_Grid' if product == 'S2_S2COR' else 'Sun_Angles_Grids'
+        # viewing_angles = 'Viewing_Incidence_Angles_Grids'
+
+        sun_zenith = np.deg2rad(DWGlintProcessor.get_grid_values_from_xml(root, f'.//{sun_angles}/Zenith'))[:-1, :-1]
+        sun_azimuth = np.deg2rad(DWGlintProcessor.get_grid_values_from_xml(root, f'.//{sun_angles}/Azimuth'))[:-1,:-1]
 
         view_zenith = np.deg2rad(
             DWGlintProcessor.get_grid_values_from_xml(root, './/Viewing_Incidence_Angles_Grids/Zenith'))[:-1, :-1]

waterdetect/Image.py

@@ -66,6 +66,12 @@ def check_necessary_bands(self, bands, bands_keys, invalid_mask):
             invalid_mask |= ndwi_mask
             bands.update({'ndwi': ndwi})
 
+        # check if the MBWI index exist
+        if 'mbwi' not in bands.keys():
+            mbwi, mbwi_mask = DWutils.calc_mbwi(bands, 3, invalid_mask)
+            invalid_mask |= ndwi_mask
+            bands.update({'mbwi': mbwi})
+
         # todo: check the band for Principal Component Analysis
 
         # check if the list contains the required bands
@@ -649,7 +655,7 @@ def apply_clustering(self):
 
         for band, value in zip(self.config.clip_band, self.config.clip_sup_value):
             if value is not None:
-                if self.config.glint_mode:
+                if self.config.glint_mode and (self.glint_processor is not None):
                     comp_array = self.glint_processor.glint_adjusted_threshold(band,
                                                                                value,
                                                                                'SUP',
@@ -663,7 +669,7 @@ def apply_clustering(self):
         # after obtaining the final labels, clip bands with inferior limit
         for band, value in zip(self.config.clip_band, self.config.clip_inf_value):
             if value is not None:
-                if self.config.glint_mode:
+                if self.config.glint_mode and (self.glint_processor is not None):
                     comp_array = self.glint_processor.glint_adjusted_threshold(band,
                                                                                value,
                                                                                'INF',

waterdetect/WaterDetect.py

@@ -137,24 +137,7 @@ def calc_mbwi(self, bands, factor=3, save_index=False):
 
         mask = self.loader.invalid_mask
 
-        # changement for negative SRE values scene
-        min_cte = np.min([np.min(bands['Green'][~mask]), np.min(bands['Red'][~mask]),
-                          np.min(bands['Nir'][~mask]), np.min(bands['Mir'][~mask]), np.min(bands['Mir2'][~mask])])
-
-        if min_cte <= 0:
-            min_cte = -min_cte + 0.001
-        else:
-            min_cte = 0
-
-        mbwi = factor * (bands['Green']+min_cte) - (bands['Red']+min_cte) - (bands['Nir']+min_cte)\
-            - (bands['Mir']+min_cte) - (bands['Mir2']+min_cte)
-
-        mbwi[~mask] = RobustScaler(copy=False).fit_transform(mbwi[~mask].reshape(-1, 1)).reshape(-1)
-        mbwi[~mask] = MinMaxScaler(feature_range=(-1, 1), copy=False).fit_transform(mbwi[~mask].reshape(-1, 1))\
-            .reshape(-1)
-
-        mask = np.isinf(mbwi) | np.isnan(mbwi) | self.loader.invalid_mask
-        mbwi = np.ma.array(mbwi, mask=mask, fill_value=-9999)
+        mbwi, mask = DWutils.calc_mbwi(bands, factor, mask)
 
         self.loader.update_mask(mask)
 
@@ -165,6 +148,7 @@ def calc_mbwi(self, bands, factor=3, save_index=False):
 
         return mbwi.filled()
 
+
     def calc_awei(self, bands, save_index=False):
         """
         Calculates the AWEI Water Index and adds it to the bands dictionary
@@ -373,12 +357,16 @@ def create_mask_report(self, image, band_combination, composite_name, pdf_merger
         # calculate the sun glint rejection and add it to the pdf report
         # the glint will be passed to the
         if self.config.calc_glint:
-            glint_processor = DWGlintProcessor(image)
-            pdf_merger_image.append(glint_processor.save_heatmap(self.saver.output_folder))
+            glint_processor = DWGlintProcessor.create(image)
+
+            # if there is a valid glint_processor, save the heatmap
+            if glint_processor is not None:
+                pdf_merger_image.append(glint_processor.save_heatmap(self.saver.output_folder))
+            else:
+                print(f'Glint_mode is On but no Glint Processor is available for this product')
 
         else:
             glint_processor = None
-            # self.calc_glint(image, self.saver.output_folder, pdf_merger_image)
 
         # create a dw_image object with the water mask and all the results
         dw_image = self.create_water_mask(band_combination, pdf_merger_image, glint_processor=glint_processor)

waterdetect/__init__.py

@@ -1,6 +1,6 @@
 # todo: Implement logging
 # import logging
-__version__ = '1.5.7'
+__version__ = '1.5.8'
 
 class DWProducts:
     Landsat8_USGS = 'L8_USGS'