Custom TMS for MosaicJSON endpoints

[kristhjo]

Hi!

First of all, thanks for the awesome work done on titiler and the related framework.

I'm investigating titiler for a project where we want to stitch together numerous COGs as a MosaicJSON. For the standard WebMercatorQuad this works great, but we also need to be able to support other TileMatrixSets. It seems like the api spec is designed for this to be possible, but that it's not implemented at the moment? If so, I would greatly appreciate any thoughts on how I might proceed on this. I guess it should be possible to follow the COG tiler as an example, but that there might be some nuances to how it's implemented in the cogeo_mosaic backend (?). Anyways, any hints poking me in the right direction or identifying possible pitfalls are welcome!

[vincentsarago]

@kristhjo Sadly create MosaicJSON in other TMS is quite complex (cogeo-mosaic won't support it because of https://github.com/developmentseed/cogeo-mosaic/blob/main/cogeo_mosaic/mosaic.py#L138-L145).

As for the Reader part, none of the default Backends will support non WebMercator TMS https://github.com/developmentseed/cogeo-mosaic/blob/main/cogeo_mosaic/backends/base.py#L64

The easiest way to use non Webmercator for Mosaic is to use titiler-pgstact (https://github.com/stac-utils/titiler-pgstac/blob/master/titiler/pgstac/mosaic.py#L93-L103) or to create a Custom Backend (if you find a way to create a MosaicJSON)

[kristhjo]

Thanks for the feedback and the great example! If I understand correctly then, the problem for cogeo_mosaic is essentially that supermercado's burntiles does not support other TMS definitions?

I will definitely check out the titiler-pgstac repo, this looks really interesting. Do you have any benchmarks or rough estimates on the difference in tile creation timings between static and dynamic mosaics though? I would expect the dynamic mosaics to be quite a bit slower as the get_assets for tile intersection needs to be calculated for every request?

[vincentsarago]

the problem for cogeo_mosaic is essentially that supermercado's burntiles does not support other TMS definitions?

Yes

Do you have any benchmarks or rough estimates on the difference in tile creation timings between static and dynamic mosaics though?

We don't have benchmark, and to be honest the value of dynamic mosaic (multiple projection, no more static files, update, ...) overcome the speed of static mosaic.

I would expect the dynamic mosaics to be quite a bit slower as the get_assets for tile intersection needs to be calculated for every request?

Depends what backend you're using but if there is a spatial index it's usually pretty fast.

[vincentsarago]

👇 is an example on how to create a Dynamic backend using GeoPandas (the mosaicJSON is replaced with a FlatGeoBuff file)

from typing import Any, Dict, List, Callable, Union

import attr
from morecantile import TileMatrixSet
from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
from rio_tiler.io import BaseReader, Reader

from cogeo_mosaic.backends.base import BaseBackend
from cogeo_mosaic.mosaic import MosaicJSON

from shapely.geometry import Polygon, Point
import geopandas as gpd
import rasterio

from pyproj import CRS, Transformer


@attr.s
class GeoPandasBackend(BaseBackend):

    input: str = attr.ib()

    tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)

    reader: BaseReader = attr.ib(default=Reader)
    reader_options: Dict = attr.ib(factory=dict)

    geographic_crs: CRS = attr.ib(default=WGS84_CRS)

    minzoom: int = attr.ib()
    maxzoom: int = attr.ib()

    # The reader is read-only, we can't pass mosaic_def to the init method
    mosaic_def: MosaicJSON = attr.ib(init=False)

    dataset: Any = attr.ib(init=False)

    _transform: Transformer = attr.ib(init=False)
    _from_wgs84: Transformer = attr.ib(init=False)

    _backend_name = "GeoPandas"

    @minzoom.default
    def _minzoom(self):
        return self.tms.minzoom

    @maxzoom.default
    def _maxzoom(self):
        return self.tms.maxzoom

    def __attrs_post_init__(self):
        """Post Init."""
        # Construct a FAKE mosaicJSON
        # mosaic_def has to be defined. As we do for the DynamoDB and SQLite backend
        # we set `tiles` to an empty list.
        self.mosaic_def = MosaicJSON(
            mosaicjson="0.0.2",
            name="it's fake but it's ok",
            minzoom=self.minzoom,
            maxzoom=self.maxzoom,
            tiles=[]
        )

        self.dataset = gpd.read_file(self.input)

        # Get Bounds
        self.bounds = tuple(self.dataset.total_bounds)

        # Get CRS
        self.crs = rasterio.crs.CRS.from_wkt(self.dataset.crs.to_wkt())

        self._transform = Transformer.from_crs(self.tms.crs, self.dataset.crs, always_xy=True)
        self._from_wgs84 = Transformer.from_crs(CRS.from_epsg(4326), self.dataset.crs, always_xy=True)

    def close(self):
        """Close all dataset."""
        pass

    def __exit__(self, exc_type, exc_value, traceback):
        """Support using with Context Managers."""
        self.close()

    def write(self, overwrite: bool = True):
        """This method is not used but is required by the abstract class."""
        pass

    def update(self):
        """We overwrite the default method."""
        pass

    def _read(self) -> MosaicJSON:
        """This method is not used but is required by the abstract class."""
        pass

    def assets_for_tile(self, x: int, y: int, z: int) -> List[str]:
        """Retrieve assets for tile."""
        bbox = self.tms.xy_bounds(x, y, z)  # BBOX in TMS coordinates
        bbox = self._transform.transform_bounds(*bbox, densify_pts=21)  # reproject bbox to dataset CRS
        return self.get_assets(Polygon.from_bounds(*bbox))

    def assets_for_point(self, lng: float, lat: float) -> List[str]:
        """Retrieve assets for point."""
        x, y = self._from_wgs84.transform(lng, lat)
        return self.get_assets(Point(x, y))

    def get_assets(self, geom: Union[Polygon, Point]) -> List[str]:
        """Find assets."""
        return [
            d.path for _, d in self.dataset.loc[self.dataset.intersects(geom)].iterrows()
        ]

    @property
    def _quadkeys(self) -> List[str]:
        return []

Example of usage with latest Open Imagery for Ian Hurricane

# Create a GeoJSON of all the COGs
# 1. list all the files
# 2. use cogeo-mosaic to get the geometry for each file
# 3. remove useless properties
$ aws s3 ls s3://noaa-eri-pds/2022_Hurricane_Ian/ --recursive  | grep ".tif" | awk '{print "s3://noaa-eri-pds/"$NF}' | cogeo-mosaic footprint - --threads 20 | jq -c '.features[] | del(.properties.bounds,.properties.datatype)' > 2022_Hurricane_Ian_features.geojsonseq

# Create a FlatGeoBuff
$ ogr2ogr -f FlatGeobuf 2022_Hurricane_Ian_features.fgb 2022_Hurricane_Ian_features.geojsonseq

import morecantile
tms = morecantile.tms.get("WGS1984Quad")

with GeoPandasBackend("2022_Hurricane_Ian_features.fgb") as mosaic:
    assets = mosaic.assets_for_point(-79.26, 33.10)
>>> ['s3://noaa-eri-pds/2022_Hurricane_Ian/20221003a_RGB/20221003aC0791545w330645n.tif']

# Using Default WebMercatorQuad TMS
with GeoPandasBackend("2022_Hurricane_Ian_features.fgb") as mosaic:
    assets = mosaic.assets_for_tile(4469, 6931, 14)
>>> ['s3://noaa-eri-pds/2022_Hurricane_Ian/20220929a_RGB/20220929aC0814845w264030n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930b_RGB/20220930bC0814845w264030n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930d_RGB/20220930dC0814845w264030n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930d_RGB/20220930dC0814845w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930b_RGB/20220930bC0814845w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220929a_RGB/20220929aC0814845w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220929a_RGB/20220929aC0814800w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930d_RGB/20220930dC0814800w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930b_RGB/20220930bC0814800w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930b_RGB/20220930bC0814715w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930d_RGB/20220930dC0814715w264115n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930d_RGB/20220930dC0814715w264030n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930b_RGB/20220930bC0814715w264030n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930d_RGB/20220930dC0814800w264030n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20220930b_RGB/20220930bC0814800w264030n.tif'] 

# Using WGS1984Quad TMS
with GeoPandasBackend("2022_Hurricane_Ian_features.fgb", tms=tms) as mosaic:
    assets = mosaic.assets_for_tile(9169, 5179, 14)
>>> ['s3://noaa-eri-pds/2022_Hurricane_Ian/20221003a_RGB/20221003aC0791545w330645n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20221002d_RGB/20221002dC0791630w330645n.tif',
 's3://noaa-eri-pds/2022_Hurricane_Ian/20221002d_RGB/20221002dC0791630w330600n.tif']