add forest filter layers

app/routes/titiler/algorithms/dist_alerts.py

@@ -1,4 +1,9 @@
 from collections import OrderedDict
+from datetime import datetime
+from typing import Optional
+
+from pydantic import ConfigDict
+from rio_tiler.models import ImageData
 
 from app.models.enumerators.titiler import AlertConfidence
 
@@ -9,6 +14,7 @@ class DISTAlerts(Alerts):
     title: str = "Land Disturbance (DIST) Alerts"
     description: str = "Decode and visualize DIST alerts"
 
+    model_config = ConfigDict(arbitrary_types_allowed=True)
     conf_colors: OrderedDict = OrderedDict(
         {
             AlertConfidence.low: AlertConfig(
@@ -21,3 +27,37 @@ class DISTAlerts(Alerts):
     )
 
     record_start_date: str = "2020-12-31"
+
+    tree_cover_density_mask: Optional[int] = None
+    tree_cover_density_data: Optional[ImageData] = None
+
+    tree_cover_height_mask: Optional[int] = None
+    tree_cover_height_data: Optional[ImageData] = None
+
+    tree_cover_loss_mask: Optional[int] = None
+    tree_cover_loss_data: Optional[ImageData] = None
+
+    def create_mask(self):
+        mask = super().create_mask()
+
+        if self.tree_cover_density_mask:
+            mask *= (
+                self.tree_cover_density_data.array[0, :, :]
+                >= self.tree_cover_density_mask
+            )
+
+        if self.tree_cover_height_mask:
+            mask *= (
+                self.tree_cover_height_data.array[0, :, :]
+                >= self.tree_cover_height_mask
+            )
+
+        if self.tree_cover_loss_mask:
+            mask *= (
+                self.tree_cover_loss_data.array[0, :, :] >= self.tree_cover_loss_mask
+                or self.tree_cover_loss_data.array[0, :, :] == 0
+                or self.tree_cover_loss_data.array[0, :, :]
+                <= datetime.strptime(self.record_start_date, "%Y-%m-%d").year
+            )
+
+        return mask

app/routes/titiler/umd_glad_dist_alerts.py

@@ -4,6 +4,7 @@
 
 from dateutil.relativedelta import relativedelta
 from fastapi import APIRouter, Depends, Query, Response
+from rio_tiler.io import COGReader
 from titiler.core.resources.enums import ImageType
 from titiler.core.utils import render_image
 
@@ -49,24 +50,62 @@ async def glad_dist_alerts_raster_tile(
         AlertConfidence.low,
         description="Show alerts that are at least of this confidence level",
     ),
+    tree_cover_density: Optional[int] = Query(
+        None,
+        ge=0,
+        le=100,
+        description="Alerts in pixels with tree cover density (in percent) below this threshold won't be displayed. `umd_tree_cover_density_2010` is used for this masking.",
+    ),
+    tree_cover_height: Optional[int] = Query(
+        None,
+        description="Alerts in pixels with tree cover height (in meters) below this threshold won't be displayed. `umd_tree_cover_height_2020` dataset in the API is used for this masking.",
+    ),
+    tree_cover_loss_cutoff: bool = Query(
+        False,
+        ge=2021,
+        description="""This filter is to be used on conjunction with `tree_cover_density` and `tree_cover_height` filters to detect only alerts in forests. """
+        """Alerts for pixels that have had tree cover loss this year or earlier (to 2021) won't be displayed.""",
+    ),
 ) -> Response:
     """UMD GLAD DIST alerts raster tiles."""
 
+    tile_x, tile_y, zoom = xyz
     bands = ["default", "intensity"]
     folder: str = f"s3://{DATA_LAKE_BUCKET}/{dataset}/{version}/raster/epsg-4326/cog"
     with AlertsReader(input=folder) as reader:
-        tile_x, tile_y, zoom = xyz
-
         # NOTE: the bands in the output `image_data` array will be in the order of
         # the input `bands` list
         image_data = reader.tile(tile_x, tile_y, zoom, bands=bands)
 
-    processed_image = DISTAlerts(
+    dist_alert = DISTAlerts(
         start_date=start_date,
         end_date=end_date,
         render_type=render_type,
         alert_confidence=alert_confidence,
-    )(image_data)
+        tree_cover_density_mask=tree_cover_density,
+        tree_cover_height_mask=tree_cover_height,
+        tree_cover_loss_mask=tree_cover_loss_cutoff,
+    )
+
+    if tree_cover_density:
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/umd_tree_cover_density_2010/v1.6/raster/epsg-4326/cog/default.tif"
+        ) as reader:
+            dist_alert.tree_cover_density_data = reader.tile(tile_x, tile_y, zoom)
+
+    if tree_cover_height:
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/umd_tree_cover_height_2020/v2022/raster/epsg-4326/cog/default.tif"
+        ) as reader:
+            dist_alert.tree_cover_height_data = reader.tile(tile_x, tile_y, zoom)
+
+    if tree_cover_loss_cutoff:
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/umd_tree_cover_loss/v1.10.1/raster/epsg-4326/cog/default.tif"
+        ) as reader:
+            dist_alert.tree_cover_loss_data = reader.tile(tile_x, tile_y, zoom)
+
+    processed_image = dist_alert(image_data)
 
     content, media_type = render_image(
         processed_image,