[Experiment] sketch algorithm support

pyproject.toml

@@ -20,6 +20,7 @@ dynamic = ["version"]
 dependencies = [
     "braceexpand",
     "rio-cogeo>=3.1",
+    "rio-tiler>=3.1.5",
     "titiler.core>=0.5,<0.8",
     "starlette-cramjam>=0.3,<0.4",
     "uvicorn",

rio_viz/algorithm/__init__.py

@@ -0,0 +1,15 @@
+"""rio_viz.algorithm."""
+
+from typing import Dict, Type
+
+from rio_viz.algorithm.base import AlgorithmMetadata, BaseAlgorithm  # noqa
+from rio_viz.algorithm.dem import Contours, HillShade, TerrainRGB, Terrarium
+from rio_viz.algorithm.index import NormalizedIndex
+
+AVAILABLE_ALGORITHM: Dict[str, Type[BaseAlgorithm]] = {
+    "hillshade": HillShade,
+    "contours": Contours,
+    "normalizedIndex": NormalizedIndex,
+    "terrarium": Terrarium,
+    "terrainrgb": TerrainRGB,
+}

rio_viz/algorithm/base.py

@@ -0,0 +1,37 @@
+"""Algorithm base class."""
+
+import abc
+from typing import Dict, Optional, Sequence
+
+from pydantic import BaseModel
+from rio_tiler.models import ImageData
+
+
+class BaseAlgorithm(BaseModel, metaclass=abc.ABCMeta):
+    """Algorithm baseclass."""
+
+    input_nbands: int
+
+    output_nbands: int
+    output_dtype: str
+    output_min: Optional[Sequence]
+    output_max: Optional[Sequence]
+
+    @abc.abstractmethod
+    def apply(self, img: ImageData) -> ImageData:
+        """Apply"""
+        ...
+
+    class Config:
+        """Config for model."""
+
+        extra = "allow"
+
+
+class AlgorithmMetadata(BaseModel):
+    """Algorithm metadata."""
+
+    name: str
+    inputs: Dict
+    outputs: Dict
+    params: Dict

rio_viz/algorithm/dem.py

@@ -0,0 +1,163 @@
+"""rio_viz.algorithm DEM."""
+
+import numpy
+from rio_tiler.colormap import apply_cmap, cmap
+from rio_tiler.models import ImageData
+from rio_tiler.utils import linear_rescale
+
+from rio_viz.algorithm.base import BaseAlgorithm
+
+
+class HillShade(BaseAlgorithm):
+    """Hillshade."""
+
+    azimuth: int = 90
+    angle_altitude: float = 90
+
+    input_nbands: int = 1
+
+    output_nbands: int = 1
+    output_dtype: str = "uint8"
+
+    def apply(self, img: ImageData) -> ImageData:
+        """Create hillshade from DEM dataset."""
+        data = img.data[0]
+        mask = img.mask
+
+        x, y = numpy.gradient(data)
+
+        slope = numpy.pi / 2.0 - numpy.arctan(numpy.sqrt(x * x + y * y))
+        aspect = numpy.arctan2(-x, y)
+        azimuthrad = self.azimuth * numpy.pi / 180.0
+        altituderad = self.angle_altitude * numpy.pi / 180.0
+        shaded = numpy.sin(altituderad) * numpy.sin(slope) + numpy.cos(
+            altituderad
+        ) * numpy.cos(slope) * numpy.cos(azimuthrad - aspect)
+        hillshade_array = 255 * (shaded + 1) / 2
+
+        # ImageData only accept image in form of (count, height, width)
+        arr = numpy.expand_dims(hillshade_array, axis=0).astype(dtype=numpy.uint8)
+
+        return ImageData(
+            arr,
+            mask,
+            assets=img.assets,
+            crs=img.crs,
+            bounds=img.bounds,
+        )
+
+
+class Contours(BaseAlgorithm):
+    """Contours.
+
+    Original idea from https://custom-scripts.sentinel-hub.com/dem/contour-lines/
+    """
+
+    increment: int = 35
+    thickness: int = 1
+    minz: int = -12000
+    maxz: int = 8000
+
+    input_nbands: int = 1
+
+    output_nbands: int = 3
+    output_dtype: str = "uint8"
+
+    def apply(self, img: ImageData) -> ImageData:
+        """Add contours."""
+        data = img.data
+
+        # Apply rescaling for minz,maxz to 1->255 and apply Terrain colormap
+        arr = linear_rescale(data, (self.minz, self.maxz), (1, 255)).astype("uint8")
+        arr, _ = apply_cmap(arr, cmap.get("terrain"))
+
+        # set black (0) for contour lines
+        arr = numpy.where(data % self.increment < self.thickness, 0, arr)
+
+        return ImageData(
+            arr,
+            img.mask,
+            assets=img.assets,
+            crs=img.crs,
+            bounds=img.bounds,
+        )
+
+
+class Terrarium(BaseAlgorithm):
+    """Encode DEM into RGB (Mapzen Terrarium)."""
+
+    input_nbands: int = 1
+
+    output_nbands: int = 3
+    output_dtype: str = "uint8"
+
+    def apply(self, img: ImageData) -> ImageData:
+        """Encode DEM into RGB."""
+        data = numpy.clip(img.data[0] + 32768.0, 0.0, 65535.0)
+        r = data / 256
+        g = data % 256
+        b = (data * 256) % 256
+        arr = numpy.stack([r, g, b]).astype(numpy.uint8)
+        print(arr.shape)
+        return ImageData(
+            arr,
+            img.mask,
+            assets=img.assets,
+            crs=img.crs,
+            bounds=img.bounds,
+        )
+
+
+class TerrainRGB(BaseAlgorithm):
+    """Encode DEM into RGB (Mapbox Terrain RGB)."""
+
+    interval: int = 1
+    baseval: int = -10000
+
+    input_nbands: int = 1
+
+    output_nbands: int = 3
+    output_dtype: str = "uint8"
+
+    def apply(self, img: ImageData) -> ImageData:
+        """Encode DEM into RGB (Mapbox Terrain RGB).
+
+        Code from https://github.com/mapbox/rio-rgbify/blob/master/rio_rgbify/encoders.py (MIT)
+
+        """
+
+        def _range_check(datarange):
+            """
+            Utility to check if data range is outside of precision for 3 digit base 256
+            """
+            maxrange = 256**3
+
+            return datarange > maxrange
+
+        round_digits = 0
+
+        data = img.data[0].astype(numpy.float64)
+        data -= self.baseval
+        data /= self.interval
+
+        data = numpy.around(data / 2**round_digits) * 2**round_digits
+
+        rows, cols = data.shape
+        datarange = data.max() - data.min()
+        if _range_check(datarange):
+            raise ValueError("Data of {} larger than 256 ** 3".format(datarange))
+
+        rgb = numpy.zeros((3, rows, cols), dtype=numpy.uint8)
+        rgb[2] = ((data / 256) - (data // 256)) * 256
+        rgb[1] = (((data // 256) / 256) - ((data // 256) // 256)) * 256
+        rgb[0] = (
+            (((data // 256) // 256) / 256) - (((data // 256) // 256) // 256)
+        ) * 256
+
+        return ImageData(
+            rgb,
+            img.mask,
+            assets=img.assets,
+            crs=img.crs,
+            bounds=img.bounds,
+        )

rio_viz/algorithm/index.py

@@ -0,0 +1,37 @@
+"""rio_viz.algorithm Normalized Index."""
+
+from typing import Sequence
+
+import numpy
+from rio_tiler.models import ImageData
+
+from rio_viz.algorithm.base import BaseAlgorithm
+
+
+class NormalizedIndex(BaseAlgorithm):
+    """Normalized Difference Index."""
+
+    input_nbands: int = 2
+
+    output_nbands: int = 1
+    output_dtype: str = "float32"
+    output_min: Sequence[float] = [-1.0]
+    output_max: Sequence[float] = [1.0]
+
+    def apply(self, img: ImageData) -> ImageData:
+        """Normalized difference."""
+        b1 = img.data[0].astype("float32")
+        b2 = img.data[1].astype("float32")
+
+        arr = numpy.where(img.mask, (b2 - b1) / (b2 + b1), 0)
+
+        # ImageData only accept image in form of (count, height, width)
+        arr = numpy.expand_dims(arr, axis=0).astype(self.output_dtype)
+
+        return ImageData(
+            arr,
+            img.mask,
+            assets=img.assets,
+            crs=img.crs,
+            bounds=img.bounds,
+        )