Fix TerrainRGB algorithm and param user-controlled nodata-height

docs/src/advanced/Algorithms.md

@@ -6,10 +6,11 @@ The algorithms are meant to overcome the limitation of `expression` (using [nume
 
 We added a set of custom algorithms:
 
-- `hillshade`: Create hillshade from elevation dataset
-- `contours`: Create contours lines (raster) from elevation dataset
-- `terrarium`: Mapzen's format to encode elevation value in RGB values (https://github.com/tilezen/joerd/blob/master/docs/formats.md#terrarium)
-- `terrainrgb`: Mapbox's format to encode elevation value in RGB values (https://docs.mapbox.com/data/tilesets/guides/access-elevation-data/)
+- `hillshade`: Create hillshade from elevation dataset (parameters: azimuth (45), angle_altitude(45))
+- `contours`: Create contours lines (raster) from elevation dataset (parameters: increment (35), thickness (1))
+- `slope`: Create degrees of slope from elevation dataset 
+- `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`
+- `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)`
 - `normalizedIndex`: Normalized Difference Index (e.g NDVI)
 - `cast`: Cast data to integer
 - `floor`: Round data to the smallest integer

src/titiler/core/tests/test_algorithms.py

@@ -236,6 +236,16 @@ def test_terrarium():
     assert out.array.dtype == "uint8"
     assert out.array[0, 0, 0] is numpy.ma.masked
 
+    # works on the above masked array img, with algo which was passed nodata_height
+    nodata_height = 10.0
+    algo = default_algorithms.get("terrarium")(nodata_height=nodata_height)
+    out = algo(img)
+    masked = out.array[:, arr.mask[0, :, :]]
+    masked_height = (masked[0] * 256 + masked[1] + masked[2] / 256) - 32768
+    numpy.testing.assert_array_equal(
+        masked_height, nodata_height * numpy.ones((100 * 100), dtype="bool")
+    )
+
 
 def test_terrainrgb():
     """test terrainrgb."""
@@ -259,6 +269,18 @@ def test_terrainrgb():
     assert out.array.dtype == "uint8"
     assert out.array[0, 0, 0] is numpy.ma.masked
 
+    # works on the above masked array img, with algo which was passed nodata_height
+    nodata_height = 10.0
+    algo = default_algorithms.get("terrainrgb")(nodata_height=nodata_height)
+    out = algo(img)
+    masked = out.array[:, arr.mask[0, :, :]]
+    masked_height = -10000 + (
+        ((masked[0] * 256 * 256) + (masked[1] * 256) + masked[2]) * 0.1
+    )
+    numpy.testing.assert_array_equal(
+        masked_height, nodata_height * numpy.ones((100 * 100), dtype="bool")
+    )
+
 
 def test_ops():
     """test ops: cast, ceil and floor."""

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

@@ -1,5 +1,7 @@
 """titiler.core.algorithm DEM."""
 
+from typing import Optional
+
 import numpy
 from pydantic import Field
 from rasterio import windows
@@ -22,6 +24,7 @@ class HillShade(BaseAlgorithm):
     azimuth: int = Field(45, ge=0, le=360)
     angle_altitude: float = Field(45.0, ge=-90.0, le=90.0)
     buffer: int = Field(3, ge=0, le=99)
+    z_exaggeration: float = Field(1.0, ge=1e-6, le=1e6)
 
     # metadata
     input_nbands: int = 1
@@ -31,6 +34,8 @@ class HillShade(BaseAlgorithm):
     def __call__(self, img: ImageData) -> ImageData:
         """Create hillshade from DEM dataset."""
         x, y = numpy.gradient(img.array[0])
+        x *= self.z_exaggeration
+        y *= self.z_exaggeration
         slope = numpy.pi / 2.0 - numpy.arctan(numpy.sqrt(x * x + y * y))
         aspect = numpy.arctan2(-x, y)
         azimuth = 360.0 - self.azimuth
@@ -71,6 +76,7 @@ class Slope(BaseAlgorithm):
 
     # parameters
     buffer: int = Field(3, ge=0, le=99, description="Buffer size for edge effects")
+    z_exaggeration: float = Field(1.0, ge=1e-6, le=1e6)
 
     # metadata
     input_nbands: int = 1
@@ -86,6 +92,8 @@ def __call__(self, img: ImageData) -> ImageData:
         pixel_size_y = abs(img.transform[4])
 
         x, y = numpy.gradient(img.array[0])
+        x *= self.z_exaggeration
+        y *= self.z_exaggeration
         dx = x / pixel_size_x
         dy = y / pixel_size_y
 
@@ -161,6 +169,7 @@ class Terrarium(BaseAlgorithm):
 
     title: str = "Terrarium"
     description: str = "Encode DEM into RGB (Mapzen Terrarium)."
+    nodata_height: Optional[float] = Field(None, ge=-99999.0, le=99999.0)
 
     # metadata
     input_nbands: int = 1
@@ -170,6 +179,10 @@ class Terrarium(BaseAlgorithm):
     def __call__(self, img: ImageData) -> ImageData:
         """Encode DEM into RGB."""
         data = numpy.clip(img.array[0] + 32768.0, 0.0, 65535.0)
+        if self.nodata_height is not None:
+            data[img.array.mask[0]] = numpy.clip(
+                self.nodata_height + 32768.0, 0.0, 65535.0
+            )
         r = data / 256
         g = data % 256
         b = (data * 256) % 256
@@ -191,6 +204,7 @@ class TerrainRGB(BaseAlgorithm):
     # parameters
     interval: float = Field(0.1, ge=0.0, le=1.0)
     baseval: float = Field(-10000.0, ge=-99999.0, le=99999.0)
+    nodata_height: Optional[float] = Field(None, ge=-99999.0, le=99999.0)
 
     # metadata
     input_nbands: int = 1
@@ -224,9 +238,15 @@ def _range_check(datarange):
         if _range_check(datarange):
             raise ValueError(f"Data of {datarange} larger than 256 ** 3")
 
-        r = ((((data // 256) // 256) / 256) - (((data // 256) // 256) // 256)) * 256
-        g = (((data // 256) / 256) - ((data // 256) // 256)) * 256
-        b = ((data / 256) - (data // 256)) * 256
+        if self.nodata_height is not None:
+            data[img.array.mask[0]] = (
+                self.nodata_height - self.baseval
+            ) / self.interval
+
+        data_int32 = data.astype(numpy.int32)
+        b = (data_int32) & 0xFF
+        g = (data_int32 >> 8) & 0xFF
+        r = (data_int32 >> 16) & 0xFF
 
         return ImageData(
             numpy.ma.stack([r, g, b]).astype(self.output_dtype),