Adding support for non-WebMercatorQuad TMS

[jlaura]

We recently added a large number of TMS to the planetcantile, that we have made available to the community at https://astrogeology.usgs.gov/titiler. In this way, we have planetary data support for COGs.

We just started testing mosaic.json and realized rather quickly that we need to have support for more than 3857. I see multiple places in the code base where comments indicate that only 3857 is currently supported and wanted to start a discussion a out how to go about adding support for (many) other TMSs.

[vincentsarago]

ref:

• Support other projection mosaicjson-spec#13
• Custom TMS for MosaicJSON endpoints titiler#518

The main issues are:

• how to reference the TMS you are using within the mosaicJSON (Support other projection mosaicjson-spec#13)
• how to advertise that a mosaic backend will only support one TMS (in the tiler)

Sadly I'm not really optimistic about the time I could spend on this :-(

[vincentsarago]

there looks to be a few layers here but would updating supermercado to support arbitrary TMSs be a step in the right direction?

oh I mis understood.

cogeo-mosaic is used to create mosaicJSON files and provides backends (used to read/write MosaicJSON into/from a datastore). You are right to create MosaicJSON we would need a supermorecado like version of supermercado (just burntiles submodule). It's unlikely that supermercado will be updated so we better create our own version in cogeo-mosaic.

As mentioned there are a lot of layers to touch in order to support non-webmercator TMS for MosaicJSON, the first one would be in MosaicJSON specification.

[AndrewAnnex]

well I may have made some progress there, no tests and assumes json input are in the crs of the TMS but it's a start https://github.com/AndrewAnnex/supermercado/tree/supermorecantile.

[vincentsarago]

nice

I was about to push this myself

"""Supermercado.burntiles but for other TMS."""

from typing import Dict, List, Tuple
import attr

import numpy as np
import morecantile
from affine import Affine

from rasterio import features
from rio_tiler.constants import WEB_MERCATOR_TMS


@attr.s
class burntiles:

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

    ## TODO Add check for quadkey support

    def project_geom(self, geom: Dict) -> Dict:
        if geom["type"] == "Polygon":
            return {
                "type": geom["type"],
                "coordinates": [
                    [self.tms.xy(*coords) for coords in part]
                    for part in geom["coordinates"]
                ],
            }

        elif geom["type"] == "LineString":
            return {
                "type": geom["type"],
                "coordinates": [self.tms.xy(*coords) for coords in geom["coordinates"]],
            }

        elif geom["type"] == "Point":
            return {
                "type": geom["type"],
                "coordinates": self.tms.xy(*geom["coordinates"]),
            }

    def _feature_extrema(self, geometry: Dict) -> Tuple[float, float, float, float]:
        if geometry["type"] == "Polygon":
            x, y = zip(*[c for part in geometry["coordinates"] for c in part])

        elif geometry["type"] == "LineString":
            x, y = zip(*[c for c in geometry["coordinates"]])

        elif geometry["type"] == "Point":
            x, y = geometry["coordinates"]
            return x, y, x, y

        return min(x), min(y), max(x), max(y)

    def find_extrema(self, features: List[Dict]) -> Tuple[float, float, float, float]:
        epsilon = 1.0e-10
        min_x, min_y, max_x, max_y = zip(
            *[self._feature_extrema(f["geometry"]) for f in features]
        )

        return (
            min(min_x) + epsilon,
            max(min(min_y) + epsilon, -85.0511287798066),
            max(max_x) - epsilon,
            min(max(max_y) - epsilon, 85.0511287798066),
        )


    def tile_extrema(self, bounds, zoom):
        minimumTile = self.tms.tile(bounds[0], bounds[3], zoom)
        maximumTile = self.tms.tile(bounds[2], bounds[1], zoom)

        return {
            "x": {"min": minimumTile.x, "max": maximumTile.x + 1},
            "y": {"min": minimumTile.y, "max": maximumTile.y + 1},
        }


    def make_transform(self, tilerange: Dict, zoom: int) -> Affine:

        ulx, uly = self.tms.xy(
            *self.tms.ul(tilerange["x"]["min"], tilerange["y"]["min"], zoom)
        )

        lrx, lry = self.tms.xy(
            *self.tms.ul(tilerange["x"]["max"], tilerange["y"]["max"], zoom)
        )

        xcell = (lrx - ulx) / float(tilerange["x"]["max"] - tilerange["x"]["min"])
        ycell = (uly - lry) / float(tilerange["y"]["max"] - tilerange["y"]["min"])

        return Affine(xcell, 0, ulx, 0, -ycell, uly)


    def burn(self, polys: List[Dict], zoom: int):
        bounds = self.find_extrema(polys)

        tilerange = self.tile_extrema(bounds, zoom)
        afftrans = self.make_transform(tilerange, zoom)

        burn = features.rasterize(
            ((self.project_geom(geom["geometry"]), 255) for geom in polys),
            out_shape=(
                (
                    tilerange["y"]["max"] - tilerange["y"]["min"],
                    tilerange["x"]["max"] - tilerange["x"]["min"],
                )
            ),
            transform=afftrans,
            all_touched=True,
        )

        xys = np.fliplr(np.dstack(np.where(burn))[0])

        xys[:, 0] += tilerange["x"]["min"]
        xys[:, 1] += tilerange["y"]["min"]

        return np.append(xys, np.zeros((xys.shape[0], 1), dtype=np.uint8) + zoom, axis=1)

(because we just need supermercado.burntiles.burn)

[AndrewAnnex]

that works too I think, but I had to do some other updates to find_extrema for the extend bounds

[vincentsarago]

@AndrewAnnex can you check the PR and see if it works?
I'll push some tests soon

[AndrewAnnex]

Hey Vincent I should be able to in an hour or two - Dr. Andrew AnnexOn Mar 10, 2023, at 11:17 AM, Vincent Sarago ***@***.***> wrote: @AndrewAnnex can you check the PR and see if it works? I'll push some tests soon —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>