Titiler implement dem algos for terrain-rgb, terrarium, quantized-mesh

[jo-chemla]

See the discussion here developmentseed/titiler#1110 (reply in thread) which mentions some wip almost working

See a quantized-mesh discussion revived by Garrett on specification of the JSON structure and CRS of the quantized-mesh spec here CesiumGS/quantized-mesh#15 (comment)

also branch algo-pymartini (could not commit because ruff or dem-algorithms

C:\Dev\Iconem\titiler\src\titiler\application\titiler\application\main.py

import time
from io import BytesIO
from typing import Dict, Literal, Tuple

import mercantile
import morecantile
import numpy
import quantized_mesh_encoder
from attrs import define
from fastapi import Depends, Path

# http://127.0.0.1:8000/cog/tiles_terrain/WebMercatorQuad/16/32975/22392.terrain?url=http%3A%2F%2F127.0.0.1%3A8081%2Fimage_6.cog.tif
from pymartini import Martini, decode_ele
from pymartini import rescale_positions as martini_rescale_positions
from pyproj import CRS, Transformer
from rio_cogeo.cogeo import cog_info
from rio_cogeo.models import Info
from starlette.responses import Response
from typing_extensions import Annotated

from titiler.core import factory
from titiler.core.resources.responses import JSONResponse


@define(kw_only=True)
class quantizedMeshExtension(factory.FactoryExtension):
    def register(self, factory: factory.TilerFactory):
        """Register endpoint to the tiler factory."""

        #     """
        @factory.router.get("/test", response_model=Info, response_class=JSONResponse)
        def test():
            return Response("test", media_type="")

        # Test with
        # http://127.0.0.1:8000/cog/tiles_terrain/WebMercatorQuad/layer.json?url=https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif
        # http://127.0.0.1:8000/cog/tiles_terrain/WebMercatorQuad/layer.json?url=WMS%3Ahttps%3A%2F%2Fdata.geopf.fr%2Fwms-r%3FSERVICE%3DWMS%26VERSION%3D1.3.0%26REQUEST%3DGetMap%26LAYERS%3DIGNF_LIDAR-HD_MNS_ELEVATION.ELEVATIONGRIDCOVERAGE.WGS84G%26STYLES%3Dnormal%26FORMAT%3Dimage%252Fgeotiff%26BBoxOrder%3DyxYX
        # http://127.0.0.1:8000/cog/tiles_terrain/WebMercatorQuad/layer.json?url=http://127.0.0.1:8081/image_6.cog.tif

        @factory.router.get(
            "/tiles_terrain/{tileMatrixSetId}/layer.json", response_class=JSONResponse
        )
        def tiles_json(
            tileMatrixSetId: Literal[tuple(factory.supported_tms.list())] = Path(  # type: ignore
                description="TileMatrixSet Name",
            ),
            # Allowed values for the projection are EPSG:3857 (Web Mercator as used by Google Maps) and EPSG:4326 (Lat/Lng coordinates in the Global-Geodetic System). It is worth noting that the EPSG3857 projection has only 1 tile at the root (zoom level 0) while EPSG:4326 has 2.
            src_path=Depends(factory.path_dependency),
        ):
            # https://github.com/developmentseed/supermorecado/blob/7b2a68fce1c25921db81dd317ee4eb5ca574171e/supermorecado/burntiles.py#L135
            WEB_MERCATOR_TMS = morecantile.tms.get("WebMercatorQuad")
            tms = WEB_MERCATOR_TMS

            def tile_extrema(
                bounds_4326: Tuple[float, float, float, float], zoom: int
            ) -> Dict:
                """Tiles min/max at the given zoom for bounds."""
                ulTile = tms.tile(bounds_4326[0], bounds_4326[3], zoom)
                lrTile = tms.tile(bounds_4326[2], bounds_4326[1], zoom)
                minx, maxx = (min(ulTile.x, lrTile.x), max(ulTile.x, lrTile.x, 2**zoom))
                miny, maxy = (min(ulTile.y, lrTile.y), max(ulTile.y, lrTile.y, 2**zoom))
                # return [ { "startX": minx, "startY": miny, "endX": maxx, "endY": maxy } ] # 3857
                return [
                    {"startX": minx * 2, "startY": miny, "endX": maxx * 2, "endY": maxy}
                ]  # 4326
                # return {
                #     "x": {"min": minx, "max": maxx + 1},
                #     "y": {"min": miny, "max": maxy + 1},
                # }

            def transform_boundingBox(
                boundingBox: Tuple[float, float, float, float],
                src_crs: CRS,
                dst_crs: CRS,
            ) -> Tuple[float, float, float, float]:
                """Transform bounds to a different CRS."""
                transformer = Transformer.from_crs(src_crs, dst_crs, always_xy=True)
                ll_ur_dst_crs = list(
                    transformer.itransform([boundingBox[0:2], boundingBox[2:4]])
                )
                bounds_dst_crs = [
                    ll_ur_dst_crs[0][0],
                    ll_ur_dst_crs[0][1],
                    ll_ur_dst_crs[1][0],
                    ll_ur_dst_crs[1][1],
                ]
                return bounds_dst_crs

            coginfo = cog_info(src_path)
            geo = coginfo.GEO
            boundingBox = geo.BoundingBox
            # project bbox to TMS coords 3857, not generic assumption is output CRS is WebMercatorQuad
            # ll = tms.xy(*boundingBox[0:2])
            # ur = tms.xy(*boundingBox[2:4])
            # crs = CRS.from_epsg(3857)
            crs = geo.CRS
            bounds_3857 = transform_boundingBox(boundingBox, crs, "EPSG:3857")
            bounds_4326 = transform_boundingBox(boundingBox, crs, "EPSG:4326")

            # bounds = [ll.x, ll.y, ur.x, ur.y]
            zoom_min = (
                0 * geo.MinZoom
            )  # has to be zero for all available tiles array list to fill properly
            zoom_max = geo.MaxZoom
            print(coginfo)
            # src_path_2 = 'http://127.0.0.1:8081/image_6.cog.tif'
            # See more on specification: https://github.com/CesiumGS/quantized-mesh/issues/15
            # and interesting post here https://bertt.wordpress.com/2018/11/26/visualizing-terrains-with-cesium-ii/
            tilejson = {
                "tilejson": "2.1.0",
                "name": "titiler quantized mesh endpoint",
                "description": "",
                "version": "1.1.0",
                "format": "quantized-mesh-1.0",
                "attribution": "",
                "schema": "tms",
                "tiles": [f"./{{z}}/{{x}}/{{y}}.terrain?url={src_path}"],
                # "projection": "EPSG:4326",
                "projection": "EPSG:4326",
                # always 4326 https://help.agi.com/TerrainServer/RESTAPIGuide.html despite PR https://github.com/CesiumGS/quantized-mesh/issues/10 https://github.com/CesiumGS/cesium/pull/8563
                # "bounds": bounds_4326,
                "bounds": bounds_4326,
                # "bounds": [ -180.00, -90.00, 180.00, 90.00 ],
                "available": list(
                    map(
                        lambda z: tile_extrema(bounds_4326, z),
                        range(zoom_min, zoom_max),
                    )
                ),
                # "bounds": [ 0.00, -90.00, 180.00, 90.00 ],
                # "available": [
                #     [ { "startX": 0, "startY": 0, "endX": 0, "endY": 0 } ],
                #     ...
                #     [ { "startX": 0, "startY": 0, "endX": 269011, "endY": 416820 } ]
                #     ]
            }
            return JSONResponse(tilejson, media_type="application/json")

        @factory.router.get(
            "/tiles_terrain/{tileMatrixSetId}/{z}/{x}/{y}.terrain",
            response_class=Response,
        )
        def tile(
            z: Annotated[
                int,
                Path(
                    description="Identifier (Z) selecting one of the scales defined in the TileMatrixSet and representing the scaleDenominator the tile.",
                ),
            ],
            x: Annotated[
                int,
                Path(
                    description="Column (X) index of the tile on the selected TileMatrix. It cannot exceed the MatrixHeight-1 for the selected TileMatrix.",
                ),
            ],
            y: Annotated[
                int,
                Path(
                    description="Row (Y) index of the tile on the selected TileMatrix. It cannot exceed the MatrixWidth-1 for the selected TileMatrix.",
                ),
            ],
            tileMatrixSetId: Literal[tuple(factory.supported_tms.list())] = Path(  # type: ignore
                description="TileMatrixSet Name",
            ),
            # Allowed values for the projection are EPSG:3857 (Web Mercator as used by Google Maps) and EPSG:4326 (Lat/Lng coordinates in the Global-Geodetic System). It is worth noting that the EPSG3857 projection has only 1 tile at the root (zoom level 0) while EPSG:4326 has 2.
            input_format: Annotated[
                Literal[tuple(["raw_altitude", "terrarium", "terrainrgb"])],
                Path(  # type: ignore
                    description="input type, either default raw altitude, or terrarium or terrainrgb",
                ),
            ] = "terrarium",
            src_path=Depends(factory.path_dependency),
            reader_params=Depends(factory.reader_dependency),
            tile_params=Depends(factory.tile_dependency),
            layer_params=Depends(factory.layer_dependency),
            dataset_params=Depends(factory.dataset_dependency),
            env=Depends(factory.environment_dependency),
        ):
            """Create map tile from a dataset."""
            # Could also use pydelatin instead of pymartini
            mesh_max_error = 10.0
            tile_size = 256
            tms = factory.supported_tms.get(tileMatrixSetId)
            t0 = time.time()
            input_format = (
                input_format if input_format is not None else "raw_altitude",
            )
            if False:
                with rasterio.Env(**env):
                    with (
                        factory.reader(  # self.reader where self extends BaseFactory in https://github.com/developmentseed/titiler/blob/0f3e3c78909f05b2bfbef21438466c3133c508fa/src/titiler/core/titiler/core/factory.py#L220
                            src_path, tms=tms, **reader_params.as_dict()
                        ) as src_dst
                    ):
                        print("src_dst", src_dst)
                        print(
                            "BEFORE src_dst.tile: ",
                            x,
                            y,
                            z,
                            **tile_params.as_dict(),
                            **layer_params.as_dict(),
                            **dataset_params.as_dict(),
                        )
                        # TODO THIS TAKES FOREVER at 0,0,0
                        image = src_dst.tile(
                            x,
                            y,
                            z,
                            tilesize=256,
                            **tile_params.as_dict(),
                            **layer_params.as_dict(),
                            **dataset_params.as_dict(),
                        )

                        print("image", image)

                t01 = time.time()
                print("duration for src_dst.tile call", t01 - t0)

                terrain = image.array[0].astype(numpy.float32)
                terrain[image.mask] = 0

                if input_format == "raw_altitude":
                    print("raw_altitude")
                elif input_format == "terrarium":
                    print("terrarium")
                    terrain = decode_ele(terrain, "terrarium")
                elif input_format == "terrainrgb":
                    print("terrainrgb")
                    terrain = decode_ele(terrain, "terrainrgb")
            else:
                # terrain = numpy.ones((tile_size, tile_size)) * 10.0
                # terrain = numpy.random.rand(tile_size, tile_size) * 100000.0
                terrain = numpy.mgrid[0:tile_size, 0:tile_size][1] / 256 * 100000
                terrain = terrain.astype(numpy.float32)
            # repeat last row and col - could be first as well
            terrain = numpy.hstack((terrain, terrain[:, [-1]]))
            terrain = numpy.vstack((terrain, terrain[[-1], :]))

            #
            flip_y = True
            tile_bounds_3857 = tms.xy_bounds(morecantile.Tile(x, y, z))
            morecantile_tms = morecantile.tms.get("WGS1984Quad")  # or 'WorldCRS84Quad'
            tile_bounds_4326 = morecantile_tms.xy_bounds(morecantile.Tile(x, y, z))
            print(z, x, y, tile_bounds_4326)
            # tile_bounds_4326  =[tile_bounds_4326[0] * 2, tile_bounds_4326[1], tile_bounds_4326[2] * 2, tile_bounds_4326[3]]
            # print(tile_bounds_4326)
            use_delatin = False
            bounds = mercantile.bounds(mercantile.Tile(x, y, z))

            mesh_max_error = float(mesh_max_error)
            if use_delatin:
                print("not supported yet")
                # tin = Delatin(terrain, max_error=mesh_max_error)
                # vertices, triangles = tin.vertices, tin.triangles.flatten()
                # rescaled = delatin_rescale_positions(vertices, bounds, flip_y=flip_y)
            else:
                martini = Martini(tile_size + 1)
                mar_tile = martini.create_tile(terrain)
                t1 = time.time()
                print("duration for martini tile creation", t1 - t0)
                vertices, triangles = mar_tile.get_mesh(
                    mesh_max_error
                )  # mesh_max_error)
                t2 = time.time()
                print("duration for martini mesh get", t2 - t1)
                rescaled = martini_rescale_positions(
                    vertices, terrain, bounds=tile_bounds_4326, flip_y=flip_y
                )
                print(rescaled)
                t3 = time.time()
                print("duration for martini rescaling", t3 - t2)

            print(
                vertices.shape,
                triangles.shape,
                rescaled.shape,
                vertices,
                triangles,
                rescaled,
            )
            with BytesIO() as f:
                quantized_mesh_encoder.encode(f, rescaled, triangles)
                t4 = time.time()
                print("duration for quantized_mesh_encoder encoding", t4 - t3)
                print("duration total", t4 - t0)
                f.seek(0)
                # return ("OK", "application/vnd.quantized-mesh", f.read())
                return Response(f.read(), media_type="application/octet-stream")
                # or vnd.quantized-mesh or application/vnd.quantized-mesh;extensions=octvertexnormals terrain

            # Each is a flat numpy array. Vertices represents the interleaved 2D
            # coordinates of each vertex, e.g. [x0, y0, x1, y1, ...]. If you need 3D
            # coordinates, you can use the rescale_positions helper function described
            # below.
            # triangles represents indices within the vertices array. So [0, 1, 3, ...] would use the first, second, and fourth vertices within the vertices array as a single triangle.

details

Eventually polish the dem-algo with a demo showcasing public terrain instances, ingesting terrain-rgb/terrarium/quantized-mesh/3dtiles for somme terrain sources listed in iconem wiki here