Replies: 2 comments 8 replies
-
|
you need to select a e.g 👇 (I needed to encode the URL to make sure it worked) Then I'm not sure the WMS can return the original data, I think they will return PNG tiles 🤷 |
Beta Was this translation helpful? Give feedback.
-
|
PS: I continued a bit playing with DEM algorithms and Titiler, and just tried plugging tiled-COG into maplibre/mapbox as terrain-source (using the Just noticed also the
|
Beta Was this translation helpful? Give feedback.
-
that's really weird. It might be on titiler algorithm side. TBH I haven't used it so the code might be wrong 😓
sure that would be lovely 🙏 |
Beta Was this translation helpful? Give feedback.
-
|
I think this elevation bug actually comes from the mapbox-gl-js and maplibre renderers (see this issue) and how they build the terrain mesh, using height discontinuities, while my geotiff has a nodata value of -9999, masked by alpha. I tried updating the algorithm to take that into account, and set the nodata and alpha=0 elevations to The terrainrgb algorithm could be rewritten to this, which might suffer from less artifacts and floating point precision error. Doing this did solve the retro-conversion in qgis, titiler/src/titiler/core/titiler/core/algorithm/dem.py Lines 227 to 229 in 57b44d6 data_int32 = data.astype(numpy.uint32)
b2 = (data_int32) & 0xff
g2 = (data_int32 >> 8) & 0xff
r2 = (data_int32 >> 16) & 0xffNote to try editing the algorithm and iterate faster, I tried disabling gdal/rasterio caches for the tiles, following this discussion, but it does not seem to work, in case you have a way to do it in mind. |
Beta Was this translation helpful? Give feedback.
-
you don't have to disable GDAL caching but you may want to disable the default titiler/src/titiler/application/titiler/application/main.py Lines 198 to 202 in 57b44d6 TITILER_API_CACHECONTROL="no-store" in the titiler env
While for the terrainrgb algorithm, if you have a patch, we could add an options to make sure it would work with mapbox/maplibre even if it means the return value is not strictly following the terrainrgb equation |
Beta Was this translation helpful? Give feedback.
-
|
Thanks again for that cache control, very useful env_var. I can indeed set the height value returned by the algorithm in the mask area of the masked-array, so that it would be equal to a user-controlled height, eventually behind a url-param flag, but the tiles that fall outside the footprint of the COG will still return empty tiles rather than that user-passed altitude. I don't think I can control with the algorithm the height value of tiles outside the cog footprint - at least at the moment, when loading the tileset in qgis, these appear empty with http message Could we edit how the algorithm applies to tiles endpoint data so that a failover response is returned in case fetched tiles are not correctly formed? Or would you see another way to handle this? At the moment, this is what the terrrain-rgb 3d-visualization looks like. images before fixing EDIT:
|
Beta Was this translation helpful? Give feedback.
-
|
Following the above discussion, I stumbled across an older repo meant to stream and act as a middleware for DEM related data, dem-tiler by @kylebarron - a fork of cogeo-mosaic-tiler which relies on pymartini or pydelatin to encode on-the-fly elevation to quantized-mesh, maybe also under the devseed umbrella. I wanted to see if it would be doable to add this as a dem-algorithm within titiler, with the most prominent difference being that the data processed by the algorithm is not image but quantized-mesh. I've added dependencies (could probably be optional dependencies instead) Quick note the output could also be vector contours based on gdal_contour and tippecanoe for contours MVT creation, the same way dem-tiler does it as well. from pymartini import Martini
import quantized_mesh_encoder
import BytesIO
class QuantizedMesh(BaseAlgorithm):
"""Encode DEM into QuantizedMesh (with pymartini)."""
title: str = "QuantizedMesh"
description: str = "Encode DEM into QuantizedMesh (Cesium Terrain format)."
# parameters
mesh_max_error: float = Field(10.0, ge=1.0, le=100.0)
# metadata
input_nbands: int = 1
# OUTPUT IS ACTUALLY A QUANTIZED-MESH RATHER THAN AN IMAGE
output_nbands: int = 3
output_dtype: str = "uint8"
def __call__(self, img: ImageData) -> ImageData:
"""Encode DEM into QuantizedMesh (with pymartini).
Code from dem-tiler (MIT) https://github.com/kylebarron/dem-tiler/blob/b10d1f61cf89c1a24af5643d0bd4ddee87c52508/dem_tiler/handlers/app.py#L235-L307
"""
# Could also use pydelatin instead of pymartini
mesh_max_error: 10.0
tile_size = 256
tile = img.array[0].astype(numpy.float64)
martini = Martini(tile_size + 1)
mar_tile = martini.create_tile(tile)
vertices, triangles = mar_tile.get_mesh(mesh_max_error) # mesh_max_error)
with BytesIO() as f:
quantized_mesh_encoder.encode(f, vertices, triangles)
f.seek(0)
return ("OK", "application/vnd.quantized-mesh", f.read())
# return ImageData(
# numpy.ma.stack([r, g, b]).astype(self.output_dtype),
# assets=img.assets,
# crs=img.crs,
# bounds=img.bounds,
# ) |
Beta Was this translation helpful? Give feedback.
-
|
🤔 there is a fine limit of what and what not we could do with algorithms. By default the algorithm should return I've played with raster -> vector with https://github.com/cogeotiff/rio-tiler-mvt and also in rio-viz https://github.com/developmentseed/rio-viz/blob/c8676343ae2f84f17322372da481f5194ad3f016/rio_viz/app.py#L567-L588 but I don't see how it could fit within titiler core lib. This could easily be an extension I think |
Beta Was this translation helpful? Give feedback.
-
|
Thanks for referencing the titiler-extensions doc, this seem indeed like the way to go! Also this discussion seem very useful, to help registering a Nice rio-viz and rio-tiler-mvt examples as well. Basing that terrain extension on pymartini could make it less bandwith-intensive than the mvt/pbf points/polygons POC, but these examples are also very useful. It seems that pymartini, quantized-mesh-encoder and pydelatin were compiled against numpy v1, With the below EDIT: It seems that I can reference specific tiles endpoints that do not end with /z/x/y.terrain via another endpoint, Attempt at /layer.json endpoint and /z/x/y.terrain, quantizedMeshExtension based on kylebarron pymartini and quantized-mesh-encoder@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=f"")
# 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))
miny, maxy = (min(ulTile.y, lrTile.y), max(ulTile.y, lrTile.y))
return [ { "startX": minx, "startY": miny, "endX": maxx + 1, "endY": maxy + 1 } ]
# 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'
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:3857",
"bounds": bounds_3857,
"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=f"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.
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."""
tms = factory.supported_tms.get(tileMatrixSetId)
t0 = time.time()
with rasterio.Env(**env):
with factory.reader(
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)
# Could also use pydelatin instead of pymartini
mesh_max_error = 10.0
tile_size = 256
terrain = image.array[0].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], :]))
martini = Martini(tile_size + 1)
mar_tile = martini.create_tile(terrain)
t1 = time.time()
print('duration for martini tile creation', t1 - t01)
vertices, triangles = mar_tile.get_mesh(mesh_max_error) # mesh_max_error)
t2 = time.time()
print('duration for martini mesh get', t2 - t1)
# 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.
rescaled = rescale_positions(
vertices,
terrain,
# bounds=bounds,
# flip_y=True,
# column_row=True
)
t3 = time.time()
print('duration for martini rescaling', t3 - t2)
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=f"application/octet-stream")
# or vnd.quantized-mesh or application/vnd.quantized-mesh;extensions=octvertexnormals terrain
screen |
Beta Was this translation helpful? Give feedback.
-
First of all, sorry for the lengthy read. Original context, the french OSM community is trying to circumvent the (probably temporary) lack of readability of some IGN-produced hillshades, by ingesting new endpoints from the geoplateforme geoservices. A potential contender would eventually be terrain elevation data and upcoming Lidar-HD derived DSM/DEM/DTM products. Long thread here.
After seeing on-the-fly titiler algorithms (really 🤯 feature), I wanted to give them a shot and try to apply custom hillshade to a TMS/XYZ endpoint rather than an existing COG tileset - to have an endpoint that would work on the whole France country rather than having to change the url for each regional assembled COG.
I know, this is a bit of a stretch, I might have been a bit ambitious following another original discussion using GDAL drivers like
GDAL_WMSas titilerurlparam (copy-pasted below) which resulted in the following valid titiler endpoints.It was not an easy task to find true altitude DEMs (most endpoints are composited colored-elevation + hillshade , or Mapzen/Terrarium encoded elevation products like AWS Open Data Terrain tiles), but I finally stumbled on the spanish DEM below. Here is the URL I've come-up with:
The request is currently failing with:
No such band index: 1, and I'm not really sure why, hence opening this discussion. Also, would it make sense to add another input param tohillshadealgorithm so it could ingest not only raw float elevation (or byte8 with z-exaggeration), but also Mapzen/Terrarium encoded altitude? Thanks as always for the impressive titiler!Beta Was this translation helpful? Give feedback.
All reactions