WQ algorithms based on OWTs

getpak/commons.py

@@ -89,14 +89,31 @@ def find_nearest(array, value):
         return idx
 
     @staticmethod
-    def walktalk(path2walk, fpattern='.nc', dir_is_file=False, badstring=None):
+    def walktalk(path2walk, fpattern='.nc', dir_is_file=False, unwanted_string=None):
+        '''
+
+        Parameters
+        ----------
+        path2walk: folder to scan for files
+        fpattern: the file pattern to look for inside the folders
+        dir_is_file: option to include the folders in the search, defaults to False
+        unwanted_string: a string to search for and remove from the list of files if present
+
+        Returns
+        -------
+        output_flist: list of files (and folders) matching the fpattern found inside the path2walk folder
+        '''
         output_flist = []
         for root, dirs, files in os.walk(path2walk, topdown=False):
             for name in (dirs if dir_is_file else files):
                 # check the file format
                 if name.endswith(fpattern):
+                    if isinstance(unwanted_string, str):
                     # check if it is not the file you want
-                    if not fnmatch.fnmatch(name, badstring):
+                        if not fnmatch.fnmatch(name, unwanted_string):
+                            f = Path(os.path.join(root, name))
+                            output_flist.append(f)
+                    else:
                         f = Path(os.path.join(root, name))
                         output_flist.append(f)
         print(f'{len(output_flist)} files found.')

getpak/inversion_functions.py

@@ -13,7 +13,52 @@
 import numpy as np
 
 
-# Below is an example extracted from Nechad et al. (2010)
+#### Chlorophyll-a
+# Gitelson
+def chl_gitelson(Red, RedEdg1, RedEdg2):
+    chl = 23.1 + 117.4 * (1 / Red - 1 / RedEdg1) * RedEdg2
+    return chl
+
+# Gitelson and Kondratyev, Dall'Olmo et al. (2003)
+def chl_gitelson2(Red, RedEdg1, a=61.324, b=-37.94):
+    chl = a * (RedEdg1 / Red) + b
+    return chl
+
+# OC2
+def chl_OC2(Blue, Green, a=0.2389, b=-1.9369, c=1.7627, d=-3.0777, e=-0.1054):
+    X = np.log10(Blue / Green)
+    chl = 10**(a + b*X + c*X**2 + d*X**3 + e*X**4)
+    return chl
+
+# 2-band ratio by Gilerson et al. (2010)
+def chl_gilerson2(Red, RedEdg1, a=35.75, b=1.124):
+    chl = (a * (RedEdg1 / Red) - 19.30)**b
+    return chl
+
+# 2-band semi-analytical by Gons et al. (2003, 2005)
+def chl_gons(Red, RedEdg1, RedEdg3, a=1.063, b=0.016, aw664=0.40, aw708=0.70):
+    bb = 1.61 * RedEdg3 / (0.082 - 0.6 * RedEdg3)
+    chl = ((RedEdg1 / Red) * (aw708 + bb) - aw664 - bb**a)/b
+    return chl
+
+# JM Hybride 1
+def chl_h1(Red, RedEdge1, RedEdge2):
+    res = RedEdge1 - Red  # B5 - B4
+    chl = np.zeros_like(res)
+    chl[res < 0] = 115.107 * RedEdge2[res < 0] * (1/Red[res < 0] - 1/RedEdge1[res < 0]) + 16.56
+    chl[res >= 0] = 115.794 * RedEdge2[res >= 0] * (1/Red[res >= 0] - 1/RedEdge1[res >= 0]) + 20.678
+    return chl
+
+# JM Hybride 2
+def chl_h2(Red, RedEdge1, RedEdge2):
+    res = RedEdge1 - Red  # B5 - B4
+    chl = np.zeros_like(res)
+    chl[res < 0] = 46.859 * RedEdge1[res < 0]/Red[res < 0] - 29.916
+    chl[res >= 0] = 115.794 * RedEdge2[res >= 0] * (1/Red[res >= 0] - 1/RedEdge1[res >= 0]) + 20.678
+    return chl
+
+#### SPM
+# Nechad et al. (2010)
 def nechad(Red, a=610.94, c=0.2324):
     spm = a * Red / (1 - (Red / c))
     return spm
@@ -47,35 +92,8 @@ def spm_s3(Red, Nir2, cutoff_value=0.027, cutoff_delta=0.007, low_params=None, h
     spm = (1 - transition_coef) * low + transition_coef * high
     return spm
 
-
-# Gitelson
-def chl_gitelson(Red, RedEdg1, RedEdg2):
-    chl = 23.1 + 117.4 * (1 / Red - 1 / RedEdg1) * RedEdg2
-    return chl
-
-
-# Gitelson and Kondratyev
-def chl_gitelson2(Red, RedEdg1):
-    chl = 61.324 * (RedEdg1 / Red) - 37.94
-    return chl
-
-# JM Hybride 1
-def chl_h1(Red, RedEdge1, RedEdge2):
-    res = RedEdge1 - Red  # B5 - B4
-    chl = np.zeros_like(res)
-    chl[res < 0] = 115.107 * RedEdge2[res < 0] * (1/Red[res < 0] - 1/RedEdge1[res < 0]) + 16.56
-    chl[res >= 0] = 115.794 * RedEdge2[res >= 0] * (1/Red[res >= 0] - 1/RedEdge1[res >= 0]) + 20.678
-    return chl
-
-# JM Hybride 2
-def chl_h2(Red, RedEdge1, RedEdge2):
-    res = RedEdge1 - Red  # B5 - B4
-    chl = np.zeros_like(res)
-    chl[res < 0] = 46.859 * RedEdge1[res < 0]/Red[res < 0] - 29.916
-    chl[res >= 0] = 115.794 * RedEdge2[res >= 0] * (1/Red[res >= 0] - 1/RedEdge1[res >= 0]) + 20.678
-    return chl
-
-# Turbidity (FNU) Dogliotti
+#### Turbidity (FNU)
+# Dogliotti
 def turb_dogliotti(Red, Nir2):
     """Switching semi-analytical-algorithm computes turbidity from red and NIR band
 
@@ -100,7 +118,8 @@ def turb_dogliotti(Red, Nir2):
     t_low[t_low > 4000] = 0
     return t_low
 
-# CDOM Brezonik et al. 2005
+#### CDOM
+# Brezonik et al. 2005
 def cdom_brezonik(Blue, RedEdg2):
 
 
@@ -129,7 +148,22 @@ def cdom_brezonik(Blue, RedEdg2):
         'function': chl_gitelson2,
         'units': 'mg/m³'
     },
-
+
+    'CHL_OC2': {
+        'function': chl_OC2,
+        'units': 'mg/m³'
+    },
+
+    'CHL_Gilerson': {
+        'function': chl_gilerson2,
+        'units': 'mg/m³'
+    },
+
+    'CHL_Gons': {
+    'function': chl_gons,
+    'units': 'mg/m³'
+    },
+
     'CHL_Hybrid1': {
         'function': chl_h1,
         'units': 'mg/m³'

getpak/raster.py

@@ -323,8 +323,8 @@ def classify_owt(self, rasterio_rast, shapefiles, rrs_dict, bands=['Rrs_B2', 'Rr
         class_spt: an array, with the same size as the input bands, with the classified pixels
         class_shp: an array with the same length as the shapefiles, with a OWT class for each polygon
         """
-        class_spt = np.zeros(rrs_dict[bands[0]].shape, dtype='uint8')
-        class_shp = np.zeros((len(shapefiles)), dtype='uint8')
+        class_spt = np.zeros(rrs_dict[bands[0]].shape, dtype='int32')
+        class_shp = np.zeros((len(shapefiles)), dtype='int32')
         for i, shape in enumerate(shapefiles):
             values, slices, mask = self.extract_px(rasterio_rast, shape, rrs_dict, bands)
             # Verifying if there are more pixels than the minimum
@@ -341,7 +341,202 @@ def classify_owt(self, rasterio_rast, shapefiles, rrs_dict, bands=['Rrs_B2', 'Rr
             # classification by polygon
             class_shp[i] = angle
 
-        return class_spt, class_shp
+        return class_spt.astype('uint8'), class_shp.astype('uint8')
+
+    def water_colour(self, rasterio_rast, shapefiles, rrs_dict, bands=['Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5'],  min_px=6):
+        """
+        Function to calculate the water colour based on the Forel-Ule scale, using the bands of Sentinel-2 MSI, with the
+        coefficients derived by linear correlation by van der Woerd and Wernand (2018). The different combinations of S2
+        bands (10, 20 or 60 m resolution) in the visible spectrum can be used, with the default being at 20 m.
+
+        Parameters
+        ----------
+        rasterio_rast: a rasterio raster open with rasterio.open
+        shapefiles: a polygon (or set of polygons), usually of waterbodies to be classified, opened as geometry using fiona
+        rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands
+        bands: an array containing the bands of the rrs_dict to be extracted, in order from blue to red edge
+        min_px: minimum number of pixels in each polygon to operate the classification
+
+        Returns
+        -------
+        colour: an array, with the same size as the input bands, with the classified pixels
+        colour_spt: an array with the same length as the shapefiles, with the mean colour (angle) for each polygon
+        """
+        colour_spt = np.zeros(rrs_dict[bands[0]].shape, dtype='float16')
+        colour_shp = np.zeros((len(shapefiles)), dtype='float16')
+
+        for i, shape in enumerate(shapefiles):
+            values, slices, mask = self.extract_px(rasterio_rast, shape, rrs_dict, bands)
+            # Verifying if there are more pixels than the minimum
+            valid_pixels = np.isnan(values[0]) == False
+            if np.count_nonzero(valid_pixels) >= min_px:
+                med = np.nanmean(values, axis=1)
+                # calculation of the CIE tristimulus (no intercepts):
+                if len(bands) == 3:
+                    X = 12.040 * med[0] + 53.696 * med[1] + 32.087 * med[2]
+                    Y = 23.122 * med[0] + 65.702 * med[1] + 16.830 * med[2]
+                    Z = 61.055 * med[0] + 1.778 * med[1] + 0.015 * med[2]
+                    delta = lambda x: -164.83 * (alpha / 100) ** 5 + 1139.90 * (alpha / 100) ** 4 - 3006.04 * (
+                            alpha / 100) ** 3 + 3677.75 * (alpha / 100) ** 2 - 1979.71 * (alpha / 100) + 371.38
+                else if len(bands) == 4:
+                    X = 12.040 * med[0] + 53.696 * med[1] + 32.028 * med[2] + 0.529 * med[3]
+                    Y = 23.122 * med[0] + 65.702 * med[1] + 16.808 * med[2] + 0.192 * med[3]
+                    Z = 61.055 * med[0] + 1.778 * med[1] + 0.015 * med[2] + 0.000 * med[3]
+                    delta = lambda x : -161.23 * (alpha / 100) ** 5 + 1117.08 * (alpha / 100) ** 4 - 2950.14 * (
+                            alpha / 100) ** 3 + 3612.17 * (alpha / 100) ** 2 - 1943.57 * (alpha / 100) + 364.28
+                else if len(bands) == 5:
+                    X = 11.756 * med[0] + 6.423 * med[1] + 53.696 * med[2] + 32.028 * med[3] + 0.529 * med[4]
+                    Y = 1.744 * med[0] + 22.289 * med[1] + 65.702 * med[2] + 16.808 * med[3] + 0.192 * med[4]
+                    Z = 62.696 * med[0] + 31.101 * med[1] + 1.778 * med[2] + 0.015 * med[3] + 0.000 * med[4]
+                    delta = lambda x: -65.74 * (alpha / 100) ** 5 + 477.16 * (alpha / 100) ** 4 - 1279.99 * (
+                            alpha / 100) ** 3 + 1524.96 * (alpha / 100) ** 2 - 751.59 * (alpha / 100) + 116.56
+                else:
+                    print("Error in the number of bands provided")
+                    colour_spt, colour_shp = None, None
+                # normalisation of the tristimulus in 2 coordinates
+                x = X / (X + Y + Z)
+                y = Y / (X + Y + Z)
+                # hue angle:
+                alpha = (np.arctan2(y - 1 / 3, x - 1 / 3)) * 180 / np.pi % 360
+                # correction for multispectral information
+                alpha_corrected = alpha + delta(alpha)
+                cls = int(Raster._Forel_Ule_scale(alpha_corrected))
+            else:
+                cls = int(0)
+
+            # classifying only the valid pixels inside the polygon
+            values = np.where(valid_pixels, cls, 0)
+            # addition to avoid replacing values by cropping by other polygons
+            colour_spt[slices[0], slices[1]] += values.reshape((mask.shape))
+            # classification by polygon
+            colour_shp[i] = cls
+
+        return colour_spt, colour_shp
+
+    @staticmethod
+    def _Forel_Ule_scale(angle):
+        """
+        Calculates the Forel-Ule water colour scale, depending on the hue angle calculated from the Sentinel-2 MSI data,
+        based on the classification by Novoa et al. (2013)
+        Parameters
+        ----------
+        angles: the hue angle, in degrees
+
+        Returns
+        -------
+        The water colour class (1-21)
+        """
+        mapping = [(21.0471, 21), (24.4487, 20), (28.2408, 19), (32.6477, 18), (37.1698, 17), (42.3707, 16),
+                   (47.8847, 15), (53.4431, 14), (59.4234, 13), (64.9378, 12), (70.9617, 11), (78.1648, 10),
+                   (88.5017, 9), (99.5371, 8), (118.5208, 7), (147.4148, 6), (178.7020, 5), (202.8305, 4),
+                   (217.1473, 3), (224.8037, 2)]
+        score = lambda s: next((L for x, L in mapping if s < x), 1)
+
+        return(score(angle))
+
+
+    def chlorophylla(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0):
+        """
+        Function to calculate the chlorophyll-a concentration (chla) based on the optical water type (OWT)
+        The functions are the ones recomended by Carrea et al. (2023) and Neil et al. (2019, 2020), and are coded in
+        inversion_functions.py
+
+        Parameters
+        ----------
+        rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands
+        class_owt_spt: an array, with the same size as the input bands, with the classified pixels
+
+        Returns
+        -------
+        chla: an array, with the same size as the input bands, with the classified pixels
+        """
+        import getpak.inversion_functions as ifunc
+        chla = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float16')
+
+        # chla functions for each OWT
+        classes = [1, 4, 5, 6]
+        index = np.where(np.isin(class_owt_spt, classes))
+        if len(index[0]>0):
+            chla[index] = ifunc.functions['CHL_Gons']['function'](Red=rrs_dict['Rrs_B4'].values[index],
+                                                                          RedEdg1=rrs_dict['Rrs_B5'].values[index],
+                                                                           RedEdg3=rrs_dict['Rrs_B7'].values[index])
+        classes = [2, 7, 8, 11, 12]
+        index = np.where(np.isin(class_owt_spt, classes))
+        if len(index[0] > 0):
+            chla[index] = ifunc.functions['CHL_Gilerson']['function'](Red=rrs_dict['Rrs_B4'].values[index],
+                                                                          RedEdg1=rrs_dict['Rrs_B5'].values[index])
+        classes = [3, 9, 10, 13]
+        index = np.where(np.isin(class_owt_spt, classes))
+        if len(index[0] > 0):
+            chla[index] = ifunc.functions['CHL_OC2']['function'](Blue=rrs_dict['Rrs_B2'].values[index],
+                                                                          Green=rrs_dict['Rrs_B3'].values[index])
+
+        # removing espurious values
+        if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)):
+            chla[chla >= upper_lim] = np.nan
+            chla[chla <= lower_lim] = np.nan
+
+        return chla
+
+    def turbidity(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0):
+        """
+        Function to calculate the turbidity (turb) based on the optical water type (OWT)
+        The functions are the ones recomended by Carrea et al. (2023)
+        Parameters
+        ----------
+        rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands
+        class_owt_spt: an array, with the same size as the input bands, with the classified pixels
+
+        Returns
+        -------
+        turb: an array, with the same size as the input bands, with the classified pixels
+        """
+        import getpak.inversion_functions as ifunc
+        turb = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float16')
+
+        # chla functions for each OWT
+        classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13]
+        index = np.where(np.isin(class_owt_spt, classes))
+        if len(index[0] > 0):
+            turb[index] = ifunc.functions['TURB_Dogliotti']['function'](Red=rrs_dict['Rrs_B4'].values[index],
+                                                                  Nir2=rrs_dict['Rrs_B8A'].values[index])
+
+        # removing espurious values
+        if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)):
+            turb[turb >= upper_lim] = np.nan
+            turb[turb <= lower_lim] = np.nan
+
+        return turb
+
+    def spm(self, rrs_dict, class_owt_spt, upper_lim=5000, lower_lim=0):
+        """
+        Function to calculate the suspended particulate matter (SPM) based on the optical water type (OWT)
+        The functions are the ones recomended by Carrea et al. (2023)
+        Parameters
+        ----------
+        rrs_dict: rrs_dict: a xarray Dataset containing the Rrs bands
+        class_owt_spt: an array, with the same size as the input bands, with the classified pixels
+
+        Returns
+        -------
+        spm: an array, with the same size as the input bands, with the classified pixels
+        """
+        import getpak.inversion_functions as ifunc
+        spm = np.zeros(rrs_dict['Rrs_B4'].shape, dtype='float16')
+
+        # chla functions for each OWT
+        classes = [1, 4, 5, 6, 2, 7, 8, 11, 12, 3, 9, 10, 13]
+        index = np.where(np.isin(class_owt_spt, classes))
+        if len(index[0] > 0):
+            spm[index] = ifunc.functions['SPM_S3']['function'](Red=rrs_dict['Rrs_B4'].values[index],
+                                                                  Nir2=rrs_dict['Rrs_B8A'].values[index])
+
+        # removing espurious values
+        if isinstance(upper_lim, (int, float)) and isinstance(lower_lim, (int, float)):
+            spm[spm >= upper_lim] = np.nan
+            spm[spm <= lower_lim] = np.nan
+
+        return spm
 
 class GRS:
     """
@@ -427,7 +622,7 @@ def metadata(grs_file_entry):
 
     @staticmethod
     def get_grs_dict(grs_nc_file, grs_version='v15'):
-        '''
+        """
         Opens the GRS netCDF files using the xarray library and dask, returning a DataArray containing only the Rrs bands
         Parameters
         ----------
@@ -437,7 +632,7 @@ def get_grs_dict(grs_nc_file, grs_version='v15'):
         Returns
         -------
         The xarray DataArray containing the 11 Rrs bands, named as 'Rrs_B*'
-        '''
+        """
         # list of bands
         bands = ['Rrs_B1', 'Rrs_B2', 'Rrs_B3', 'Rrs_B4', 'Rrs_B5', 'Rrs_B6',
                  'Rrs_B7', 'Rrs_B8', 'Rrs_B8A', 'Rrs_B11', 'Rrs_B12']