Is it possible to mask rendered mosaic tiles using a GeoJSON?

[hrodmn]

I can apply an intersects filter to my mosaic search to limit the mosaic to STAC items that intersect a geometry, but is it be possible to use a GeoJSON to mask pixels in the rendered tiles?

e.g. I want to crop the rendered output so that I can only see pixels on the left side of the purple line:

[vincentsarago]

👋 @hrodmn

That's interesting and TBH I didn't thought about this before but yeah it might make some sense to mask the images when there is a spatial filter in play.

It might be quite tricky to do it natively in titiler-pgstac but rio-tiler/GDAL has ways to do this:

• https://cogeotiff.github.io/rio-tiler/advanced/feature/
• Cropped Tiles developmentseed/titiler#376 (comment)

I think the easiest way to do this will be the apply the mask after tile creation, but it means that you'll have to rasterize the geometry on each tile call, and if the geometry is complex then it might slow down the process.

The other solution is to pass a Cutline to the COG Reader but this might involve more work.

[hrodmn]

Thanks for the suggestions, @vincentsarago! I'll them it out. I thought about it some more last night and for my application it might be simpler to mask the tiles using javascript on the client side, though. I haven't tried it out yet but will let you know what I find!

[hrodmn]

@vincentsarago I worked out a solution for this! At first I tried a geojson approach but did not like how I had to pull the whole geojson for every single tile request. The vector features that I want to use as a mask are all accessible via tipg which makes it possible to use vector tile features that exactly match the raster tiles.

Performance is actually pretty good when running on my Lambda titiler-pgstac deployment (launched using eoapi-cdk fixtures).

Right now I am just setting areas not covered by the mask to zero, but ideally those would get transparent values in the rendered image.
I figured out how to apply the mask to the ImageData.array.mask to properly mask out pixels that are not covered by the vector feature.

"""Custom MosaicTiler Factory for PgSTAC Mosaic Backend.

This module contains a custom MosaicTilerFactory class that extends the default MosaicTiler from titiler-pgstac
(https://github.com/stac-utils/titiler-pgstac/blob/main/titiler/pgstac/factory.py).
The custom factory adds the capability to apply a vector mask to the output image.
The mask is applied by rasterizing a vector tile geometry and using the resulting mask to apply a mask to the image.
"""

import os
from dataclasses import dataclass
from typing import Callable, Dict, Literal, Optional, Type

import httpx
import jinja2
import mapbox_vector_tile
import morecantile
import rasterio
from cogeo_mosaic.backends import BaseBackend
from fastapi import Depends, HTTPException, Query
from pydantic import conint
from rasterio.features import rasterize
from rio_tiler.constants import MAX_THREADS
from rio_tiler.mosaic.methods.base import MosaicMethodBase
from starlette.responses import Response
from starlette.templating import Jinja2Templates
from titiler.core.dependencies import (
    AssetsBidxExprParams,
    ColorFormulaParams,
    DefaultDependency,
    HistogramParams,
    PartFeatureParams,
    StatisticsParams,
    TileParams,
)
from titiler.core.factory import img_endpoint_params
from titiler.core.resources.enums import ImageType, OptionalHeader
from titiler.core.utils import render_image
from titiler.mosaic.factory import PixelSelectionParams
from titiler.pgstac.dependencies import BackendParams, PgSTACParams, TmsTileParams
from titiler.pgstac.factory import MosaicTilerFactory
from titiler.pgstac.mosaic import PGSTACBackend
from typing_extensions import Annotated

MOSAIC_THREADS = int(os.getenv("MOSAIC_CONCURRENCY", MAX_THREADS))
MOSAIC_STRICT_ZOOM = str(os.getenv("MOSAIC_STRICT_ZOOM", False)).lower() in [
    "true",
    "yes",
]

DEFAULT_TEMPLATES = Jinja2Templates(
    directory="",
    loader=jinja2.ChoiceLoader(
        [
            jinja2.PackageLoader(__package__, "templates"),
            jinja2.PackageLoader("titiler.core", "templates"),
        ]
    ),
)  # type:ignore

# define the TIPG endpoint to use for vector tile masks
TIPG_ENDPOINT = "https://tipg-server.com"


def VectorTileFilter(
    tipg_collection: Annotated[
        Optional[str],
        Query(
            description="tipg collection id (schema.table) to query for applying a vector mask, e.g. ncxuser.properties"
        ),
    ] = None,
    tipg_ids: Annotated[
        Optional[str],
        Query(
            description="comma separated list of feature ids from the tipg collection to query for the vector mask, e.g. 101 or 1,2,3"
        ),
    ] = None,
    tipg_filter_lang: Annotated[
        Optional[str],
        Query(
            description="language for the filter parameter (either cql2-text or cql-json)"
        ),
    ] = None,
    tipg_filter: Annotated[
        Optional[str],
        Query(
            description="filter to apply for the tipg query, e.g. property_id = 101 and project_id = 2"
        ),
    ] = None,
    tile=Depends(TmsTileParams),
) -> Optional[str]:
    # if all of the parameters are empty there will not be a vector mask
    if not (tipg_collection or tipg_ids or tipg_filter_lang or tipg_filter):
        return None

    if not tipg_collection:
        raise HTTPException(
            status_code=400,
            detail="You must provide 'tipg_collection' to apply a vector mask",
        )

    if not (tipg_ids or (tipg_filter_lang and tipg_filter)):
        raise HTTPException(
            status_code=400,
            detail="You must provide either 'tipg_ids' or 'tipg_filter_lang' and 'tipg_filter' to apply a vector mask",
        )


    vector_tile_mask_request_url = f"{TIPG_ENDPOINT}/collections/{tipg_collection}/tiles/{tile.z}/{tile.x}/{tile.y}?"

    # add the ids parameter to the request URL if it is provided
    if tipg_ids:
        vector_tile_mask_request_url += f"ids={tipg_ids}"
    # otherwise, add the filter_lang and filter parameters to the request URL
    elif tipg_filter_lang and tipg_filter:
        vector_tile_mask_request_url += (
            f"filter_lang={tipg_filter_lang}&filter={tipg_filter}"
        )
    else:
        raise HTTPException(
            status_code=400,
            detail="You must provide either 'tipg_ids' or 'tipg_filter_lang' and 'tipg_filter' to apply a vector mask",
        )

    return vector_tile_mask_request_url


@dataclass
class MaskableTilerFactory(MosaicTilerFactory):
    """Custom MosaicTiler for PgSTAC Mosaic Backend with capability to apply a vector mask to the output image."""

    path_dependency: Callable[..., str]

    reader: Type[BaseBackend] = PGSTACBackend
    layer_dependency: Type[DefaultDependency] = AssetsBidxExprParams

    # Statistics/Histogram Dependencies
    stats_dependency: Type[DefaultDependency] = StatisticsParams
    histogram_dependency: Type[DefaultDependency] = HistogramParams

    # Tile/Tilejson/WMTS Dependencies
    tile_dependency: Type[DefaultDependency] = TileParams

    # Crop endpoints Dependencies
    # WARNINGS: `/bbox` and `/feature` endpoints should be used carefully because
    # each request might need to open/read a lot of files if the user decide to
    # submit large bbox/geojson. This will also depends on the STAC Items resolution.
    img_part_dependency: Type[DefaultDependency] = PartFeatureParams

    pixel_selection_dependency: Callable[..., MosaicMethodBase] = PixelSelectionParams

    pgstac_dependency: Type[DefaultDependency] = PgSTACParams
    backend_dependency: Type[DefaultDependency] = BackendParams

    # Add/Remove some endpoints
    add_statistics: bool = True
    add_viewer: bool = True
    add_part: bool = True

    templates: Jinja2Templates = DEFAULT_TEMPLATES

    def _tiles_routes(self) -> None:
        """register tiles routes."""

        @self.router.get("/tiles/{z}/{x}/{y}", **img_endpoint_params)
        @self.router.get("/tiles/{z}/{x}/{y}.{format}", **img_endpoint_params)
        @self.router.get("/tiles/{z}/{x}/{y}@{scale}x", **img_endpoint_params)
        @self.router.get("/tiles/{z}/{x}/{y}@{scale}x.{format}", **img_endpoint_params)
        @self.router.get("/tiles/{tileMatrixSetId}/{z}/{x}/{y}", **img_endpoint_params)
        @self.router.get(
            "/tiles/{tileMatrixSetId}/{z}/{x}/{y}.{format}",
            **img_endpoint_params,
        )
        @self.router.get(
            "/tiles/{tileMatrixSetId}/{z}/{x}/{y}@{scale}x",
            **img_endpoint_params,
        )
        @self.router.get(
            "/tiles/{tileMatrixSetId}/{z}/{x}/{y}@{scale}x.{format}",
            **img_endpoint_params,
        )
        def tile(
            search_id=Depends(self.path_dependency),
            tile=Depends(TmsTileParams),
            tileMatrixSetId: Annotated[  # type: ignore
                Literal[tuple(self.supported_tms.list())],
                f"Identifier selecting one of the TileMatrixSetId supported (default: '{self.default_tms}')",
            ] = self.default_tms,
            scale: Annotated[  # type: ignore
                Optional[conint(gt=0, le=4)],
                "Tile size scale. 1=256x256, 2=512x512...",
            ] = None,
            format: Annotated[
                Optional[ImageType],
                "Default will be automatically defined if the output image needs a mask (png) or not (jpeg).",
            ] = None,
            vector_tile_mask_request_url=Depends(VectorTileFilter),
            layer_params=Depends(self.layer_dependency),
            dataset_params=Depends(self.dataset_dependency),
            pixel_selection=Depends(self.pixel_selection_dependency),
            tile_params=Depends(self.tile_dependency),
            post_process=Depends(self.process_dependency),
            rescale=Depends(self.rescale_dependency),
            color_formula=Depends(ColorFormulaParams),
            colormap=Depends(self.colormap_dependency),
            render_params=Depends(self.render_dependency),
            pgstac_params=Depends(self.pgstac_dependency),
            backend_params=Depends(self.backend_dependency),
            reader_params=Depends(self.reader_dependency),
            env=Depends(self.environment_dependency),
        ):
            """Create map tile."""
            scale = scale or 1

            tms = self.supported_tms.get(tileMatrixSetId)
            with rasterio.Env(**env):
                with self.reader(
                    search_id,
                    tms=tms,
                    reader_options={**reader_params},
                    **backend_params,
                ) as src_dst:

                    if MOSAIC_STRICT_ZOOM and (
                        tile.z < src_dst.minzoom or tile.z > src_dst.maxzoom
                    ):
                        raise HTTPException(
                            400,
                            f"Invalid ZOOM level {tile.z}. Should be between {src_dst.minzoom} and {src_dst.maxzoom}",
                        )

                    image, assets = src_dst.tile(
                        tile.x,
                        tile.y,
                        tile.z,
                        tilesize=scale * 256,
                        pixel_selection=pixel_selection,
                        threads=MOSAIC_THREADS,
                        **tile_params,
                        **layer_params,
                        **dataset_params,
                        **pgstac_params,
                        # vrt_options=vrt_options,
                    )

                    if vector_tile_mask_request_url:
                        # fetch the vector tile from the TiPG endpoint
                        response = httpx.get(vector_tile_mask_request_url, timeout=10)
                        response.raise_for_status()

                        # get the bounding box of the tile in Web Mercator coordinates
                        xmin, ymin, xmax, ymax = tms.xy_bounds(
                            morecantile.Tile(tile.x, tile.y, tile.z)
                        )

                        # compute the scale and translation factors
                        scale_x = (xmax - xmin) / 4096.0
                        scale_y = (ymax - ymin) / 4096.0
                        translate_x = xmin
                        translate_y = ymax

                        # decode the vector tile with transformed coordinates
                        tile_data = mapbox_vector_tile.decode(
                            tile=response.content,
                            y_coord_down=True,
                            transformer=lambda x, y: (
                                x * scale_x + translate_x,
                                y * -scale_y + translate_y,
                            ),
                        )

                        # if there aren't any vector features, the mask is empty
                        if not tile_data:
                            raise HTTPException(
                                status_code=204,
                                detail="Vector tile mask is empty",
                            )

                        # rasterize the vector tile geometry
                        vector_mask = rasterize(
                            shapes=[
                                f["geometry"] for f in tile_data["default"]["features"]
                            ],
                            out_shape=(image.height, image.width),
                            transform=image.transform,
                            all_touched=True,
                            default_value=0,  # areas inside the vector geometry will be 0
                            fill=1,  # areas outside the vector geometry will be 1 (to be masked out)
                            dtype="uint8",
                        ).astype(bool)

                        # set image.array.mask to True where image.array.mask is True OR where vector_mask is True
                        image.array.mask = image.array.mask | vector_mask

            if post_process:
                image = post_process(image)

            if rescale:
                image.rescale(rescale)

            if color_formula:
                image.apply_color_formula(color_formula)

            content, media_type = render_image(
                image,
                output_format=format,
                colormap=colormap,
                **render_params,
            )

            headers: Dict[str, str] = {}
            if OptionalHeader.x_assets in self.optional_headers:
                ids = [x["id"] for x in assets]
                headers["X-Assets"] = ",".join(ids)

            return Response(content, media_type=media_type, headers=headers)

[hrodmn]

Here is a view of some raster tiles masked by California, Oregon, and Washington:

[hrodmn]

The answer is, yes, you can mask raster tiles using a geometry source. It is possible to do with a geojson, but that is not ideal since the geojson must be fetched for every tile request!