feat: flood detection algorithm

Dockerfile

@@ -14,7 +14,7 @@ COPY src/cogserver cogserver
 
 ENV HOST=0.0.0.0
 ENV PORT=8000
-ENV LOG_LEVEL=info
+ENV LOG_LEVEL=debug
 ENV RELOAD=--reload
 ENV WORKERS=1
 ENV THREADS=1

requirements.txt

@@ -5,4 +5,5 @@ postgis==1.0.4
 uvicorn==0.29.0
 boto3
 pyyaml
-gunicorn
+gunicorn
+scikit-image

src/cogserver/algorithms/__init__.py

@@ -1,8 +1,10 @@
 from titiler.core.algorithm import Algorithms, algorithms as default_algorithms
 from .rca import RapidChangeAssessment
+from .flood_detection import DetectFlood
 
 algorithms: Algorithms = default_algorithms.register(
     {
-        "rca": RapidChangeAssessment
+        "rca": RapidChangeAssessment,
+        "flooding": DetectFlood,
     }
 )

src/cogserver/algorithms/flood_detection.py

@@ -0,0 +1,65 @@
+from typing import List, Sequence
+
+import numpy as np
+from titiler.core.algorithm import BaseAlgorithm
+from rio_tiler.models import ImageData
+from skimage.filters import threshold_otsu
+
+
+## Credit: Sashka Warner (https://github.com/sashkaw)
+
+class DetectFlood(BaseAlgorithm):
+    title: str = "Flood detection "
+    description: str = "Algorithm to calculate Modified Normalized Difference Water Index (MNDWI), and apply Otsu thresholding algorithm to identify surface water"
+
+    """
+    Desc: Algorithm to calculate Modified Normalized Difference Water Index (MNDWI),
+    and apply Otsu thresholding algorithm to identify surface water.
+    """
+
+    input_bands: List = [
+        {'title': 'Green band', 'description': 'The green band with the wavelength between 0.53µm - 0.59µm',
+         'required': True,
+         'keywords': ['Green band']},
+        {'title': 'Short wave infrared band', 'description': 'The SWIR band with wavelength between 0.9μ – 1.7μm',
+         'required': True,
+         'keywords': ['Shortwave infrared band']},
+    ]
+    input_description: str = "The bands that will be used to make this calculation"
+
+    # Metadata
+    input_nbands: int = 2
+    output_nbands: int = 1
+    output_min: Sequence[int] = [-1]
+    output_max: Sequence[int] = [1]
+    output_colormap_name: str = 'viridis'
+
+    def __call__(self, img: ImageData, *args, **kwargs):
+        # Extract bands of interest
+        green_band = img.data[0].astype("float32")
+        swir_band = img.data[1].astype("float32")
+
+        # Calculate Modified Normalized Difference Water Index (MNDWI)
+        numerator = (green_band - swir_band)
+        denominator = (green_band + swir_band)
+        # Use np.divide to avoid divide by zero errors
+        mndwi_arr = np.divide(numerator, denominator, np.zeros_like(numerator), where=denominator != 0)
+
+        # Apply Otsu thresholding method
+        otsu_threshold = threshold_otsu(mndwi_arr)
+
+        # Use Otsu threshold to classify the computed MNDWI
+        classified_arr = mndwi_arr >= otsu_threshold
+
+        # Reshape data -> ImageData only accepts image in form of (count, height, width)
+        # classified_arr = np.around(classified_arr).astype(int)
+        # classified_arr = np.expand_dims(classified_arr, axis=0).astype(self.output_dtype)
+        classified_arr = np.expand_dims(classified_arr, axis=0).astype(int)
+
+        return ImageData(
+            classified_arr,
+            img.mask,
+            assets=img.assets,
+            crs=img.crs,
+            bounds=img.bounds,
+        )

src/cogserver/dependencies.py

@@ -3,17 +3,17 @@
 from typing import List
 import base64
 
-def parse_signed_url(url:str=None):
 
+def parse_signed_url(url: str = None):
     if '?' in url:
         furl, b64token = url.split('?')
         try:
             decoded_token = base64.b64decode(b64token).decode()
         except Exception:
             decoded_token = b64token
-        decoded_url =  f'{furl}?{decoded_token}'
+        decoded_url = f'{furl}?{decoded_token}'
     else:
-        decoded_url =  f'{url}'
+        decoded_url = f'{url}'
     return decoded_url
 
 
@@ -42,6 +42,7 @@ def SignedDatasetPath(url: Annotated[str, Query(description="Unsigned/signed dat
     """
     return parse_signed_url(url=url)
 
+
 def SignedDatasetPaths(url: Annotated[List[str], Query(description="Unsigned/signed dataset URLs")]) -> str:
     """
         FastAPI dependency function that enables
@@ -61,8 +62,8 @@ def SignedDatasetPaths(url: Annotated[List[str], Query(description="Unsigned/sig
         The returned value is a str representing a RAM stored GDAL VRT file
         which Titiler will use to resolve the request
 
-        Obviously the rasters need to spatially overlap. Additionaly, the VRT can be created with various params
-        (spatial align, resolution, resmapling) that, to some extent can influence the performace of the server
+        Obviously the rasters need to spatially overlap. Additionally, the VRT can be created with various params
+        (spatial align, resolution, resampling) that, to some extent can influence the performance of the server
 
     """
     decoded_urls = list()