Merge pull request #170 from wri/gtc-3025/filter-alerts-by-forest

Gtc 3025/filter alerts by forest

app/routes/titiler/algorithms/dist_alerts.py

@@ -1,4 +1,8 @@
 from collections import OrderedDict
+from typing import Optional
+
+from pydantic import ConfigDict
+from rio_tiler.models import ImageData
 
 from app.models.enumerators.titiler import AlertConfidence
 
@@ -9,6 +13,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 +26,41 @@ 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
+
+    # the highest loss year that is used to exclude alerts for
+    # the purpose of showing only alerts in forests
+    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:
+            # Tree cover loss data before 2020 can't be used to filter out pixels as not forest.
+            # Instead, we use tree cover height taken that year as source of truth.
+            # For example, if a pixel had tree cover loss in 2018, but has tree cover
+            # height (2020) that meets the forest threshold, the pixel meets
+            # the forest criteria for alerts and is not masked out.
+            mask *= (
+                self.tree_cover_loss_data.array[0, :, :] > self.tree_cover_loss_mask
+            ) | (self.tree_cover_loss_data.array[0, :, :] <= 2020)
+
+        return mask

app/routes/titiler/umd_glad_dist_alerts.py

@@ -4,10 +4,12 @@
 
 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
 
 from ...models.enumerators.titiler import AlertConfidence, RenderType
+from ...settings.globals import GLOBALS
 from .. import DATE_REGEX, raster_xyz
 from .algorithms.dist_alerts import DISTAlerts
 from .readers import AlertsReader
@@ -17,13 +19,13 @@
 router = APIRouter()
 
 # TODO: update to the actual dataset when ready
-dataset = "dan_test"
+DATASET = "dan_test"
 
 today = date.today()
 
 
 @router.get(
-    f"/{dataset}/{{version}}/dynamic/{{z}}/{{x}}/{{y}}.png",
+    f"/{DATASET}/{{version}}/dynamic/{{z}}/{{x}}/{{y}}.png",
     response_class=Response,
     tags=["Raster Tiles"],
     response_description="PNG Raster Tile",
@@ -49,24 +51,65 @@ 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 in conjunction with `tree_cover_density` and `tree_cover_height` filters to detect only alerts in forests, by masking out pixels that have had tree cover loss prior to the alert.""",
+    ),
 ) -> 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"
+    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,
+    )
+
+    filter_datasets = GLOBALS.dist_alerts_forest_filters
+    if tree_cover_density:
+        dataset = filter_datasets["tree_cover_density"]
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/{dataset['dataset']}/{dataset['version']}/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:
+        dataset = filter_datasets["tree_cover_height"]
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/{dataset['dataset']}/{dataset['version']}/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:
+        dataset = filter_datasets["tree_cover_loss"]
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/{dataset['dataset']}/{dataset['version']}/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,

app/settings/globals.py

@@ -1,14 +1,22 @@
 import json
 from json import JSONDecodeError
-from typing import Optional
+from typing import Dict, Optional
 
 from fastapi.logger import logger
 from pydantic import Field, field_validator
-from pydantic_settings import SettingsConfigDict, BaseSettings
+from pydantic_settings import BaseSettings, SettingsConfigDict
 from starlette.datastructures import Secret
 
 from ..models.pydantic.database import DatabaseURL
 
+# NOTE: For non-local deployment, need to update dataset/version here until the terraform
+# JSON decode issue preventing setting this env variable is resolved.
+dist_alerts_forest_filters = {
+    "tree_cover_loss": {"dataset": "umd_tree_cover_loss", "version": "v1.10.1"},
+    "tree_cover_height": {"dataset": "umd_tree_cover_height_2020", "version": "v2022"},
+    "tree_cover_density": {"dataset": "umd_tree_cover_density_2010", "version": "v1.6"},
+}
+
 
 class Globals(BaseSettings):
     env: str = Field("dev", description="Environment name.")
@@ -68,6 +76,11 @@ class Globals(BaseSettings):
     )
     api_key_name: str = Field("x-api-key", description="Header key name for API key.")
 
+    dist_alerts_forest_filters: Dict = Field(
+        dist_alerts_forest_filters,
+        description="Datasets that are used as forest filters for DIST alerts",
+    )
+
     @field_validator("token", mode="before")
     def get_token(cls, v: Optional[str]) -> Optional[str]:
         if v:
@@ -122,6 +135,15 @@ def set_lambda_host(cls, v, values, **kwargs) -> str:
             v = f"https://lambda.{aws_region}.amazonaws.com"
         return v
 
+    @field_validator("dist_alerts_forest_filters", mode="before")
+    def parse_dist_alerts_forest_filters(cls, v: str | Dict) -> Dict:
+        if isinstance(v, dict):
+            return v
+        try:
+            return json.loads(v)
+        except json.JSONDecodeError as e:
+            raise ValueError(f"Invalid JSON for dist_alerts_forest_filters: {e}")
+
     model_config = SettingsConfigDict(case_sensitive=False, validate_assignment=True)
 
 

docker-compose.dev.yml

@@ -26,6 +26,7 @@ services:
       - ENV=dev
       - PLANET_API_KEY
       - AWS_REQUEST_PAYER=requester
+      - 'DIST_ALERTS_FOREST_FILTERS={"tree_cover_loss": {"dataset": "umd_tree_cover_loss", "version": "v1.10.1"}, "tree_cover_height": {"dataset": "umd_tree_cover_height_2020", "version": "v2022"}, "tree_cover_density": {"dataset": "umd_tree_cover_density_2010", "version": "v1.6"}}'
     ports:
       - 8088:80
     entrypoint: wait_for_postgres.sh /start-reload.sh

tests/titiler/test_alerts_algo.py

@@ -5,6 +5,7 @@
 from rio_tiler.models import ImageData
 
 from app.models.enumerators.titiler import IntegratedAlertConfidence, RenderType
+from app.routes.titiler.algorithms.dist_alerts import DISTAlerts
 from app.routes.titiler.algorithms.integrated_alerts import IntegratedAlerts
 from tests.conftest import DATE_CONF_TIF, INTENSITY_TIF
 
@@ -24,6 +25,50 @@ def get_tile_data():
     return ImageData(data)
 
 
+def get_tcl_data():
+    """Tree Cover Loss test data."""
+    with rasterio.open(DATE_CONF_TIF) as date_conf_file:
+        date_conf = date_conf_file.read(1)
+
+    data = np.zeros_like(date_conf)
+
+    data[122, 109] = 2019
+    data[120, 109] = 2022
+    data[154, 71] = 2023
+
+    return ImageData(data)
+
+
+def get_tch_data():
+    """Tree Cover Height test data."""
+
+    with rasterio.open(DATE_CONF_TIF) as date_conf_file:
+        date_conf = date_conf_file.read(1)
+
+    data = np.zeros_like(date_conf)
+
+    data[122, 109] = 2
+    data[120, 109] = 5
+    data[154, 71] = 4
+
+    return ImageData(data)
+
+
+def get_tcd_data():
+    """Tree Cover Density test data."""
+
+    with rasterio.open(DATE_CONF_TIF) as date_conf_file:
+        date_conf = date_conf_file.read(1)
+
+    data = np.zeros_like(date_conf)
+
+    data[122, 109] = 40
+    data[120, 109] = 30
+    data[154, 71] = 20
+
+    return ImageData(data)
+
+
 def test_integrated_alerts_defaults():
     """Test default values of the Alerts class."""
     alerts = IntegratedAlerts()
@@ -104,3 +149,33 @@ def test_encoded_rgba():
 
     # test high confidence in alpha channel
     assert rgba.array[3, 154, 71] == 8
+
+
+def test_forest_mask():
+    """Test that only alerts that match forest criteria are shown."""
+
+    alerts_data = get_tile_data()
+    tcl = get_tcl_data()
+    tch = get_tch_data()
+    tcd = get_tcd_data()
+
+    alerts = DISTAlerts(
+        tree_cover_density_mask=30,
+        tree_cover_height_mask=3,
+        tree_cover_loss_mask=2021,
+        render_type=RenderType.true_color,
+    )
+
+    alerts.tree_cover_density_data = tcd
+    alerts.tree_cover_height_data = tch
+    alerts.tree_cover_loss_data = tcl
+
+    rgba = alerts(alerts_data)
+
+    np.testing.assert_array_equal(rgba.array[:, 122, 109], np.array([220, 102, 153, 0]))
+
+    np.testing.assert_array_equal(rgba.array[:, 154, 71], np.array([220, 102, 153, 0]))
+
+    np.testing.assert_array_equal(
+        rgba.array[:, 120, 109], np.array([220, 102, 153, 255])
+    )