feat: flood detection algorithm

requirements.txt

@@ -5,4 +5,5 @@ postgis==1.0.4
 uvicorn==0.30.1
 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,
+        "flood_detection": DetectFlood,
     }
 )

src/cogserver/algorithms/flood_detection.py

@@ -0,0 +1,89 @@
+# Credit: Sashka Warner (https://github.com/sashkaw/flood-data-api)
+"""
+MIT License
+
+Copyright (c) 2023 Sashka Warner
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+"""
+
+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
+
+
+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', 'b3']},
+        {'title': 'Short wave infrared band', 'description': 'The SWIR band with wavelength between 0.9μ – 1.7μm',
+         'required': True,
+         'keywords': ['swir', 'b6']},
+    ]
+    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'
+    output_description: str = "The output is a binary image where 1 represents water and 0 represents non-water"
+
+    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()