diff --git a/docs/src/advanced/Algorithms.md b/docs/src/advanced/Algorithms.md index 8f65c7860..505ce57bb 100644 --- a/docs/src/advanced/Algorithms.md +++ b/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 diff --git a/src/titiler/core/tests/test_algorithms.py b/src/titiler/core/tests/test_algorithms.py index 3628e33f4..4be542cf5 100644 --- a/src/titiler/core/tests/test_algorithms.py +++ b/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.""" diff --git a/src/titiler/core/titiler/core/algorithm/dem.py b/src/titiler/core/titiler/core/algorithm/dem.py index 8c5a2ed2c..571315af5 100644 --- a/src/titiler/core/titiler/core/algorithm/dem.py +++ b/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),