GTC-3144 Create carbon emissions titiler endpoint

This is using the colors/symbology from the existing GFW raster tile
cache.

It includes one filter which is tree cover density, using
umd_tree_cover_density_2010/v1.6

app/routes/titiler/algorithms/carbon_gross_emissions.py

@@ -0,0 +1,122 @@
+from collections import OrderedDict, namedtuple
+from typing import Optional
+
+import numpy as np
+from rio_tiler.models import ImageData
+from titiler.core.algorithm import BaseAlgorithm
+
+from app.models.enumerators.titiler import RenderType
+
+Colors: namedtuple = namedtuple("Colors", ["red", "green", "blue"])
+
+
+class CarbonGrossEmissions(BaseAlgorithm):
+    """Visualize carbon gross emissions"""
+
+    title: str = "Carbon gross emissions"
+    description: str = "Visualize carbon gross emissions"
+
+    conf_colors: OrderedDict[float, tuple] = OrderedDict(
+        {
+            0.0: Colors(254, 246, 249),
+            1.0: Colors(245, 237, 242),
+            6.0: Colors(236, 228, 236),
+            14.0: Colors(227, 220, 231),
+            24.0: Colors(217, 210, 228),
+            39.0: Colors(209, 201, 227),
+            56.0: Colors(202, 192, 228),
+            76.0: Colors(194, 183, 229),
+            99.0: Colors(187, 173, 231),
+            126.0: Colors(182, 164, 232),
+            156.0: Colors(177, 154, 231),
+            188.0: Colors(173, 145, 229),
+            224.0: Colors(169, 134, 225),
+            263.0: Colors(165, 126, 218),
+            305.0: Colors(162, 117, 211),
+            351.0: Colors(158, 109, 202),
+            399.0: Colors(153, 100, 191),
+            451.0: Colors(149, 93, 181),
+            505.0: Colors(144, 86, 171),
+            563.0: Colors(140, 79, 160),
+            624.0: Colors(134, 71, 148),
+            688.0: Colors(128, 65, 138),
+            755.0: Colors(123, 59, 127),
+            825.0: Colors(116, 53, 117),
+            899.0: Colors(109, 47, 105),
+            975.0: Colors(102, 42, 95),
+            1055.0: Colors(95, 37, 85),
+            1137.0: Colors(88, 32, 76),
+            1223.0: Colors(80, 26, 66),
+            1312.0: Colors(72, 21, 57),
+            1404.0: Colors(64, 15, 50),
+            1500.0: Colors(57, 8, 42)
+        }
+    )
+
+    tree_cover_density_mask: Optional[int] = None
+    tree_cover_density_data: Optional[ImageData] = None
+    render_type: RenderType = RenderType.encoded
+
+    # metadata
+    input_nbands: int = 2
+    output_nbands: int = 4
+    output_dtype: str = "uint8"
+
+    def __call__(self, img: ImageData) -> ImageData:
+
+        self.emissions = img.data[0]
+        self.intensity = img.data[1]
+        self.no_data = img.array.mask[0]
+
+        self.mask = self.create_mask()
+
+        if self.render_type == RenderType.true_color:
+            rgb = self.create_true_color_rgb()
+            alpha = self.create_true_color_alpha()
+        # else:  # encoded
+        #    rgb = self.create_encoded_rgb()
+        #    alpha = self.create_encoded_alpha()
+
+        data = np.vstack([rgb, alpha[np.newaxis, ...]]).astype(self.output_dtype)
+        data = np.ma.MaskedArray(data, mask=False)
+
+        return ImageData(data, assets=img.assets, crs=img.crs, bounds=img.bounds)
+
+    def create_mask(self):
+        mask = ~self.no_data
+
+        if self.tree_cover_density_mask:
+            mask *= (
+                self.tree_cover_density_data.array[0, :, :]
+                >= self.tree_cover_density_mask
+            )
+
+        return mask
+
+    # def create_encoded_rgb(self):
+    #    r, g, b = self._rgb_zeros_array()
+
+    def create_true_color_rgb(self):
+        r, g, b = self._rgb_zeros_array()
+
+        for k, v in reversed(self.conf_colors.items()):
+            if self.emissions >= k:
+                r = v.red
+                g = v.green
+                b = v.blue
+                break
+
+        return np.stack([r, g, b], axis=0)
+
+    def create_true_color_alpha(self):
+        """Set the transparency (alpha) channel based on intensity input. The
+        intensity multiplier is used to control how isolated pixels fade out at low
+        zoom levels, matching the rendering behavior in Flagship.
+
+        Returns:
+            np.ndarray: Array representing the alpha (transparency) channel, where pixel
+            visibility is adjusted by intensity.
+
+        """
+        alpha = np.where(self.mask, self.intensity * 150, 0)
+        return np.minimum(255, alpha)

app/routes/titiler/gfw_forest_carbon_gross_emissions.py

@@ -0,0 +1,81 @@
+import os
+from typing import Optional, Tuple
+
+from aenum import Enum, extend_enum
+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 ...crud.sync_db.tile_cache_assets import get_versions
+from ...models.enumerators.tile_caches import TileCacheType
+from ...models.enumerators.titiler import RenderType
+from .. import raster_xyz
+from ...settings.globals import GLOBALS
+from .algorithms.carbon_gross_emissions import CarbonGrossEmissions
+from .readers import AlertsReader
+
+DATA_LAKE_BUCKET = os.environ.get("DATA_LAKE_BUCKET")
+
+router = APIRouter()
+
+dataset = "gfw_forest_carbon_gross_emissions"
+
+
+class GfwForestCarbonGrossEmissions(str, Enum):
+    latest = "v20240402"
+
+
+_versions = get_versions(dataset, TileCacheType.cog)
+for _version in _versions:
+    extend_enum(GfwForestCarbonGrossEmissions, _version, _version)
+
+
+@router.get(
+    f"/{dataset}/{{version}}/dynamic/{{z}}/{{x}}/{{y}}.png",
+    response_class=Response,
+    tags=["Raster Tiles"],
+    response_description="PNG Raster Tile",
+)
+async def global_forest_carbon_gross_emissions_raster_tile(
+    *,
+    version: GfwForestCarbonGrossEmissions,
+    xyz: Tuple[int, int, int] = Depends(raster_xyz),
+    tree_cover_density_threshold: Optional[int] = Query(
+        None,
+        ge=0,
+        le=100,
+        description="Show alerts in pixels with tree cover density (in percent) greater than or equal to this threshold. `umd_tree_cover_density_2010` is used for this masking.",
+    ),
+    render_type: RenderType = Query(RenderType.encoded, description="Render true color or encoded tiles")
+) -> Response:
+    """Forest Carbon Gross Emissions 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:
+        # 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)
+
+    carbon_gross_emissions = CarbonGrossEmissions(
+        tree_cover_density_mask=tree_cover_density_threshold,
+        render_type=render_type,
+    )(image_data)
+
+    filter_datasets = GLOBALS.carbon_gross_emissions_filters
+    if tree_cover_density_threshold:
+        filter_dataset = filter_datasets["tree_cover_density"]
+        with COGReader(
+            f"s3://{DATA_LAKE_BUCKET}/{filter_dataset['dataset']}/{filter_dataset['version']}/raster/epsg-4326/cog/default.tif"
+        ) as reader:
+            carbon_gross_emissions.tree_cover_density_data = reader.tile(tile_x, tile_y, zoom)
+
+    content, media_type = render_image(
+        carbon_gross_emissions,
+        output_format=ImageType("png"),
+        add_mask=False,
+    )
+
+    return Response(content, media_type=media_type)

app/settings/globals.py

@@ -16,6 +16,10 @@
     "tree_cover_height": {"dataset": "umd_tree_cover_height_2020", "version": "v2022"},
     "tree_cover_density": {"dataset": "umd_tree_cover_density_2010", "version": "v1.6"},
 }
+carbon_gross_emissions_filters = {
+    "tree_cover_density": {"dataset": "umd_tree_cover_density_2010", "version": "v1.6"},
+    # "tree_cover_density": {"dataset": "umd_tree_cover_density_2000", "version": "v1.8"},
+}
 
 
 class Globals(BaseSettings):
@@ -80,6 +84,10 @@ class Globals(BaseSettings):
         dist_alerts_forest_filters,
         description="Datasets that are used as forest filters for DIST alerts",
     )
+    carbon_gross_emissions_filters: Dict = Field(
+        carbon_gross_emissions_filters,
+        description="Datasets that are used as filters for carbon gross emissions"
+    )
 
     @field_validator("token", mode="before")
     def get_token(cls, v: Optional[str]) -> Optional[str]: