diff --git a/app/routes/titiler/algorithms/carbon_gross_emissions.py b/app/routes/titiler/algorithms/carbon_gross_emissions.py new file mode 100644 index 00000000..53d297b4 --- /dev/null +++ b/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) diff --git a/app/routes/titiler/gfw_forest_carbon_gross_emissions.py b/app/routes/titiler/gfw_forest_carbon_gross_emissions.py new file mode 100644 index 00000000..90914a30 --- /dev/null +++ b/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) diff --git a/app/settings/globals.py b/app/settings/globals.py index 87d9de5b..eab00c7d 100644 --- a/app/settings/globals.py +++ b/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]: