Skip to content

Commit 790fc97

Browse files
Fix TerrainRGB algorithm and param user-controlled nodata-height (#1116)
* Fix TerrainRGB algorithm and param user-controlled nodata-height Added two params `use_nodata_height and nodata_height` (or could have used a field which cold be either undefined or that height value? * Make use of user-controlled height in terrainrgb * Add user-controlled nodata-height for terrarium as well * Add z_exaggeration parameter to hillshade/slope algorithms defaults to 1, applied to gradients directly * Add slope in algorithms doc * Make terrainRGB/terrarium docs links instead of plain urls * Add parameters hints to hillshade/contours * nodata_height optional Co-authored-by: Vincent Sarago <vincent.sarago@gmail.com> * check nodata_height not None Co-authored-by: Vincent Sarago <vincent.sarago@gmail.com> * Update src/titiler/core/titiler/core/algorithm/dem.py Co-authored-by: Vincent Sarago <vincent.sarago@gmail.com> * Update src/titiler/core/titiler/core/algorithm/dem.py Co-authored-by: Vincent Sarago <vincent.sarago@gmail.com> * from typing import Optional * WIP /cog/viewer algorithms inputs Adding algorithm params dynamically based on the /algorithms endpoint, stored in scope, and params updated on change of selected algorithm * Final touch-ups to /cog/viewer for algorithms visualization Uses number inputs if param is integer or number, otherwise text input (eg nodata-height which can be either null or number) Updates tilejson url when algorithm or its params are changed * Revert "Final touch-ups to /cog/viewer for algorithms visualization" This reverts commit bf094ae. * Revert "WIP /cog/viewer algorithms inputs" This reverts commit 31d5277. * Add nodata_height coverage to tests * Discard solution 1 for testing terrainrgb and terrarium * pre-commit run --all-files --------- Co-authored-by: Vincent Sarago <vincent.sarago@gmail.com>
1 parent df1233c commit 790fc97

File tree

3 files changed

+50
-7
lines changed

3 files changed

+50
-7
lines changed

docs/src/advanced/Algorithms.md

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,11 @@ The algorithms are meant to overcome the limitation of `expression` (using [nume
66

77
We added a set of custom algorithms:
88

9-
- `hillshade`: Create hillshade from elevation dataset
10-
- `contours`: Create contours lines (raster) from elevation dataset
11-
- `terrarium`: Mapzen's format to encode elevation value in RGB values (https://github.com/tilezen/joerd/blob/master/docs/formats.md#terrarium)
12-
- `terrainrgb`: Mapbox's format to encode elevation value in RGB values (https://docs.mapbox.com/data/tilesets/guides/access-elevation-data/)
9+
- `hillshade`: Create hillshade from elevation dataset (parameters: azimuth (45), angle_altitude(45))
10+
- `contours`: Create contours lines (raster) from elevation dataset (parameters: increment (35), thickness (1))
11+
- `slope`: Create degrees of slope from elevation dataset
12+
- `terrarium`: [Mapzen's format]((https://github.com/tilezen/joerd/blob/master/docs/formats.md#terrarium)) to encode elevation value in RGB values `elevation = (red * 256 + green + blue / 256) - 32768`
13+
- `terrainrgb`: [Mapbox](https://docs.mapbox.com/data/tilesets/guides/access-elevation-data/)/[Maptiler](https://docs.maptiler.com/guides/map-tilling-hosting/data-hosting/rgb-terrain-by-maptiler/)'s format to encode elevation value in RGB values `elevation = -10000 + ((red * 256 * 256 + green * 256 + blue) * 0.1)`
1314
- `normalizedIndex`: Normalized Difference Index (e.g NDVI)
1415
- `cast`: Cast data to integer
1516
- `floor`: Round data to the smallest integer

src/titiler/core/tests/test_algorithms.py

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -236,6 +236,16 @@ def test_terrarium():
236236
assert out.array.dtype == "uint8"
237237
assert out.array[0, 0, 0] is numpy.ma.masked
238238

239+
# works on the above masked array img, with algo which was passed nodata_height
240+
nodata_height = 10.0
241+
algo = default_algorithms.get("terrarium")(nodata_height=nodata_height)
242+
out = algo(img)
243+
masked = out.array[:, arr.mask[0, :, :]]
244+
masked_height = (masked[0] * 256 + masked[1] + masked[2] / 256) - 32768
245+
numpy.testing.assert_array_equal(
246+
masked_height, nodata_height * numpy.ones((100 * 100), dtype="bool")
247+
)
248+
239249

240250
def test_terrainrgb():
241251
"""test terrainrgb."""
@@ -259,6 +269,18 @@ def test_terrainrgb():
259269
assert out.array.dtype == "uint8"
260270
assert out.array[0, 0, 0] is numpy.ma.masked
261271

272+
# works on the above masked array img, with algo which was passed nodata_height
273+
nodata_height = 10.0
274+
algo = default_algorithms.get("terrainrgb")(nodata_height=nodata_height)
275+
out = algo(img)
276+
masked = out.array[:, arr.mask[0, :, :]]
277+
masked_height = -10000 + (
278+
((masked[0] * 256 * 256) + (masked[1] * 256) + masked[2]) * 0.1
279+
)
280+
numpy.testing.assert_array_equal(
281+
masked_height, nodata_height * numpy.ones((100 * 100), dtype="bool")
282+
)
283+
262284

263285
def test_ops():
264286
"""test ops: cast, ceil and floor."""

src/titiler/core/titiler/core/algorithm/dem.py

Lines changed: 23 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
11
"""titiler.core.algorithm DEM."""
22

3+
from typing import Optional
4+
35
import numpy
46
from pydantic import Field
57
from rasterio import windows
@@ -22,6 +24,7 @@ class HillShade(BaseAlgorithm):
2224
azimuth: int = Field(45, ge=0, le=360)
2325
angle_altitude: float = Field(45.0, ge=-90.0, le=90.0)
2426
buffer: int = Field(3, ge=0, le=99)
27+
z_exaggeration: float = Field(1.0, ge=1e-6, le=1e6)
2528

2629
# metadata
2730
input_nbands: int = 1
@@ -31,6 +34,8 @@ class HillShade(BaseAlgorithm):
3134
def __call__(self, img: ImageData) -> ImageData:
3235
"""Create hillshade from DEM dataset."""
3336
x, y = numpy.gradient(img.array[0])
37+
x *= self.z_exaggeration
38+
y *= self.z_exaggeration
3439
slope = numpy.pi / 2.0 - numpy.arctan(numpy.sqrt(x * x + y * y))
3540
aspect = numpy.arctan2(-x, y)
3641
azimuth = 360.0 - self.azimuth
@@ -71,6 +76,7 @@ class Slope(BaseAlgorithm):
7176

7277
# parameters
7378
buffer: int = Field(3, ge=0, le=99, description="Buffer size for edge effects")
79+
z_exaggeration: float = Field(1.0, ge=1e-6, le=1e6)
7480

7581
# metadata
7682
input_nbands: int = 1
@@ -86,6 +92,8 @@ def __call__(self, img: ImageData) -> ImageData:
8692
pixel_size_y = abs(img.transform[4])
8793

8894
x, y = numpy.gradient(img.array[0])
95+
x *= self.z_exaggeration
96+
y *= self.z_exaggeration
8997
dx = x / pixel_size_x
9098
dy = y / pixel_size_y
9199

@@ -161,6 +169,7 @@ class Terrarium(BaseAlgorithm):
161169

162170
title: str = "Terrarium"
163171
description: str = "Encode DEM into RGB (Mapzen Terrarium)."
172+
nodata_height: Optional[float] = Field(None, ge=-99999.0, le=99999.0)
164173

165174
# metadata
166175
input_nbands: int = 1
@@ -170,6 +179,10 @@ class Terrarium(BaseAlgorithm):
170179
def __call__(self, img: ImageData) -> ImageData:
171180
"""Encode DEM into RGB."""
172181
data = numpy.clip(img.array[0] + 32768.0, 0.0, 65535.0)
182+
if self.nodata_height is not None:
183+
data[img.array.mask[0]] = numpy.clip(
184+
self.nodata_height + 32768.0, 0.0, 65535.0
185+
)
173186
r = data / 256
174187
g = data % 256
175188
b = (data * 256) % 256
@@ -191,6 +204,7 @@ class TerrainRGB(BaseAlgorithm):
191204
# parameters
192205
interval: float = Field(0.1, ge=0.0, le=1.0)
193206
baseval: float = Field(-10000.0, ge=-99999.0, le=99999.0)
207+
nodata_height: Optional[float] = Field(None, ge=-99999.0, le=99999.0)
194208

195209
# metadata
196210
input_nbands: int = 1
@@ -224,9 +238,15 @@ def _range_check(datarange):
224238
if _range_check(datarange):
225239
raise ValueError(f"Data of {datarange} larger than 256 ** 3")
226240

227-
r = ((((data // 256) // 256) / 256) - (((data // 256) // 256) // 256)) * 256
228-
g = (((data // 256) / 256) - ((data // 256) // 256)) * 256
229-
b = ((data / 256) - (data // 256)) * 256
241+
if self.nodata_height is not None:
242+
data[img.array.mask[0]] = (
243+
self.nodata_height - self.baseval
244+
) / self.interval
245+
246+
data_int32 = data.astype(numpy.int32)
247+
b = (data_int32) & 0xFF
248+
g = (data_int32 >> 8) & 0xFF
249+
r = (data_int32 >> 16) & 0xFF
230250

231251
return ImageData(
232252
numpy.ma.stack([r, g, b]).astype(self.output_dtype),

0 commit comments

Comments
 (0)