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.
1.
For context, see discussion developmentseed/titiler#1110
Side-note: while my TiTiler PR http://github.com/developmentseed/titiler/pull/1119 to implement
hillshade, algorithms and custom expressionsis not merged, indices can be computed on the fly via apps described in this thread TerriaJS/terriajs#7665 (custom raster symbology requested on terriamap)2.
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
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