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/titiler/core/algorithm/dem.py

@@ -22,6 +22,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., ge=1e-6, le=1e6)
 
     # metadata
     input_nbands: int = 1
@@ -31,6 +32,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 +74,7 @@ class Slope(BaseAlgorithm):
 
     # parameters
     buffer: int = Field(3, ge=0, le=99, description="Buffer size for edge effects")
+    z_exaggeration: float = Field(1., ge=1e-6, le=1e6)
 
     # metadata
     input_nbands: int = 1
@@ -86,6 +90,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 +167,8 @@ class Terrarium(BaseAlgorithm):
 
     title: str = "Terrarium"
     description: str = "Encode DEM into RGB (Mapzen Terrarium)."
+    use_nodata_height: bool = Field(False)
+    nodata_height: float = Field(0.0, ge=-99999.0, le=99999.0)
 
     # metadata
     input_nbands: int = 1
@@ -170,6 +178,8 @@ 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.use_nodata_height: 
+            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 +201,8 @@ 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)
+    use_nodata_height: bool = Field(False)
+    nodata_height: float = Field(0.0, ge=-99999.0, le=99999.0)
 
     # metadata
     input_nbands: int = 1
@@ -224,9 +236,13 @@ 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.use_nodata_height: 
+            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),