Fix TerrainRGB algorithm and param user-controlled nodata-height

[jo-chemla]

• Fixes Terrain-RGB algorithm
• Added two params use_nodata_height and nodata_height (or could have used a field which could be either undefined or that height value?).
• Also added user-controlled nodata-height parametrization for terrarium

Can be tested via

'tiles': ['http://localhost:8000/cog/tiles/WebMercatorQuad/{z}/{x}/{y}@1x.png?algorithm=terrainrgb&algorithm_params={"use_nodata_height":true,"nodata_height":0.0}&nodata=-9999&url=https://3d.iconem.com/greece/amphipolis/Amphipolis_general_dsm_10cm_cog_cropped_alpha.tif']

TODO: Also should take a moment to add suggested updates from discussion #1110

• terrarium/terrain-rgb intepretation for dem-related algorithms like hillshade & slope
• z-exaggeration algorithm-param
• algorithms doc update
• Just noticed also the /cog/viewer endpoint to help users play with url-params, bands expressions etc. Extend single-band user-parametrization to also allow for algorithms and a few other url-params like rescale etc? other PR Add algorithms selection to /cog/viewer and custom expression support #1119

[vincentsarago]

we will need to update the tests in

titiler/src/titiler/core/tests/test_algorithms.py

Lines 69 to 107 in 57b44d6

	def test_terrain_algo():
	"""test terrain algorithm deps."""
	# Add the `Multiply` algorithm to the default ones
	app = FastAPI()
	arr = numpy.random.randint(0, 3000, (1, 256, 256))
	@app.get("/", response_class=Response)
	def main(algorithm=Depends(default_algorithms.dependency)):
	"""endpoint."""
	img = ImageData(arr)
	if algorithm:
	img = algorithm(img)
	return Response(img.render(img_format="PNG"), media_type="image/png")
	client = TestClient(app)
	# MAPBOX Terrain RGB
	response = client.get("/", params={"algorithm": "terrainrgb"})
	assert response.status_code == 200
	with MemoryFile(response.content) as mem:
	with mem.open() as dst:
	data = dst.read().astype(numpy.float64)
	# https://docs.mapbox.com/data/tilesets/guides/access-elevation-data/
	elevation = -10000 + (((data[0] * 256 * 256) + (data[1] * 256) + data[2]) * 0.1)
	numpy.testing.assert_array_equal(elevation, arr[0])
	# TILEZEN Terrarium
	response = client.get("/", params={"algorithm": "terrarium"})
	assert response.status_code == 200
	with MemoryFile(response.content) as mem:
	with mem.open() as dst:
	data = dst.read().astype(numpy.float64)
	# https://github.com/tilezen/joerd/blob/master/docs/formats.md#terrarium
	elevation = (data[0] * 256 + data[1] + data[2] / 256) - 32768
	numpy.testing.assert_array_equal(elevation, arr[0])

to account for the new options

[vincentsarago]

I'm usually testing with https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif but I cannot get it to work 🤔

[jo-chemla]

Weird, will need to check. I've also updated the /cog/viewer endpoint to also support algorithm select and parameters edition, should be good to go!

I'll try to have a look at your feedback tomorrow!

[vincentsarago]

I think we should decouple the viewer part from this PR.

e.g algorithm might need more than one band, so I don't think they should be in the 1b viz

[vincentsarago]

FYI I've started a PR in rio-viz to support raster-dem viz developmentseed/rio-viz#64 ;-)

[jo-chemla]

Wow that was fast, 3d maplibre support looking nice! I can revert the changes specific to the viewer in that titiler PR, should I port the algorithms selection to another PR in Rio-viz instead? Still a wip, it can indeed be useful to other bands as well, will need to strengthen that wip.

[vincentsarago]

should I port the algorithms selection to another PR in Rio-viz instead? Still a wip, it can indeed be useful to other bands as well, will need to strengthen that wip.

@jo-chemla I've personally don't have a lot of time to work on this (I've already spend too much time 🙈) so I would say you can but don't expect too much 😅

[jo-chemla]

Alright, thanks again for getting back! I've reverted commits dedicated to /cog/viewer.

Trying to gather what's missing before merging:

• testing with https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif and understand the issue
• update tests in

  titiler/src/titiler/core/tests/test_algorithms.py

  Lines 69 to 107 in 57b44d6

  	def test_terrain_algo():
  	"""test terrain algorithm deps."""
  	# Add the `Multiply` algorithm to the default ones
  	app = FastAPI()
  	arr = numpy.random.randint(0, 3000, (1, 256, 256))
  	@app.get("/", response_class=Response)
  	def main(algorithm=Depends(default_algorithms.dependency)):
  	"""endpoint."""
  	img = ImageData(arr)
  	if algorithm:
  	img = algorithm(img)
  	return Response(img.render(img_format="PNG"), media_type="image/png")
  	client = TestClient(app)
  	# MAPBOX Terrain RGB
  	response = client.get("/", params={"algorithm": "terrainrgb"})
  	assert response.status_code == 200
  	with MemoryFile(response.content) as mem:
  	with mem.open() as dst:
  	data = dst.read().astype(numpy.float64)
  	# https://docs.mapbox.com/data/tilesets/guides/access-elevation-data/
  	elevation = -10000 + (((data[0] * 256 * 256) + (data[1] * 256) + data[2]) * 0.1)
  	numpy.testing.assert_array_equal(elevation, arr[0])
  	# TILEZEN Terrarium
  	response = client.get("/", params={"algorithm": "terrarium"})
  	assert response.status_code == 200
  	with MemoryFile(response.content) as mem:
  	with mem.open() as dst:
  	data = dst.read().astype(numpy.float64)
  	# https://github.com/tilezen/joerd/blob/master/docs/formats.md#terrarium
  	elevation = (data[0] * 256 + data[1] + data[2] / 256) - 32768
  	numpy.testing.assert_array_equal(elevation, arr[0])

  to account for new options
• pass the tests

I think I can PR the /cog/viewer in a new PR almost as-is, and

• activate that algorithm-section in only update-1b first, but filter only algos that apply to single-band data
• then I can also push the algorithms-section to multi-bands, but there might be some duplicated code
• after your comment on terrain in rio-viz, I thought rio-viz was the module/template behind the cog_viewer.html endpoint, but after inspection they seem unrelated - or at least code and changes made in one have to be ported to the other, right? In that case I might focus on the titiler /cog/viewer only, and in the future we might port the changes to rio-viz, is that correct?

[vincentsarago]

🙏 @jo-chemla

after your comment on terrain in rio-viz, I thought rio-viz was the module/template behind the cog_viewer.html endpoint, but after inspection they seem unrelated - or at least code and changes made in one have to be ported to the other, right? In that case I might focus on the titiler /cog/viewer only, and in the future we might port the changes to rio-viz, is that correct?

Yeah they are different templates (because at the beginning they were quite different but with time they evolved to be a bit too close now 😓

let's not worry about rio-viz for now, you can focus on titiler cog viewer 🙏

[jo-chemla]

Just had a quick look,

• Seems that the error with the swisstopo dem was just a trailing space, it works for me with the /cog/preview & tiles endpoints, with or without algorithm specified as url-param
• current tests do pass with python -m pytest src/titiler/core --cov=titiler.core --cov-report=xml --cov-append --cov-report=term-missing, with a coverage that can be improved here

Name                                                   Stmts   Miss Branch BrPart  Cover   Missing
src\titiler\core\titiler\core\algorithm\dem.py           123      3     10      5    94%   50->61, 102->113, 182, 236, 239

The missing tests are:

• on pre-existing missing-branches if self.buffer:
• and newly, non-covered branches are the 2 added tests if self.nodata_height is not None:. I've added a test where nodata_height url-param is passed, with a mask, and check that the values of the output in the mask area equate the provided height.

The thing is I'm not sure which test is best among the 2 solutions below. Either edit test_terrain_algo with a client.get response where I add the url-param, or rather edit the test_terrarium methods passing the option to the python method. Just committed, and test do pass, you tell me.

# SOLUTION 1
def test_terrain_algo():
    #...
    #def main(algorithm=Depends(default_algorithms.dependency)):
    # client = TestClient(app)
    #   ...

    # SOLUTION 1
    # test nodata_height 
    # MAPBOX Terrain RGB
    response = client.get("/", params={"algorithm": "terrainrgb", "algorithm_params":json.dumps({"nodata_height":10.0})})
    assert response.status_code == 200
    with MemoryFile(response.content) as mem:
        with mem.open() as dst:
            data = dst.read().astype(numpy.float64)
    elevation = -10000 + (((data[0] * 256 * 256) + (data[1] * 256) + data[2]) * 0.1)
    numpy.testing.assert_array_equal(elevation, arr[0])

    # TILEZEN Terrarium
    response = client.get("/", params={"algorithm": "terrarium", "algorithm_params":json.dumps({"nodata_height":10.0})})
    assert response.status_code == 200
    with MemoryFile(response.content) as mem:
        with mem.open() as dst:
            data = dst.read().astype(numpy.float64)
    elevation = (data[0] * 256 + data[1] + data[2] / 256) - 32768
    numpy.testing.assert_array_equal(elevation, arr[0])

# SOLUTION 2
def test_terrarium():
    # ...

    # nodata_height
    nodata_height = 10.0
    algo = default_algorithms.get("terrarium")(nodata_height=nodata_height)
    arr = numpy.ma.MaskedArray(
        numpy.random.randint(0, 5000, (1, 256, 256), dtype="uint16"),
        mask=numpy.zeros((1, 256, 256), dtype="bool"),
    )
    arr.mask[0, 0:100, 0:100] = True
    img = ImageData(arr)
    out = algo(img)
    masked = out.array[:, arr.mask[0,:,:]]
    masked_height = (masked[0] * 256 + masked[1] + masked[2] / 256) - 32768
    print(masked_height.shape)
    numpy.testing.assert_array_equal(masked_height, nodata_height * numpy.ones((100 * 100), dtype="bool"))

[vincentsarago]

@jo-chemla Solution 2 is fine (which is how the tests are written for now) 🙏

Yeah I think both solution are fine. I'm ok with what we have right now 👍

we don't need to worry to much about the coverage %

[jo-chemla]

Thanks, from my point of view we should be ready to go then!
In parallel I've finalized the touchups for the PR dedicated to /cog/viewer, which now works for both 1b and 3b afaict, please do tell me if things are not done correctly on the other thread #1119.

Also, the last item from the original discussion was to generalize algorithms to ingest not only 1b DEMs, but also allow transformation of the input data from terrarium/terrainrgb into altitude/height data (or algorithms chaining) so it could be used on-the-fly by algorithms (hillshade etc). I'm not sure how this could/should work and it can probably be the subject of another PR if deemed useful.

Thanks for all the help as always!

[vincentsarago]

@jo-chemla could you run pre-commit on your branch?

python -m pip install pre-commit 
# in titiler root 
pre-commit install 
pre-commit run --all-files

Also, the last item from the original discussion was to generalize algorithms to ingest not only 1b DEMs, but also allow transformation of the input data from terrarium/terrainrgb into altitude/height data (or algorithms chaining) so it could be used on-the-fly by algorithms (hillshade etc). I'm not sure how this could/should work and it can probably be the subject of another PR if deemed useful.

Do you mean go from RGBA (terranium) to Float? (I think you could use band expression for this)

algorithms chaining

we could allow algorithm to be a list but then it might be a bit messy when wanting to pass options to each algorithm.

FYI We've recently worked on OpenEO https://github.com/sentinel-hub/titiler-openeo and the custom processing language is quite interesting. We might bring this to titiler (if any customer need it 🙈 )

[jo-chemla]

Thanks again. Just ran the pre-commit and committed, so should be good to go on this one!

Nice to see the custom Processing language for OpenEO, first time I hear about it 👀. It does look a lot like Node editors for 3D-artist related workflows, eg in case any inspiration could be useful:

• chaining Blender nodes for material or geometry,
• threejs shading language TSL NodeMaterial,
• or Alicevision meshroom processing
• PDAL pipelines is another interesting take for more linear workflows.

And for algorithms chaining, 🤯 indeed the band expression is the way to prepare the data before algorithm execution (eg the existing aws open-data terrain tiles which are stored as terrarium with url https://s3.amazonaws.com/elevation-tiles-prod/terrarium/${z}/${x}/${y}.png can be wrapped into a GDAL_WMS TMS service passed as url, with expression=(b1 * 256 + b2 + b3 / 256) - 32768, and algorithm eg hillshade. 🙈 Full url here

Makes me want to add expression support before merging #1119

[vincentsarago]

Again I don't want to add too much stuff into the viewer (mostly because it's written in pure htm/js I did and I'm not a frontend dev 😅)

If you really want to play/visualize the data I often use the /map endpoint directly

https://titiler.xyz/cog/WebMercatorQuad/map?url=%3CGDAL_WMS%3E%3CService%20name%3D%27TMS%27%3E%3CServerUrl%3Ehttps%3A%2F%2Fs3.amazonaws.com%2Felevation-tiles-prod%2Fterrarium%2F%24%7Bz%7D%2F%24%7Bx%7D%2F%24%7By%7D.png%3C%2FServerUrl%3E%3C%2FService%3E%3CDataWindow%3E%3CUpperLeftX%3E-20037508.34%3C%2FUpperLeftX%3E%3CUpperLeftY%3E20037508.34%3C%2FUpperLeftY%3E%3CLowerRightX%3E20037508.34%3C%2FLowerRightX%3E%3CLowerRightY%3E-20037508.34%3C%2FLowerRightY%3E%3CTileLevel%3E18%3C%2FTileLevel%3E%3CTileCountX%3E1%3C%2FTileCountX%3E%3CTileCountY%3E1%3C%2FTileCountY%3E%3CYOrigin%3Etop%3C%2FYOrigin%3E%3C%2FDataWindow%3E%3CProjection%3EEPSG%3A3857%3C%2FProjection%3E%3CBlockSizeX%3E256%3C%2FBlockSizeX%3E%3CBlockSizeY%3E256%3C%2FBlockSizeY%3E%3CBandsCount%3E3%3C%2FBandsCount%3E%3CCache%20%2F%3E%3C%2FGDAL_WMS%3E&bidx=1&expression=%28b1%20%2A%20256%20%2B%20b2%20%2B%20b3%20%2F%20256%29%20-%2032768&max_size=1024&algorithm=hillshade

[jo-chemla]

Alright, just tried adding the expression support to /cog/viewer, only to 1b section since I could only test it on this. I'll let you decide whether or not you want that addition merged!