Terrain

Terrain visualisation — color relief, hill shade, slope, and aspect via GDAL DEMProcessing. Subclasses pyramids.dataset.Dataset, so all pyramids methods are inherited.

digitalrivers.terrain.Terrain

Bases: Dataset

Terrain analysis tools built on GDAL DEMProcessing.

Wraps a single- or multi-band raster and exposes convenience methods for visualisation (color_relief, hill_shade) and analysis: slope, aspect, the ruggedness derivatives (roughness, tpi, tri), and line-of-sight viewshed. Every method returns a pyramids.dataset.Dataset so results compose with the rest of the stack.

Parameters:

Name	Type	Description	Default
raster	Dataset	A gdal.Dataset to wrap. To open a terrain raster from a file path, use the inherited Terrain.read_file(path) classmethod.	required
access	str	"read_only" (default) or "write".	'read_only'

Examples:

• Wrap an in-memory DEM and compute a ruggedness derivative:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> arr = np.array(
  ...     [[10, 11, 12], [10, 9, 8], [5, 6, 30]], dtype=np.float32
  ... )
  >>> ds = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=1.0, epsg=32636,
  ...     no_data_value=-9999.0,
  ... )
  >>> terrain = Terrain(ds.raster)
  >>> terrain.roughness().read_array().shape
  (3, 3)
  

• Derive slope and read back its single float32 band:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> ramp = np.arange(9, dtype=np.float32).reshape(3, 3)
  >>> ds = Dataset.create_from_array(
  ...     ramp, top_left_corner=(0, 0), cell_size=1.0, epsg=32636,
  ...     no_data_value=-9999.0,
  ... )
  >>> slope = Terrain(ds.raster).slope()
  >>> slope.dtype
  ['float32']
  

See Also

digitalrivers.dem.DEM: Hydrological DEM processing (fill, flow direction, accumulation) with native window-configurable tpi / ruggedness.

Source code in src/digitalrivers/terrain.py

class Terrain(Dataset):
    """Terrain analysis tools built on GDAL `DEMProcessing`.

    Wraps a single- or multi-band raster and exposes convenience methods for
    visualisation (`color_relief`, `hill_shade`) and analysis: `slope`,
    `aspect`, the ruggedness derivatives (`roughness`, `tpi`, `tri`), and
    line-of-sight `viewshed`. Every method returns a `pyramids.dataset.Dataset`
    so results compose with the rest of the stack.

    Args:
        raster: A `gdal.Dataset` to wrap. To open a terrain raster from a file
            path, use the inherited `Terrain.read_file(path)` classmethod.
        access: `"read_only"` (default) or `"write"`.

    Examples:
        - Wrap an in-memory DEM and compute a ruggedness derivative:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> arr = np.array(
            ...     [[10, 11, 12], [10, 9, 8], [5, 6, 30]], dtype=np.float32
            ... )
            >>> ds = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=1.0, epsg=32636,
            ...     no_data_value=-9999.0,
            ... )
            >>> terrain = Terrain(ds.raster)
            >>> terrain.roughness().read_array().shape
            (3, 3)

            ```
        - Derive slope and read back its single float32 band:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> ramp = np.arange(9, dtype=np.float32).reshape(3, 3)
            >>> ds = Dataset.create_from_array(
            ...     ramp, top_left_corner=(0, 0), cell_size=1.0, epsg=32636,
            ...     no_data_value=-9999.0,
            ... )
            >>> slope = Terrain(ds.raster).slope()
            >>> slope.dtype
            ['float32']

            ```

    See Also:
        digitalrivers.dem.DEM: Hydrological DEM processing (fill, flow direction,
            accumulation) with native window-configurable `tpi` / `ruggedness`.
    """

    def __init__(self, raster: gdal.Dataset, access: str = "read_only"):
        """Wrap a GDAL dataset for terrain analysis.

        Args:
            raster: A `gdal.Dataset` to wrap. To open a terrain raster from a
                file path use the inherited `Terrain.read_file(path)`
                classmethod; the bare constructor accepts a `gdal.Dataset`
                only.
            access: `"read_only"` (default) or `"write"`.

        Examples:
            - Wrap an in-memory raster and read its grid dimensions:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from digitalrivers import Terrain
                >>> arr = np.ones((4, 5), dtype=np.float32)
                >>> ds = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=1.0, epsg=4326,
                ... )
                >>> terrain = Terrain(ds.raster)
                >>> terrain.shape[-2:]
                (4, 5)

                ```
        """
        super().__init__(raster, access)

    def color_relief(
        self,
        band: int = 0,
        path: str | None = None,
        color_table: DataFrame | None = None,
        **kwargs,
    ) -> "Dataset":
        """Create a color relief for a band in the Dataset.

        A color relief raster is a raster image where each pixel's value is mapped to a specific color based on a
        predefined color palette or color table.

        Args:
            band: int, default is 0.
                band index.
            path: str, default is None.
                path to save the color relief raster.
            color_table: DataFrame, default is None.
                DataFrame with columns: band, values, color
                    ```text
                      values    color
                    0      1  #709959
                    1      2  #F2EEA2
                    2      3  #F2CE85
                    3      1  #C28C7C
                    4      2  #D6C19C
                    5      3  #D6C19C
                    ```
                or DataFrame with columns: values, red, green, blue, alpha, (the alpha column is optional)
                    ```text
                      values    red  green   blue  alpha
                    0      1    112    153     89    255
                    1      2    242    238    162    255
                    2      3    242    206    133    255
                    3      1    194    140    124    255
                    4      2    214    193    156    255
                    5      3    214    193    156    255
                    ```
        Returns:
            Dataset:
                Dataset with the color relief with four bands read, green, blue, and alpha.

        Examples:
            - First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.
                ```python
                >>> import numpy as np
                >>> arr = np.random.randint(0, 15, size=(10, 10))
                >>> dataset = Dataset.create_from_array(arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326)

                ```
            - Now let's create the color table using hex colors.
                ```python
                >>> import pandas as pd
                >>> color_hex = ["#709959", "#F2EEA2", "#F2CE85", "#C28C7C", "#D6C19C"]
                >>> values = [1, 3, 5, 7, 9]
                >>> df = pd.DataFrame(columns=["values", "color"])
                >>> df.loc[:, "values"] = values
                >>> df.loc[:, "color"] = color_hex

                ```
            - Now let's create the color relief for the dataset using the color table `DataFrame`.
                ```python
                >>> color_relief = Terrain(dataset.raster).color_relief(band=0, color_table=df)
                >>> print(color_relief) # doctest: +SKIP
                <BLANKLINE>
                            Cell size: 0.05
                            Dimension: 10 * 10
                            EPSG: 4326
                            Number of Bands: 4
                            Band names: ['Band_1', 'Band_2', 'Band_3', 'Band_4']
                            Mask: None
                            Data type: byte
                            projection: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
                            Metadata: {}
                            File: ...
                <BLANKLINE>
                >>> print(color_relief.band_color)
                {0: 'red', 1: 'green', 2: 'blue', 3: 'alpha'}

                ```
            - The result color relief dataset will have 4 bands red, green, blue, and alpha. with values from 0 to 255.
            - To plot the color relief dataset, you can use the `plot` method. but you need to provide the the rgb indices
                with the alpha index as the fourth index, otherwise the alpha band will be missing.
                ```python
                >>> fig, ax = color_relief.plot(rgb=[0, 1, 2, 3]) # doctest: +SKIP

                ```
            ![color-relief](./../_images/dataset/color-relief.png)

        See Also:
            Dataset.hill_shade: create a hill-shade for a band in the Dataset.
        """
        if path is None:
            driver = "MEM"
            path = ""
        else:
            driver = "GTiff"
        color_df = self.analysis._process_color_table(color_table)

        with tempfile.TemporaryDirectory() as temp_dir:
            color_table_path = os.path.join(temp_dir, f"{uuid.uuid1()}.txt")
            color_df.to_csv(color_table_path, index=False, header=False)

            options = gdal.DEMProcessingOptions(
                band=band + 1,
                format=driver,
                colorFilename=color_table_path,
                addAlpha=True,
                creationOptions=CREATION_OPTIONS,
                **kwargs,
            )
            dst = gdal.DEMProcessing(path, self.raster, "color-relief", options=options)

        color_relief = Dataset(dst, access="write")
        color_relief.band_color = {0: "red", 1: "green", 2: "blue", 3: "alpha"}
        return color_relief

    def hill_shade(
        self,
        band: int = 0,
        azimuth: int | float | list[int] = 315,
        altitude: int | float | list[int] = 45,
        vertical_exaggeration: int | float | list[int] = 1,
        scale: int | float | list[int] = 1,
        path: str | None = None,
        weights: list[int] | None = None,
        **kwargs,
    ) -> "Dataset":
        """Create hill-shade.

        Hillshade is a technique used in digital elevation modeling (DEM) to create a grayscale representation of a
        terrain's surface that simulates the effect of sunlight falling across the landscape.
        This technique helps to visualize the shape and features of the terrain by highlighting the variations in
        elevation and the slope of the surface.

        Hillshade calculates the illumination of each pixel based on the slope (gradient) and aspect (direction) of the
        terrain surface relative to a specified light source.

        The main parameters influencing the hillshade effect are:
        - Light source direction (Azimuth): the azimuth angle of the light source, which is the angle between the light
            source
        - Light source elevation (altitude): the source of light elevation, it is measured in degrees from the horizon.
        - Vertical exaggeration (Z-factor): the vertical exaggeration is used to emphasize the vertical features of the
            terrain.

        Notes:
            if the `hill_shade` parameters are given as lists then the hill shade will be calculated for each set
            of parameter and then the average will be returned.

        Args:
            band: int
                band index.
            azimuth: int | float | list[int]
                The source of light direction, it is measured clockwise from the north. zero means from north to south.
                45 degrees means from the northeast to the southwest.
            altitude: int | float | list[int]
                The source of light elevation, it is measured in degrees from the horizon. zero means from the horizon.
                90 degrees means from the zenith.
                the overall image gets brighter as the light source gets closer to the zenith. The brightest slopes/DEM
                features will be perpendicular to the light source, and the darkest will be angled 90˚ or more away.
            vertical_exaggeration: int | float | list[int]
                Vertical exaggeration, the vertical exaggeration It is used to emphasize the
                vertical features of the terrain.
            scale: int | float | list[int]
                the scale is the ratio of vertical units to horizontal. If the horizontal unit of the source DEM is
                degrees (e.g Lat/Long WGS84 projection), you can use scale=111120 if the vertical units are meters
                (or scale=370400 if they are in feet).
            path: str, optional, default is None
                path to save the hill-shade raster.
            weights: list[int], default is None.
                list of weights to combine the hill-shades if the other parameters are given as lists, an average hill
                shade will be calculated based on the weights. if None, the weights will be equal.
            **kwargs:
                multi_directional: bool
                    if True, the hill shade will be calculated for multiple azimuth values [225, 270, 315, 360] each with a
                    altitude of 30 degrees, and then the average will be returned. with multi_directional = True any given
                    azimuth will be ignored.
                    For more details visit: https://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
                combined: bool
                    combined shading, a combination of slope and oblique shading.
                igor: bool
                    shading which tries to minimize effects on other map features beneath. with `igor=True` the altitude
                    will be calculated ignored.
                    For more details visit: https://maperitive.net/docs/Commands/GenerateReliefImageIgor.html

        Returns:
            Dataset: 8-bit
                Dataset with the hill-shade created.

        Examples:
            - First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.
                ```python
                >>> import numpy as np
                >>> arr = np.random.randint(0, 15, size=(100, 100))
                >>> dataset = Dataset.create_from_array(arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326)

                >>> hill_shade = Terrain(dataset.raster).hill_shade(
                ...     band=0, altitude=45, azimuth=315, vertical_exaggeration=1, scale=1
                ... )

                >>> print(hill_shade.dtype) # doctest: +SKIP
                ['byte']
                >>> hill_shade.plot() # doctest: +SKIP

                ```
                ![hill-shade](./../_images/dataset/hill-shade.png)
                ```python
                >>> hill_shade.stats() # doctest: +SKIP
                        min    max       mean        std
                Band_1  1.0  223.0  58.880951  71.079056

                ```
            - You can also provide the function with a list os values for each parameter, then the functions will
                calculate the hill shade for each set of parameters and then the average will be returned.
                ```python
                >>> hill_shade = Terrain(dataset.raster).hill_shade(
                ...     band=0, azimuth=[315, 45], altitude=[45, 45], vertical_exaggeration=[1, 1], scale=[1, 1]
                ... )

                >>> hill_shade.plot() # doctest: +SKIP

                ```
                ![hill-shade-multi](./../_images/dataset/hill-shade-multi.png)

        See Also:
            Dataset.color_relief: create a color relief for a band in the Dataset.
            Dataset.slope: create a slope for a band in the Dataset.
        """
        if "multi_directional" in kwargs:
            if not isinstance(kwargs["multi_directional"], bool):
                raise ValueError("The multi_directional parameter must be a boolean.")
            if kwargs["multi_directional"]:
                multi_directional = True
                azimuth = None
                # altitude, vertical_exaggeration, scale = None, None, None,
            else:
                multi_directional = False

            kwargs.pop("multi_directional")
            kwargs["multiDirectional"] = multi_directional
        if "igor" in kwargs:
            if not isinstance(kwargs["igor"], bool):
                raise ValueError("The igor parameter must be a boolean.")
            if kwargs["igor"]:
                altitude = None

        # if not (
        #     type(azimuth)
        #     is type(altitude)
        #     is type(vertical_exaggeration)
        #     is type(scale)
        # ):
        #     raise ValueError(
        #         f"The azimuth, altitude, vertical_exaggeration, and scale parameter must be of the same type. Given"
        #         f" azimuth: {type(azimuth)}, altitude: {type(altitude)}, vertical_exaggeration: {type(vertical_exaggeration)},"
        #         f"scale: {type(scale)}"
        #     )

        if path is None:
            driver = "MEM"
            path = ""
        else:
            driver = "GTiff"

        wrap = lambda v: v if isinstance(v, list) else [v]
        azimuth, altitude, vertical_exaggeration, scale = (
            wrap(azimuth),
            wrap(altitude),
            wrap(vertical_exaggeration),
            wrap(scale),
        )
        if not (
            len(azimuth) == len(altitude) == len(vertical_exaggeration) == len(scale)
        ):
            raise ValueError("All list parameters must have the same length.")

        # get the hill shade for all the parameters
        hill_shades: list[gdal.Dataset] = []
        for az, alt, ver_ex, scale_1 in zip(
            azimuth, altitude, vertical_exaggeration, scale
        ):
            dst = self._create_hill_shade(
                band, driver, az, alt, ver_ex, scale_1, path, **kwargs
            )
            hill_shades.append(dst)

        if len(hill_shades) > 1:
            if weights is None:
                weights = np.ones(len(azimuth))
            weights = np.array(weights) / np.sum(weights)
            hill_shades_arr: list[np.ndarray] = [
                hill_shade.ReadAsArray() for hill_shade in hill_shades
            ]
            combined_hillshade = np.average(hill_shades_arr, axis=0, weights=weights)
            combined_hillshade = np.clip(combined_hillshade, 0, 255).astype(np.uint8)
            hill_shade = Dataset.dataset_like(
                Dataset(hill_shades[0]), combined_hillshade
            )
        else:
            hill_shade = Dataset(hill_shades[0], access="write")

        hill_shade.band_color = {0: "gray_index"}

        return hill_shade

    def _create_hill_shade(
        self,
        band: int,
        driver: str,
        azimuth: int | float = 315,
        altitude: int | float = 45,
        vertical_exaggeration: int | float = 1,
        scale: int | float = 1,
        path: str | None = None,
        **kwargs,
    ) -> gdal.Dataset:
        """Run a single GDAL `DEMProcessing("hillshade")` call.

        Args:
            band: Zero-based band index.
            driver: GDAL driver name (`"MEM"` or `"GTiff"`).
            azimuth: Light-source azimuth in degrees clockwise from
                north.
            altitude: Light-source elevation in degrees above horizon.
            vertical_exaggeration: Z-factor for vertical emphasis.
            scale: Ratio of vertical to horizontal units.
            path: Output file path (empty string for in-memory).
            **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

        Returns:
            gdal.Dataset: Raw GDAL dataset with the computed hill shade.
        """
        options = gdal.DEMProcessingOptions(
            band=band + 1,
            format=driver,
            azimuth=azimuth,
            altitude=altitude,
            zFactor=vertical_exaggeration,
            scale=scale,
            creationOptions=["COMPRESS=LZW"],
            **kwargs,
        )
        dst = gdal.DEMProcessing(path, self.raster, "hillshade", options=options)

        return dst

    def slope(
        self,
        band: int = 0,
        scale: int | float | list[int] = 1,
        slope_format: str = "degree",
        path: str | None = None,
        algorithm: str | None = None,
        creation_options: list[str] | None = None,
        **kwargs,
    ) -> "Dataset":
        """Compute the slope of the terrain surface.

        Uses GDAL `DEMProcessing` to calculate the slope (rate of
        elevation change) for every cell.

        Args:
            band: Zero-based band index. Defaults to 0.
            scale: Ratio of vertical to horizontal units.  Use
                `111120` when the horizontal CRS is in degrees and
                vertical units are metres.  Defaults to 1.
            slope_format: Output format — `"degree"` (default) or
                `"percent"`.
            algorithm: Slope algorithm.  One of `"Horn"`,
                `"ZevenbergenThorne"`, or `None` (GDAL default).
                Zevenbergen-Thorne suits smooth landscapes; Horn
                performs better on rough terrain.
            path: If given, write the result to this GeoTIFF path.
                Otherwise the raster is created in memory.
            creation_options: GDAL creation options.  Defaults to
                `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
            **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

        Returns:
            Dataset: Single-band `float32` raster with slope values.
                No-data value is `-9999.0`.

        Examples:
            - First create a one band dataset, consisting of 10 columns
                and 10 rows, with random values between 0 and 15.
                ```python
                >>> import numpy as np
                >>> arr = np.random.randint(0, 15, size=(10, 10))
                >>> dataset = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326
                ... )

                ```
            - Now let's create the slope for the dataset.
                ```python
                >>> slope = Terrain(dataset.raster).slope()
                >>> fig, ax = slope.plot() # doctest: +SKIP

                ```
                ![slope](./../_images/dataset/slope.png)

        See Also:
            Terrain.hill_shade: Create a hill-shade for a band in the
                Dataset.
            Terrain.color_relief: Create a color relief for a band in
                the Dataset.
        """
        if path is None:
            driver = "MEM"
            path = ""
        else:
            driver = "GTiff"

        if creation_options is None:
            creation_options = CREATION_OPTIONS.copy()

        options = gdal.DEMProcessingOptions(
            band=band + 1,
            format=driver,
            alg=algorithm,
            slopeFormat=slope_format,
            scale=scale,
            creationOptions=creation_options,
            **kwargs,
        )
        dst = gdal.DEMProcessing(path, self.raster, "slope", options=options)
        src = Dataset(dst, access="write")

        return src

    def aspect(
        self,
        band: int = 0,
        scale: int | float | list[int] = 1,
        vertical_exaggeration: int | float | list[int] = 1,
        zero_flat_surface: bool = False,
        algorithm: str | None = None,
        path: str | None = None,
        creation_options: list[str] | None = None,
        **kwargs,
    ) -> "Dataset":
        """Compute the aspect (slope direction) of the terrain surface.

        Uses GDAL `DEMProcessing` to calculate the compass direction
        of the steepest downhill slope for every cell.  Values range
        from 0° (north) clockwise to 360°.

        Args:
            band: Zero-based band index. Defaults to 0.
            scale: Accepted for signature symmetry with `slope` but **not
                used** — GDAL's aspect mode is scale-invariant (it reports a
                direction) and rejects the `-s` option. Defaults to 1.
            vertical_exaggeration: Accepted for signature symmetry but **not
                used** for the same reason (aspect rejects the `-z` option).
                Defaults to 1.
            zero_flat_surface: If `True` flat areas get an aspect of
                0°.  If `False` (default) flat areas receive the
                no-data value.
            algorithm: Aspect algorithm.  One of `"Horn"`,
                `"ZevenbergenThorne"`, or `None` (GDAL default).
            path: If given, write the result to this GeoTIFF path.
                Otherwise the raster is created in memory.
            creation_options: GDAL creation options.  Defaults to
                `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
            **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

        Returns:
            Dataset: Single-band `float32` raster with aspect values
                in degrees (0–360).  No-data value is `-9999.0`.

        Examples:
            - Create a small raster and compute its aspect.
                ```python
                >>> import numpy as np
                >>> arr = np.random.randint(0, 15, size=(10, 10))
                >>> dataset = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326
                ... )

                ```
            - Compute the aspect raster.
                ```python
                >>> aspect = Terrain(dataset.raster).aspect()
                >>> fig, ax = aspect.plot() # doctest: +SKIP

                ```
                ![aspect](./../_images/dataset/aspect.png)

        See Also:
            Terrain.hill_shade: Create a hill-shade for a band in the
                Dataset.
            Terrain.slope: Compute the slope of the terrain surface.
        """
        if path is None:
            driver = "MEM"
            path = ""
        else:
            driver = "GTiff"

        if creation_options is None:
            creation_options = CREATION_OPTIONS.copy()

        # `scale` / `vertical_exaggeration` are NOT forwarded: GDAL's aspect mode
        # does not accept the slope/hillshade `-s` / `-z` options (aspect is a
        # scale-invariant direction). Passing them is silently tolerated only
        # when `zeroForFlat` is unset, and otherwise raises "Argument(s) are not
        # valid with any processing mode" — so they are dropped here.
        del scale, vertical_exaggeration
        options = gdal.DEMProcessingOptions(
            band=band + 1,
            format=driver,
            alg=algorithm,
            zeroForFlat=zero_flat_surface,
            creationOptions=creation_options,
            **kwargs,
        )
        dst = gdal.DEMProcessing(path, self.raster, "aspect", options=options)
        src = Dataset(dst, access="write")

        return src

    def _ruggedness(
        self,
        mode: str,
        band: int,
        path: str | None,
        compute_edges: bool,
        creation_options: list[str] | None,
        **kwargs,
    ) -> "Dataset":
        """Run a single GDAL `DEMProcessing` ruggedness mode.

        Shared backend for `roughness`, `tpi`, and `tri` — the three
        `gdaldem` ruggedness derivatives differ only by their GDAL mode
        string, so the driver / options plumbing lives here once.

        Args:
            mode: GDAL `DEMProcessing` mode — `"Roughness"`, `"TPI"`, or
                `"TRI"`.
            band: Zero-based band index.
            path: Output GeoTIFF path, or `None` for an in-memory raster.
            compute_edges: When `True` GDAL also computes values for the
                raster's edge cells (using the available partial window)
                instead of leaving them as no-data.
            creation_options: GDAL creation options. `None` falls back to
                `CREATION_OPTIONS`.
            **kwargs: Forwarded to `gdal.DEMProcessingOptions` (e.g.
                `alg` for `TRI`).

        Returns:
            Dataset: Single-band `float32` raster. No-data value is
                `-9999.0`.
        """
        if path is None:
            driver = "MEM"
            path = ""
        else:
            driver = "GTiff"

        if creation_options is None:
            creation_options = CREATION_OPTIONS.copy()

        options = gdal.DEMProcessingOptions(
            band=band + 1,
            format=driver,
            computeEdges=compute_edges,
            creationOptions=creation_options,
            **kwargs,
        )
        dst = gdal.DEMProcessing(path, self.raster, mode, options=options)
        return Dataset(dst, access="write")

    def roughness(
        self,
        band: int = 0,
        path: str | None = None,
        compute_edges: bool = False,
        creation_options: list[str] | None = None,
        **kwargs,
    ) -> "Dataset":
        """Compute terrain roughness — the largest elevation difference in a 3x3 window.

        Roughness (Wilson et al., 2007) is the maximum absolute difference
        between a cell and its eight neighbours. It is the simplest
        ruggedness measure and reacts strongly to local relief: flat
        surfaces score `0`, cliffs and noisy LiDAR returns score high.

        Args:
            band: Zero-based band index. Defaults to 0.
            path: If given, write the result to this GeoTIFF path.
                Otherwise the raster is created in memory.
            compute_edges: If `True`, compute values for the edge cells
                from the available partial window instead of leaving them
                no-data. Defaults to `False`.
            creation_options: GDAL creation options. Defaults to
                `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
            **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

        Returns:
            Dataset: Single-band `float32` raster of roughness values in
                the DEM's vertical units. No-data value is `-9999.0`.

        Examples:
            - Compute roughness for a small elevation raster.
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from digitalrivers import Terrain
                >>> arr = np.array(
                ...     [[10, 11, 12, 40], [10, 9, 8, 7],
                ...      [5, 6, 30, 6], [4, 3, 2, 1]],
                ...     dtype=np.float32,
                ... )
                >>> ds = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=1.0,
                ...     epsg=32636, no_data_value=-9999.0,
                ... )
                >>> roughness = Terrain(ds.raster).roughness()
                >>> roughness.read_array().shape
                (4, 4)

                ```

        See Also:
            Terrain.tpi: Topographic Position Index.
            Terrain.tri: Terrain Ruggedness Index.
        """
        return self._ruggedness(
            "Roughness", band, path, compute_edges, creation_options, **kwargs
        )

    def tpi(
        self,
        band: int = 0,
        path: str | None = None,
        compute_edges: bool = False,
        creation_options: list[str] | None = None,
        **kwargs,
    ) -> "Dataset":
        """Compute the Topographic Position Index (TPI).

        TPI (Weiss, 2001) is each cell's elevation minus the mean
        elevation of its eight neighbours. Positive values mark local
        highs (ridges, peaks), negative values mark local lows (valleys,
        channels), and values near zero mark flat areas or constant
        slopes. It is widely used for landform classification.

        Note:
            This is the GDAL formulation — the focal mean is taken over
            the eight neighbours **excluding** the centre cell, on a fixed
            3x3 window. `DEM.tpi` is a native alternative whose focal mean
            **includes** the centre cell and accepts an arbitrary
            `window` size, so the two return slightly different values.

        Args:
            band: Zero-based band index. Defaults to 0.
            path: If given, write the result to this GeoTIFF path.
                Otherwise the raster is created in memory.
            compute_edges: If `True`, compute values for the edge cells
                from the available partial window instead of leaving them
                no-data. Defaults to `False`.
            creation_options: GDAL creation options. Defaults to
                `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
            **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

        Returns:
            Dataset: Single-band `float32` raster of TPI values (signed,
                in the DEM's vertical units). No-data value is `-9999.0`.

        Examples:
            - Compute TPI for a small elevation raster.
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from digitalrivers import Terrain
                >>> arr = np.array(
                ...     [[10, 11, 12, 40], [10, 9, 8, 7],
                ...      [5, 6, 30, 6], [4, 3, 2, 1]],
                ...     dtype=np.float32,
                ... )
                >>> ds = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=1.0,
                ...     epsg=32636, no_data_value=-9999.0,
                ... )
                >>> tpi = Terrain(ds.raster).tpi()
                >>> tpi.read_array().shape
                (4, 4)

                ```

        See Also:
            Terrain.roughness: Maximum 3x3 elevation difference.
            Terrain.tri: Terrain Ruggedness Index.
            DEM.tpi: Native, window-configurable TPI (includes the centre
                cell in the focal mean).
        """
        return self._ruggedness(
            "TPI", band, path, compute_edges, creation_options, **kwargs
        )

    def tri(
        self,
        band: int = 0,
        algorithm: str | None = None,
        path: str | None = None,
        compute_edges: bool = False,
        creation_options: list[str] | None = None,
        **kwargs,
    ) -> "Dataset":
        """Compute the Terrain Ruggedness Index (TRI).

        TRI is the mean absolute difference between a cell and its eight
        neighbours. Two formulations are available via `algorithm`:

        * `"Riley"` (Riley et al., 1999) — square-root of the summed
          squared differences; the original TRI.
        * `"Wilson"` (Wilson et al., 2007) — the mean absolute
          difference; better suited to bathymetric / continuous data.

        Note:
            With `algorithm=None` GDAL uses the Riley root-sum-square
            form. The native `DEM.ruggedness` computes the Wilson
            mean-absolute-difference form, so it corresponds to
            `tri(algorithm="Wilson")` (on a 3x3 window) rather than the
            default here.

        Args:
            band: Zero-based band index. Defaults to 0.
            algorithm: TRI formulation — `"Riley"`, `"Wilson"`, or
                `None` (GDAL default). Defaults to `None`.
            path: If given, write the result to this GeoTIFF path.
                Otherwise the raster is created in memory.
            compute_edges: If `True`, compute values for the edge cells
                from the available partial window instead of leaving them
                no-data. Defaults to `False`.
            creation_options: GDAL creation options. Defaults to
                `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
            **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

        Returns:
            Dataset: Single-band `float32` raster of TRI values in the
                DEM's vertical units. No-data value is `-9999.0`.

        Examples:
            - Compute TRI for a small elevation raster.
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from digitalrivers import Terrain
                >>> arr = np.array(
                ...     [[10, 11, 12, 40], [10, 9, 8, 7],
                ...      [5, 6, 30, 6], [4, 3, 2, 1]],
                ...     dtype=np.float32,
                ... )
                >>> ds = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=1.0,
                ...     epsg=32636, no_data_value=-9999.0,
                ... )
                >>> tri = Terrain(ds.raster).tri()
                >>> tri.read_array().shape
                (4, 4)

                ```

        See Also:
            Terrain.roughness: Maximum 3x3 elevation difference.
            Terrain.tpi: Topographic Position Index.
            DEM.ruggedness: Native, window-configurable Wilson-form TRI.
        """
        if algorithm is not None:
            kwargs["alg"] = algorithm
        return self._ruggedness(
            "TRI", band, path, compute_edges, creation_options, **kwargs
        )

    def viewshed(
        self,
        observer_x: float,
        observer_y: float,
        band: int = 0,
        observer_height: float = 1.75,
        target_height: float = 0.0,
        max_distance: float = 0.0,
        mode: str = "max",
        visible_value: float = 255.0,
        invisible_value: float = 0.0,
        out_of_range_value: float = 0.0,
        no_data_value: float = -1.0,
        curvature_coefficient: float = 0.85714,
        path: str | None = None,
        creation_options: list[str] | None = None,
    ) -> "Dataset":
        """Compute the viewshed (line-of-sight visibility) from an observer point.

        Wraps GDAL `ViewshedGenerate` to flag, for every cell, whether it
        is visible from an observer standing at `(observer_x, observer_y)`,
        accounting for the intervening terrain. The observer and target
        heights are added above the DEM surface, and Earth curvature /
        atmospheric refraction can be modelled via `curvature_coefficient`.

        Args:
            observer_x: Observer X coordinate, in the DEM's CRS.
            observer_y: Observer Y coordinate, in the DEM's CRS.
            band: Zero-based band index of the elevation band. Defaults
                to 0.
            observer_height: Observer height above the DEM surface, in the
                DEM's vertical units. Defaults to 1.75 (roughly eye level).
            target_height: Target height above the DEM surface that must be
                visible. Defaults to 0.0 (ground level).
            max_distance: Maximum line-of-sight distance in CRS units.
                `0.0` (default) means unlimited.
            mode: Cell-evaluation method — `"max"` (default), `"min"`,
                `"edge"`, or `"diagonal"` — mapping to the GDAL
                `GVM_Max` / `GVM_Min` / `GVM_Edge` / `GVM_Diagonal`
                viewshed modes.
            visible_value: Output value written to visible cells. Defaults
                to 255.0.
            invisible_value: Output value written to hidden cells. Defaults
                to 0.0.
            out_of_range_value: Output value for cells beyond
                `max_distance`. Defaults to 0.0.
            no_data_value: Output no-data value. Defaults to -1.0.
            curvature_coefficient: Earth-curvature / refraction
                coefficient. Defaults to 0.85714 (GDAL's standard
                atmospheric value); use 1.0 to ignore curvature.
            path: If given, write the result to this GeoTIFF path.
                Otherwise the raster is created in memory.
            creation_options: GDAL creation options. Defaults to
                `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.

        Returns:
            Dataset: Single-band raster encoding visibility (`visible_value`
                for visible cells, `invisible_value` otherwise).

        Raises:
            ValueError: If `mode` is not one of `"max"`, `"min"`,
                `"edge"`, `"diagonal"`.

        Examples:
            - Compute the viewshed from the top-left corner of a small DEM.
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from digitalrivers import Terrain
                >>> arr = np.array(
                ...     [[10, 11, 12, 40], [10, 9, 8, 7],
                ...      [5, 6, 30, 6], [4, 3, 2, 1]],
                ...     dtype=np.float32,
                ... )
                >>> ds = Dataset.create_from_array(
                ...     arr, top_left_corner=(0, 0), cell_size=1.0,
                ...     epsg=32636, no_data_value=-9999.0,
                ... )
                >>> vs = Terrain(ds.raster).viewshed(
                ...     observer_x=0.5, observer_y=-0.5,
                ... )
                >>> vs.read_array().shape
                (4, 4)

                ```

        See Also:
            Terrain.hill_shade: Shaded-relief visualisation of the surface.
        """
        modes = {
            "diagonal": gdal.GVM_Diagonal,
            "edge": gdal.GVM_Edge,
            "max": gdal.GVM_Max,
            "min": gdal.GVM_Min,
        }
        if mode not in modes:
            raise ValueError(f"mode must be one of {sorted(modes)}; got {mode!r}")

        if path is None:
            driver = "MEM"
            path = ""
        else:
            driver = "GTiff"

        if creation_options is None:
            creation_options = CREATION_OPTIONS.copy()

        src_band = self.raster.GetRasterBand(band + 1)
        # `gdal.ViewshedGenerate` + the `GVM_*` enums require a recent GDAL; that
        # floor is supplied transitively by the `pyramids-gis >=0.29.0` pin (see
        # the dependency note in pyproject.toml), not a direct GDAL dependency.
        dst = gdal.ViewshedGenerate(
            src_band,
            driver,
            path,
            creation_options,
            observer_x,
            observer_y,
            observer_height,
            target_height,
            visible_value,
            invisible_value,
            out_of_range_value,
            no_data_value,
            curvature_coefficient,
            modes[mode],
            max_distance,
            heightMode=gdal.GVOT_NORMAL,
        )
        return Dataset(dst, access="write")

__init__(raster, access='read_only')

Wrap a GDAL dataset for terrain analysis.

Parameters:

Name	Type	Description	Default
raster	Dataset	A gdal.Dataset to wrap. To open a terrain raster from a file path use the inherited Terrain.read_file(path) classmethod; the bare constructor accepts a gdal.Dataset only.	required
access	str	"read_only" (default) or "write".	'read_only'

Examples:

• Wrap an in-memory raster and read its grid dimensions:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> arr = np.ones((4, 5), dtype=np.float32)
  >>> ds = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=1.0, epsg=4326,
  ... )
  >>> terrain = Terrain(ds.raster)
  >>> terrain.shape[-2:]
  (4, 5)
  

Source code in src/digitalrivers/terrain.py

def __init__(self, raster: gdal.Dataset, access: str = "read_only"):
    """Wrap a GDAL dataset for terrain analysis.

    Args:
        raster: A `gdal.Dataset` to wrap. To open a terrain raster from a
            file path use the inherited `Terrain.read_file(path)`
            classmethod; the bare constructor accepts a `gdal.Dataset`
            only.
        access: `"read_only"` (default) or `"write"`.

    Examples:
        - Wrap an in-memory raster and read its grid dimensions:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> arr = np.ones((4, 5), dtype=np.float32)
            >>> ds = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=1.0, epsg=4326,
            ... )
            >>> terrain = Terrain(ds.raster)
            >>> terrain.shape[-2:]
            (4, 5)

            ```
    """
    super().__init__(raster, access)

color_relief(band=0, path=None, color_table=None, **kwargs)

Create a color relief for a band in the Dataset.

A color relief raster is a raster image where each pixel's value is mapped to a specific color based on a predefined color palette or color table.

Parameters:

Name	Type	Description	Default
band	int	int, default is 0. band index.	0
path	str | None	str, default is None. path to save the color relief raster.	None
color_table	DataFrame | None	DataFrame, default is None. DataFrame with columns: band, values, color values color 0 1 #709959 1 2 #F2EEA2 2 3 #F2CE85 3 1 #C28C7C 4 2 #D6C19C 5 3 #D6C19C or DataFrame with columns: values, red, green, blue, alpha, (the alpha column is optional) values red green blue alpha 0 1 112 153 89 255 1 2 242 238 162 255 2 3 242 206 133 255 3 1 194 140 124 255 4 2 214 193 156 255 5 3 214 193 156 255	None

Returns: Dataset: Dataset with the color relief with four bands read, green, blue, and alpha.

Examples:

• First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.

  >>> import numpy as np
  >>> arr = np.random.randint(0, 15, size=(10, 10))
  >>> dataset = Dataset.create_from_array(arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326)
  

• Now let's create the color table using hex colors.

  >>> import pandas as pd
  >>> color_hex = ["#709959", "#F2EEA2", "#F2CE85", "#C28C7C", "#D6C19C"]
  >>> values = [1, 3, 5, 7, 9]
  >>> df = pd.DataFrame(columns=["values", "color"])
  >>> df.loc[:, "values"] = values
  >>> df.loc[:, "color"] = color_hex
  

• Now let's create the color relief for the dataset using the color table DataFrame.

  >>> color_relief = Terrain(dataset.raster).color_relief(band=0, color_table=df)
  >>> print(color_relief) # doctest: +SKIP
  <BLANKLINE>
              Cell size: 0.05
              Dimension: 10 * 10
              EPSG: 4326
              Number of Bands: 4
              Band names: ['Band_1', 'Band_2', 'Band_3', 'Band_4']
              Mask: None
              Data type: byte
              projection: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
              Metadata: {}
              File: ...
  <BLANKLINE>
  >>> print(color_relief.band_color)
  {0: 'red', 1: 'green', 2: 'blue', 3: 'alpha'}
  

• The result color relief dataset will have 4 bands red, green, blue, and alpha. with values from 0 to 255.
• To plot the color relief dataset, you can use the plot method. but you need to provide the the rgb indices with the alpha index as the fourth index, otherwise the alpha band will be missing.

  >>> fig, ax = color_relief.plot(rgb=[0, 1, 2, 3]) # doctest: +SKIP
  

  

See Also

Dataset.hill_shade: create a hill-shade for a band in the Dataset.

Source code in src/digitalrivers/terrain.py

def color_relief(
    self,
    band: int = 0,
    path: str | None = None,
    color_table: DataFrame | None = None,
    **kwargs,
) -> "Dataset":
    """Create a color relief for a band in the Dataset.

    A color relief raster is a raster image where each pixel's value is mapped to a specific color based on a
    predefined color palette or color table.

    Args:
        band: int, default is 0.
            band index.
        path: str, default is None.
            path to save the color relief raster.
        color_table: DataFrame, default is None.
            DataFrame with columns: band, values, color
                ```text
                  values    color
                0      1  #709959
                1      2  #F2EEA2
                2      3  #F2CE85
                3      1  #C28C7C
                4      2  #D6C19C
                5      3  #D6C19C
                ```
            or DataFrame with columns: values, red, green, blue, alpha, (the alpha column is optional)
                ```text
                  values    red  green   blue  alpha
                0      1    112    153     89    255
                1      2    242    238    162    255
                2      3    242    206    133    255
                3      1    194    140    124    255
                4      2    214    193    156    255
                5      3    214    193    156    255
                ```
    Returns:
        Dataset:
            Dataset with the color relief with four bands read, green, blue, and alpha.

    Examples:
        - First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.
            ```python
            >>> import numpy as np
            >>> arr = np.random.randint(0, 15, size=(10, 10))
            >>> dataset = Dataset.create_from_array(arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326)

            ```
        - Now let's create the color table using hex colors.
            ```python
            >>> import pandas as pd
            >>> color_hex = ["#709959", "#F2EEA2", "#F2CE85", "#C28C7C", "#D6C19C"]
            >>> values = [1, 3, 5, 7, 9]
            >>> df = pd.DataFrame(columns=["values", "color"])
            >>> df.loc[:, "values"] = values
            >>> df.loc[:, "color"] = color_hex

            ```
        - Now let's create the color relief for the dataset using the color table `DataFrame`.
            ```python
            >>> color_relief = Terrain(dataset.raster).color_relief(band=0, color_table=df)
            >>> print(color_relief) # doctest: +SKIP
            <BLANKLINE>
                        Cell size: 0.05
                        Dimension: 10 * 10
                        EPSG: 4326
                        Number of Bands: 4
                        Band names: ['Band_1', 'Band_2', 'Band_3', 'Band_4']
                        Mask: None
                        Data type: byte
                        projection: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
                        Metadata: {}
                        File: ...
            <BLANKLINE>
            >>> print(color_relief.band_color)
            {0: 'red', 1: 'green', 2: 'blue', 3: 'alpha'}

            ```
        - The result color relief dataset will have 4 bands red, green, blue, and alpha. with values from 0 to 255.
        - To plot the color relief dataset, you can use the `plot` method. but you need to provide the the rgb indices
            with the alpha index as the fourth index, otherwise the alpha band will be missing.
            ```python
            >>> fig, ax = color_relief.plot(rgb=[0, 1, 2, 3]) # doctest: +SKIP

            ```
        ![color-relief](./../_images/dataset/color-relief.png)

    See Also:
        Dataset.hill_shade: create a hill-shade for a band in the Dataset.
    """
    if path is None:
        driver = "MEM"
        path = ""
    else:
        driver = "GTiff"
    color_df = self.analysis._process_color_table(color_table)

    with tempfile.TemporaryDirectory() as temp_dir:
        color_table_path = os.path.join(temp_dir, f"{uuid.uuid1()}.txt")
        color_df.to_csv(color_table_path, index=False, header=False)

        options = gdal.DEMProcessingOptions(
            band=band + 1,
            format=driver,
            colorFilename=color_table_path,
            addAlpha=True,
            creationOptions=CREATION_OPTIONS,
            **kwargs,
        )
        dst = gdal.DEMProcessing(path, self.raster, "color-relief", options=options)

    color_relief = Dataset(dst, access="write")
    color_relief.band_color = {0: "red", 1: "green", 2: "blue", 3: "alpha"}
    return color_relief

hill_shade(band=0, azimuth=315, altitude=45, vertical_exaggeration=1, scale=1, path=None, weights=None, **kwargs)

Create hill-shade.

Hillshade is a technique used in digital elevation modeling (DEM) to create a grayscale representation of a terrain's surface that simulates the effect of sunlight falling across the landscape. This technique helps to visualize the shape and features of the terrain by highlighting the variations in elevation and the slope of the surface.

Hillshade calculates the illumination of each pixel based on the slope (gradient) and aspect (direction) of the terrain surface relative to a specified light source.

The main parameters influencing the hillshade effect are: - Light source direction (Azimuth): the azimuth angle of the light source, which is the angle between the light source - Light source elevation (altitude): the source of light elevation, it is measured in degrees from the horizon. - Vertical exaggeration (Z-factor): the vertical exaggeration is used to emphasize the vertical features of the terrain.

Notes

if the hill_shade parameters are given as lists then the hill shade will be calculated for each set of parameter and then the average will be returned.

Parameters:

Name	Type	Description	Default
band	int	int band index.	0
azimuth	int | float | list[int]	int | float | list[int] The source of light direction, it is measured clockwise from the north. zero means from north to south. 45 degrees means from the northeast to the southwest.	315
altitude	int | float | list[int]	int | float | list[int] The source of light elevation, it is measured in degrees from the horizon. zero means from the horizon. 90 degrees means from the zenith. the overall image gets brighter as the light source gets closer to the zenith. The brightest slopes/DEM features will be perpendicular to the light source, and the darkest will be angled 90˚ or more away.	45
vertical_exaggeration	int | float | list[int]	int | float | list[int] Vertical exaggeration, the vertical exaggeration It is used to emphasize the vertical features of the terrain.	1
scale	int | float | list[int]	int | float | list[int] the scale is the ratio of vertical units to horizontal. If the horizontal unit of the source DEM is degrees (e.g Lat/Long WGS84 projection), you can use scale=111120 if the vertical units are meters (or scale=370400 if they are in feet).	1
path	str | None	str, optional, default is None path to save the hill-shade raster.	None
weights	list[int] | None	list[int], default is None. list of weights to combine the hill-shades if the other parameters are given as lists, an average hill shade will be calculated based on the weights. if None, the weights will be equal.	None
**kwargs		multi_directional: bool if True, the hill shade will be calculated for multiple azimuth values [225, 270, 315, 360] each with a altitude of 30 degrees, and then the average will be returned. with multi_directional = True any given azimuth will be ignored. For more details visit: https://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf combined: bool combined shading, a combination of slope and oblique shading. igor: bool shading which tries to minimize effects on other map features beneath. with igor=True the altitude will be calculated ignored. For more details visit: https://maperitive.net/docs/Commands/GenerateReliefImageIgor.html	{}

Returns:

Name	Type	Description
Dataset	'Dataset'	8-bit Dataset with the hill-shade created.

Examples:

• First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.

  >>> import numpy as np
  >>> arr = np.random.randint(0, 15, size=(100, 100))
  >>> dataset = Dataset.create_from_array(arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326)
  
  >>> hill_shade = Terrain(dataset.raster).hill_shade(
  ...     band=0, altitude=45, azimuth=315, vertical_exaggeration=1, scale=1
  ... )
  
  >>> print(hill_shade.dtype) # doctest: +SKIP
  ['byte']
  >>> hill_shade.plot() # doctest: +SKIP
  

  

  >>> hill_shade.stats() # doctest: +SKIP
          min    max       mean        std
  Band_1  1.0  223.0  58.880951  71.079056
  

• You can also provide the function with a list os values for each parameter, then the functions will calculate the hill shade for each set of parameters and then the average will be returned.

  >>> hill_shade = Terrain(dataset.raster).hill_shade(
  ...     band=0, azimuth=[315, 45], altitude=[45, 45], vertical_exaggeration=[1, 1], scale=[1, 1]
  ... )
  
  >>> hill_shade.plot() # doctest: +SKIP
  

  

See Also

Dataset.color_relief: create a color relief for a band in the Dataset. Dataset.slope: create a slope for a band in the Dataset.

Source code in src/digitalrivers/terrain.py

def hill_shade(
    self,
    band: int = 0,
    azimuth: int | float | list[int] = 315,
    altitude: int | float | list[int] = 45,
    vertical_exaggeration: int | float | list[int] = 1,
    scale: int | float | list[int] = 1,
    path: str | None = None,
    weights: list[int] | None = None,
    **kwargs,
) -> "Dataset":
    """Create hill-shade.

    Hillshade is a technique used in digital elevation modeling (DEM) to create a grayscale representation of a
    terrain's surface that simulates the effect of sunlight falling across the landscape.
    This technique helps to visualize the shape and features of the terrain by highlighting the variations in
    elevation and the slope of the surface.

    Hillshade calculates the illumination of each pixel based on the slope (gradient) and aspect (direction) of the
    terrain surface relative to a specified light source.

    The main parameters influencing the hillshade effect are:
    - Light source direction (Azimuth): the azimuth angle of the light source, which is the angle between the light
        source
    - Light source elevation (altitude): the source of light elevation, it is measured in degrees from the horizon.
    - Vertical exaggeration (Z-factor): the vertical exaggeration is used to emphasize the vertical features of the
        terrain.

    Notes:
        if the `hill_shade` parameters are given as lists then the hill shade will be calculated for each set
        of parameter and then the average will be returned.

    Args:
        band: int
            band index.
        azimuth: int | float | list[int]
            The source of light direction, it is measured clockwise from the north. zero means from north to south.
            45 degrees means from the northeast to the southwest.
        altitude: int | float | list[int]
            The source of light elevation, it is measured in degrees from the horizon. zero means from the horizon.
            90 degrees means from the zenith.
            the overall image gets brighter as the light source gets closer to the zenith. The brightest slopes/DEM
            features will be perpendicular to the light source, and the darkest will be angled 90˚ or more away.
        vertical_exaggeration: int | float | list[int]
            Vertical exaggeration, the vertical exaggeration It is used to emphasize the
            vertical features of the terrain.
        scale: int | float | list[int]
            the scale is the ratio of vertical units to horizontal. If the horizontal unit of the source DEM is
            degrees (e.g Lat/Long WGS84 projection), you can use scale=111120 if the vertical units are meters
            (or scale=370400 if they are in feet).
        path: str, optional, default is None
            path to save the hill-shade raster.
        weights: list[int], default is None.
            list of weights to combine the hill-shades if the other parameters are given as lists, an average hill
            shade will be calculated based on the weights. if None, the weights will be equal.
        **kwargs:
            multi_directional: bool
                if True, the hill shade will be calculated for multiple azimuth values [225, 270, 315, 360] each with a
                altitude of 30 degrees, and then the average will be returned. with multi_directional = True any given
                azimuth will be ignored.
                For more details visit: https://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf
            combined: bool
                combined shading, a combination of slope and oblique shading.
            igor: bool
                shading which tries to minimize effects on other map features beneath. with `igor=True` the altitude
                will be calculated ignored.
                For more details visit: https://maperitive.net/docs/Commands/GenerateReliefImageIgor.html

    Returns:
        Dataset: 8-bit
            Dataset with the hill-shade created.

    Examples:
        - First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.
            ```python
            >>> import numpy as np
            >>> arr = np.random.randint(0, 15, size=(100, 100))
            >>> dataset = Dataset.create_from_array(arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326)

            >>> hill_shade = Terrain(dataset.raster).hill_shade(
            ...     band=0, altitude=45, azimuth=315, vertical_exaggeration=1, scale=1
            ... )

            >>> print(hill_shade.dtype) # doctest: +SKIP
            ['byte']
            >>> hill_shade.plot() # doctest: +SKIP

            ```
            ![hill-shade](./../_images/dataset/hill-shade.png)
            ```python
            >>> hill_shade.stats() # doctest: +SKIP
                    min    max       mean        std
            Band_1  1.0  223.0  58.880951  71.079056

            ```
        - You can also provide the function with a list os values for each parameter, then the functions will
            calculate the hill shade for each set of parameters and then the average will be returned.
            ```python
            >>> hill_shade = Terrain(dataset.raster).hill_shade(
            ...     band=0, azimuth=[315, 45], altitude=[45, 45], vertical_exaggeration=[1, 1], scale=[1, 1]
            ... )

            >>> hill_shade.plot() # doctest: +SKIP

            ```
            ![hill-shade-multi](./../_images/dataset/hill-shade-multi.png)

    See Also:
        Dataset.color_relief: create a color relief for a band in the Dataset.
        Dataset.slope: create a slope for a band in the Dataset.
    """
    if "multi_directional" in kwargs:
        if not isinstance(kwargs["multi_directional"], bool):
            raise ValueError("The multi_directional parameter must be a boolean.")
        if kwargs["multi_directional"]:
            multi_directional = True
            azimuth = None
            # altitude, vertical_exaggeration, scale = None, None, None,
        else:
            multi_directional = False

        kwargs.pop("multi_directional")
        kwargs["multiDirectional"] = multi_directional
    if "igor" in kwargs:
        if not isinstance(kwargs["igor"], bool):
            raise ValueError("The igor parameter must be a boolean.")
        if kwargs["igor"]:
            altitude = None

    # if not (
    #     type(azimuth)
    #     is type(altitude)
    #     is type(vertical_exaggeration)
    #     is type(scale)
    # ):
    #     raise ValueError(
    #         f"The azimuth, altitude, vertical_exaggeration, and scale parameter must be of the same type. Given"
    #         f" azimuth: {type(azimuth)}, altitude: {type(altitude)}, vertical_exaggeration: {type(vertical_exaggeration)},"
    #         f"scale: {type(scale)}"
    #     )

    if path is None:
        driver = "MEM"
        path = ""
    else:
        driver = "GTiff"

    wrap = lambda v: v if isinstance(v, list) else [v]
    azimuth, altitude, vertical_exaggeration, scale = (
        wrap(azimuth),
        wrap(altitude),
        wrap(vertical_exaggeration),
        wrap(scale),
    )
    if not (
        len(azimuth) == len(altitude) == len(vertical_exaggeration) == len(scale)
    ):
        raise ValueError("All list parameters must have the same length.")

    # get the hill shade for all the parameters
    hill_shades: list[gdal.Dataset] = []
    for az, alt, ver_ex, scale_1 in zip(
        azimuth, altitude, vertical_exaggeration, scale
    ):
        dst = self._create_hill_shade(
            band, driver, az, alt, ver_ex, scale_1, path, **kwargs
        )
        hill_shades.append(dst)

    if len(hill_shades) > 1:
        if weights is None:
            weights = np.ones(len(azimuth))
        weights = np.array(weights) / np.sum(weights)
        hill_shades_arr: list[np.ndarray] = [
            hill_shade.ReadAsArray() for hill_shade in hill_shades
        ]
        combined_hillshade = np.average(hill_shades_arr, axis=0, weights=weights)
        combined_hillshade = np.clip(combined_hillshade, 0, 255).astype(np.uint8)
        hill_shade = Dataset.dataset_like(
            Dataset(hill_shades[0]), combined_hillshade
        )
    else:
        hill_shade = Dataset(hill_shades[0], access="write")

    hill_shade.band_color = {0: "gray_index"}

    return hill_shade

slope(band=0, scale=1, slope_format='degree', path=None, algorithm=None, creation_options=None, **kwargs)

Compute the slope of the terrain surface.

Uses GDAL DEMProcessing to calculate the slope (rate of elevation change) for every cell.

Parameters:

Name	Type	Description	Default
band	int	Zero-based band index. Defaults to 0.	0
scale	int | float | list[int]	Ratio of vertical to horizontal units. Use 111120 when the horizontal CRS is in degrees and vertical units are metres. Defaults to 1.	1
slope_format	str	Output format — "degree" (default) or "percent".	'degree'
algorithm	str | None	Slope algorithm. One of "Horn", "ZevenbergenThorne", or None (GDAL default). Zevenbergen-Thorne suits smooth landscapes; Horn performs better on rough terrain.	None
path	str | None	If given, write the result to this GeoTIFF path. Otherwise the raster is created in memory.	None
creation_options	list[str] | None	GDAL creation options. Defaults to ['COMPRESS=DEFLATE', 'PREDICTOR=2'].	None
**kwargs		Forwarded to gdal.DEMProcessingOptions.	{}

Returns:

Name	Type	Description
Dataset	'Dataset'	Single-band float32 raster with slope values. No-data value is -9999.0.

Examples:

• First create a one band dataset, consisting of 10 columns and 10 rows, with random values between 0 and 15.

  >>> import numpy as np
  >>> arr = np.random.randint(0, 15, size=(10, 10))
  >>> dataset = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326
  ... )
  

• Now let's create the slope for the dataset.

  >>> slope = Terrain(dataset.raster).slope()
  >>> fig, ax = slope.plot() # doctest: +SKIP
  

  

See Also

Terrain.hill_shade: Create a hill-shade for a band in the Dataset. Terrain.color_relief: Create a color relief for a band in the Dataset.

Source code in src/digitalrivers/terrain.py

def slope(
    self,
    band: int = 0,
    scale: int | float | list[int] = 1,
    slope_format: str = "degree",
    path: str | None = None,
    algorithm: str | None = None,
    creation_options: list[str] | None = None,
    **kwargs,
) -> "Dataset":
    """Compute the slope of the terrain surface.

    Uses GDAL `DEMProcessing` to calculate the slope (rate of
    elevation change) for every cell.

    Args:
        band: Zero-based band index. Defaults to 0.
        scale: Ratio of vertical to horizontal units.  Use
            `111120` when the horizontal CRS is in degrees and
            vertical units are metres.  Defaults to 1.
        slope_format: Output format — `"degree"` (default) or
            `"percent"`.
        algorithm: Slope algorithm.  One of `"Horn"`,
            `"ZevenbergenThorne"`, or `None` (GDAL default).
            Zevenbergen-Thorne suits smooth landscapes; Horn
            performs better on rough terrain.
        path: If given, write the result to this GeoTIFF path.
            Otherwise the raster is created in memory.
        creation_options: GDAL creation options.  Defaults to
            `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
        **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

    Returns:
        Dataset: Single-band `float32` raster with slope values.
            No-data value is `-9999.0`.

    Examples:
        - First create a one band dataset, consisting of 10 columns
            and 10 rows, with random values between 0 and 15.
            ```python
            >>> import numpy as np
            >>> arr = np.random.randint(0, 15, size=(10, 10))
            >>> dataset = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326
            ... )

            ```
        - Now let's create the slope for the dataset.
            ```python
            >>> slope = Terrain(dataset.raster).slope()
            >>> fig, ax = slope.plot() # doctest: +SKIP

            ```
            ![slope](./../_images/dataset/slope.png)

    See Also:
        Terrain.hill_shade: Create a hill-shade for a band in the
            Dataset.
        Terrain.color_relief: Create a color relief for a band in
            the Dataset.
    """
    if path is None:
        driver = "MEM"
        path = ""
    else:
        driver = "GTiff"

    if creation_options is None:
        creation_options = CREATION_OPTIONS.copy()

    options = gdal.DEMProcessingOptions(
        band=band + 1,
        format=driver,
        alg=algorithm,
        slopeFormat=slope_format,
        scale=scale,
        creationOptions=creation_options,
        **kwargs,
    )
    dst = gdal.DEMProcessing(path, self.raster, "slope", options=options)
    src = Dataset(dst, access="write")

    return src

aspect(band=0, scale=1, vertical_exaggeration=1, zero_flat_surface=False, algorithm=None, path=None, creation_options=None, **kwargs)

Compute the aspect (slope direction) of the terrain surface.

Uses GDAL DEMProcessing to calculate the compass direction of the steepest downhill slope for every cell. Values range from 0° (north) clockwise to 360°.

Parameters:

Name	Type	Description	Default
band	int	Zero-based band index. Defaults to 0.	0
scale	int | float | list[int]	Accepted for signature symmetry with slope but not used — GDAL's aspect mode is scale-invariant (it reports a direction) and rejects the -s option. Defaults to 1.	1
vertical_exaggeration	int | float | list[int]	Accepted for signature symmetry but not used for the same reason (aspect rejects the -z option). Defaults to 1.	1
zero_flat_surface	bool	If True flat areas get an aspect of 0°. If False (default) flat areas receive the no-data value.	False
algorithm	str | None	Aspect algorithm. One of "Horn", "ZevenbergenThorne", or None (GDAL default).	None
path	str | None	If given, write the result to this GeoTIFF path. Otherwise the raster is created in memory.	None
creation_options	list[str] | None	GDAL creation options. Defaults to ['COMPRESS=DEFLATE', 'PREDICTOR=2'].	None
**kwargs		Forwarded to gdal.DEMProcessingOptions.	{}

Returns:

Name	Type	Description
Dataset	'Dataset'	Single-band float32 raster with aspect values in degrees (0–360). No-data value is -9999.0.

Examples:

• Create a small raster and compute its aspect.

  >>> import numpy as np
  >>> arr = np.random.randint(0, 15, size=(10, 10))
  >>> dataset = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326
  ... )
  

• Compute the aspect raster.

  >>> aspect = Terrain(dataset.raster).aspect()
  >>> fig, ax = aspect.plot() # doctest: +SKIP
  

  

See Also

Terrain.hill_shade: Create a hill-shade for a band in the Dataset. Terrain.slope: Compute the slope of the terrain surface.

Source code in src/digitalrivers/terrain.py

def aspect(
    self,
    band: int = 0,
    scale: int | float | list[int] = 1,
    vertical_exaggeration: int | float | list[int] = 1,
    zero_flat_surface: bool = False,
    algorithm: str | None = None,
    path: str | None = None,
    creation_options: list[str] | None = None,
    **kwargs,
) -> "Dataset":
    """Compute the aspect (slope direction) of the terrain surface.

    Uses GDAL `DEMProcessing` to calculate the compass direction
    of the steepest downhill slope for every cell.  Values range
    from 0° (north) clockwise to 360°.

    Args:
        band: Zero-based band index. Defaults to 0.
        scale: Accepted for signature symmetry with `slope` but **not
            used** — GDAL's aspect mode is scale-invariant (it reports a
            direction) and rejects the `-s` option. Defaults to 1.
        vertical_exaggeration: Accepted for signature symmetry but **not
            used** for the same reason (aspect rejects the `-z` option).
            Defaults to 1.
        zero_flat_surface: If `True` flat areas get an aspect of
            0°.  If `False` (default) flat areas receive the
            no-data value.
        algorithm: Aspect algorithm.  One of `"Horn"`,
            `"ZevenbergenThorne"`, or `None` (GDAL default).
        path: If given, write the result to this GeoTIFF path.
            Otherwise the raster is created in memory.
        creation_options: GDAL creation options.  Defaults to
            `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
        **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

    Returns:
        Dataset: Single-band `float32` raster with aspect values
            in degrees (0–360).  No-data value is `-9999.0`.

    Examples:
        - Create a small raster and compute its aspect.
            ```python
            >>> import numpy as np
            >>> arr = np.random.randint(0, 15, size=(10, 10))
            >>> dataset = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=0.05, epsg=4326
            ... )

            ```
        - Compute the aspect raster.
            ```python
            >>> aspect = Terrain(dataset.raster).aspect()
            >>> fig, ax = aspect.plot() # doctest: +SKIP

            ```
            ![aspect](./../_images/dataset/aspect.png)

    See Also:
        Terrain.hill_shade: Create a hill-shade for a band in the
            Dataset.
        Terrain.slope: Compute the slope of the terrain surface.
    """
    if path is None:
        driver = "MEM"
        path = ""
    else:
        driver = "GTiff"

    if creation_options is None:
        creation_options = CREATION_OPTIONS.copy()

    # `scale` / `vertical_exaggeration` are NOT forwarded: GDAL's aspect mode
    # does not accept the slope/hillshade `-s` / `-z` options (aspect is a
    # scale-invariant direction). Passing them is silently tolerated only
    # when `zeroForFlat` is unset, and otherwise raises "Argument(s) are not
    # valid with any processing mode" — so they are dropped here.
    del scale, vertical_exaggeration
    options = gdal.DEMProcessingOptions(
        band=band + 1,
        format=driver,
        alg=algorithm,
        zeroForFlat=zero_flat_surface,
        creationOptions=creation_options,
        **kwargs,
    )
    dst = gdal.DEMProcessing(path, self.raster, "aspect", options=options)
    src = Dataset(dst, access="write")

    return src

roughness(band=0, path=None, compute_edges=False, creation_options=None, **kwargs)

Compute terrain roughness — the largest elevation difference in a 3x3 window.

Roughness (Wilson et al., 2007) is the maximum absolute difference between a cell and its eight neighbours. It is the simplest ruggedness measure and reacts strongly to local relief: flat surfaces score 0, cliffs and noisy LiDAR returns score high.

Parameters:

Name	Type	Description	Default
band	int	Zero-based band index. Defaults to 0.	0
path	str | None	If given, write the result to this GeoTIFF path. Otherwise the raster is created in memory.	None
compute_edges	bool	If True, compute values for the edge cells from the available partial window instead of leaving them no-data. Defaults to False.	False
creation_options	list[str] | None	GDAL creation options. Defaults to ['COMPRESS=DEFLATE', 'PREDICTOR=2'].	None
**kwargs		Forwarded to gdal.DEMProcessingOptions.	{}

Returns:

Name	Type	Description
Dataset	'Dataset'	Single-band float32 raster of roughness values in the DEM's vertical units. No-data value is -9999.0.

Examples:

• Compute roughness for a small elevation raster.

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> arr = np.array(
  ...     [[10, 11, 12, 40], [10, 9, 8, 7],
  ...      [5, 6, 30, 6], [4, 3, 2, 1]],
  ...     dtype=np.float32,
  ... )
  >>> ds = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=1.0,
  ...     epsg=32636, no_data_value=-9999.0,
  ... )
  >>> roughness = Terrain(ds.raster).roughness()
  >>> roughness.read_array().shape
  (4, 4)
  

See Also

Terrain.tpi: Topographic Position Index. Terrain.tri: Terrain Ruggedness Index.

Source code in src/digitalrivers/terrain.py

def roughness(
    self,
    band: int = 0,
    path: str | None = None,
    compute_edges: bool = False,
    creation_options: list[str] | None = None,
    **kwargs,
) -> "Dataset":
    """Compute terrain roughness — the largest elevation difference in a 3x3 window.

    Roughness (Wilson et al., 2007) is the maximum absolute difference
    between a cell and its eight neighbours. It is the simplest
    ruggedness measure and reacts strongly to local relief: flat
    surfaces score `0`, cliffs and noisy LiDAR returns score high.

    Args:
        band: Zero-based band index. Defaults to 0.
        path: If given, write the result to this GeoTIFF path.
            Otherwise the raster is created in memory.
        compute_edges: If `True`, compute values for the edge cells
            from the available partial window instead of leaving them
            no-data. Defaults to `False`.
        creation_options: GDAL creation options. Defaults to
            `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
        **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

    Returns:
        Dataset: Single-band `float32` raster of roughness values in
            the DEM's vertical units. No-data value is `-9999.0`.

    Examples:
        - Compute roughness for a small elevation raster.
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> arr = np.array(
            ...     [[10, 11, 12, 40], [10, 9, 8, 7],
            ...      [5, 6, 30, 6], [4, 3, 2, 1]],
            ...     dtype=np.float32,
            ... )
            >>> ds = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=1.0,
            ...     epsg=32636, no_data_value=-9999.0,
            ... )
            >>> roughness = Terrain(ds.raster).roughness()
            >>> roughness.read_array().shape
            (4, 4)

            ```

    See Also:
        Terrain.tpi: Topographic Position Index.
        Terrain.tri: Terrain Ruggedness Index.
    """
    return self._ruggedness(
        "Roughness", band, path, compute_edges, creation_options, **kwargs
    )

tpi(band=0, path=None, compute_edges=False, creation_options=None, **kwargs)

Compute the Topographic Position Index (TPI).

TPI (Weiss, 2001) is each cell's elevation minus the mean elevation of its eight neighbours. Positive values mark local highs (ridges, peaks), negative values mark local lows (valleys, channels), and values near zero mark flat areas or constant slopes. It is widely used for landform classification.

Note

This is the GDAL formulation — the focal mean is taken over the eight neighbours excluding the centre cell, on a fixed 3x3 window. DEM.tpi is a native alternative whose focal mean includes the centre cell and accepts an arbitrary window size, so the two return slightly different values.

Parameters:

Name	Type	Description	Default
band	int	Zero-based band index. Defaults to 0.	0
path	str | None	If given, write the result to this GeoTIFF path. Otherwise the raster is created in memory.	None
compute_edges	bool	If True, compute values for the edge cells from the available partial window instead of leaving them no-data. Defaults to False.	False
creation_options	list[str] | None	GDAL creation options. Defaults to ['COMPRESS=DEFLATE', 'PREDICTOR=2'].	None
**kwargs		Forwarded to gdal.DEMProcessingOptions.	{}

Returns:

Name	Type	Description
Dataset	'Dataset'	Single-band float32 raster of TPI values (signed, in the DEM's vertical units). No-data value is -9999.0.

Examples:

• Compute TPI for a small elevation raster.

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> arr = np.array(
  ...     [[10, 11, 12, 40], [10, 9, 8, 7],
  ...      [5, 6, 30, 6], [4, 3, 2, 1]],
  ...     dtype=np.float32,
  ... )
  >>> ds = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=1.0,
  ...     epsg=32636, no_data_value=-9999.0,
  ... )
  >>> tpi = Terrain(ds.raster).tpi()
  >>> tpi.read_array().shape
  (4, 4)
  

See Also

Terrain.roughness: Maximum 3x3 elevation difference. Terrain.tri: Terrain Ruggedness Index. DEM.tpi: Native, window-configurable TPI (includes the centre cell in the focal mean).

Source code in src/digitalrivers/terrain.py

def tpi(
    self,
    band: int = 0,
    path: str | None = None,
    compute_edges: bool = False,
    creation_options: list[str] | None = None,
    **kwargs,
) -> "Dataset":
    """Compute the Topographic Position Index (TPI).

    TPI (Weiss, 2001) is each cell's elevation minus the mean
    elevation of its eight neighbours. Positive values mark local
    highs (ridges, peaks), negative values mark local lows (valleys,
    channels), and values near zero mark flat areas or constant
    slopes. It is widely used for landform classification.

    Note:
        This is the GDAL formulation — the focal mean is taken over
        the eight neighbours **excluding** the centre cell, on a fixed
        3x3 window. `DEM.tpi` is a native alternative whose focal mean
        **includes** the centre cell and accepts an arbitrary
        `window` size, so the two return slightly different values.

    Args:
        band: Zero-based band index. Defaults to 0.
        path: If given, write the result to this GeoTIFF path.
            Otherwise the raster is created in memory.
        compute_edges: If `True`, compute values for the edge cells
            from the available partial window instead of leaving them
            no-data. Defaults to `False`.
        creation_options: GDAL creation options. Defaults to
            `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
        **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

    Returns:
        Dataset: Single-band `float32` raster of TPI values (signed,
            in the DEM's vertical units). No-data value is `-9999.0`.

    Examples:
        - Compute TPI for a small elevation raster.
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> arr = np.array(
            ...     [[10, 11, 12, 40], [10, 9, 8, 7],
            ...      [5, 6, 30, 6], [4, 3, 2, 1]],
            ...     dtype=np.float32,
            ... )
            >>> ds = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=1.0,
            ...     epsg=32636, no_data_value=-9999.0,
            ... )
            >>> tpi = Terrain(ds.raster).tpi()
            >>> tpi.read_array().shape
            (4, 4)

            ```

    See Also:
        Terrain.roughness: Maximum 3x3 elevation difference.
        Terrain.tri: Terrain Ruggedness Index.
        DEM.tpi: Native, window-configurable TPI (includes the centre
            cell in the focal mean).
    """
    return self._ruggedness(
        "TPI", band, path, compute_edges, creation_options, **kwargs
    )

tri(band=0, algorithm=None, path=None, compute_edges=False, creation_options=None, **kwargs)

Compute the Terrain Ruggedness Index (TRI).

TRI is the mean absolute difference between a cell and its eight neighbours. Two formulations are available via algorithm:

• "Riley" (Riley et al., 1999) — square-root of the summed squared differences; the original TRI.
• "Wilson" (Wilson et al., 2007) — the mean absolute difference; better suited to bathymetric / continuous data.

Note

With algorithm=None GDAL uses the Riley root-sum-square form. The native DEM.ruggedness computes the Wilson mean-absolute-difference form, so it corresponds to tri(algorithm="Wilson") (on a 3x3 window) rather than the default here.

Parameters:

Name	Type	Description	Default
band	int	Zero-based band index. Defaults to 0.	0
algorithm	str | None	TRI formulation — "Riley", "Wilson", or None (GDAL default). Defaults to None.	None
path	str | None	If given, write the result to this GeoTIFF path. Otherwise the raster is created in memory.	None
compute_edges	bool	If True, compute values for the edge cells from the available partial window instead of leaving them no-data. Defaults to False.	False
creation_options	list[str] | None	GDAL creation options. Defaults to ['COMPRESS=DEFLATE', 'PREDICTOR=2'].	None
**kwargs		Forwarded to gdal.DEMProcessingOptions.	{}

Returns:

Name	Type	Description
Dataset	'Dataset'	Single-band float32 raster of TRI values in the DEM's vertical units. No-data value is -9999.0.

Examples:

• Compute TRI for a small elevation raster.

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> arr = np.array(
  ...     [[10, 11, 12, 40], [10, 9, 8, 7],
  ...      [5, 6, 30, 6], [4, 3, 2, 1]],
  ...     dtype=np.float32,
  ... )
  >>> ds = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=1.0,
  ...     epsg=32636, no_data_value=-9999.0,
  ... )
  >>> tri = Terrain(ds.raster).tri()
  >>> tri.read_array().shape
  (4, 4)
  

See Also

Terrain.roughness: Maximum 3x3 elevation difference. Terrain.tpi: Topographic Position Index. DEM.ruggedness: Native, window-configurable Wilson-form TRI.

Source code in src/digitalrivers/terrain.py

def tri(
    self,
    band: int = 0,
    algorithm: str | None = None,
    path: str | None = None,
    compute_edges: bool = False,
    creation_options: list[str] | None = None,
    **kwargs,
) -> "Dataset":
    """Compute the Terrain Ruggedness Index (TRI).

    TRI is the mean absolute difference between a cell and its eight
    neighbours. Two formulations are available via `algorithm`:

    * `"Riley"` (Riley et al., 1999) — square-root of the summed
      squared differences; the original TRI.
    * `"Wilson"` (Wilson et al., 2007) — the mean absolute
      difference; better suited to bathymetric / continuous data.

    Note:
        With `algorithm=None` GDAL uses the Riley root-sum-square
        form. The native `DEM.ruggedness` computes the Wilson
        mean-absolute-difference form, so it corresponds to
        `tri(algorithm="Wilson")` (on a 3x3 window) rather than the
        default here.

    Args:
        band: Zero-based band index. Defaults to 0.
        algorithm: TRI formulation — `"Riley"`, `"Wilson"`, or
            `None` (GDAL default). Defaults to `None`.
        path: If given, write the result to this GeoTIFF path.
            Otherwise the raster is created in memory.
        compute_edges: If `True`, compute values for the edge cells
            from the available partial window instead of leaving them
            no-data. Defaults to `False`.
        creation_options: GDAL creation options. Defaults to
            `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.
        **kwargs: Forwarded to `gdal.DEMProcessingOptions`.

    Returns:
        Dataset: Single-band `float32` raster of TRI values in the
            DEM's vertical units. No-data value is `-9999.0`.

    Examples:
        - Compute TRI for a small elevation raster.
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> arr = np.array(
            ...     [[10, 11, 12, 40], [10, 9, 8, 7],
            ...      [5, 6, 30, 6], [4, 3, 2, 1]],
            ...     dtype=np.float32,
            ... )
            >>> ds = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=1.0,
            ...     epsg=32636, no_data_value=-9999.0,
            ... )
            >>> tri = Terrain(ds.raster).tri()
            >>> tri.read_array().shape
            (4, 4)

            ```

    See Also:
        Terrain.roughness: Maximum 3x3 elevation difference.
        Terrain.tpi: Topographic Position Index.
        DEM.ruggedness: Native, window-configurable Wilson-form TRI.
    """
    if algorithm is not None:
        kwargs["alg"] = algorithm
    return self._ruggedness(
        "TRI", band, path, compute_edges, creation_options, **kwargs
    )

viewshed(observer_x, observer_y, band=0, observer_height=1.75, target_height=0.0, max_distance=0.0, mode='max', visible_value=255.0, invisible_value=0.0, out_of_range_value=0.0, no_data_value=-1.0, curvature_coefficient=0.85714, path=None, creation_options=None)

Compute the viewshed (line-of-sight visibility) from an observer point.

Wraps GDAL ViewshedGenerate to flag, for every cell, whether it is visible from an observer standing at (observer_x, observer_y), accounting for the intervening terrain. The observer and target heights are added above the DEM surface, and Earth curvature / atmospheric refraction can be modelled via curvature_coefficient.

Parameters:

Name	Type	Description	Default
observer_x	float	Observer X coordinate, in the DEM's CRS.	required
observer_y	float	Observer Y coordinate, in the DEM's CRS.	required
band	int	Zero-based band index of the elevation band. Defaults to 0.	0
observer_height	float	Observer height above the DEM surface, in the DEM's vertical units. Defaults to 1.75 (roughly eye level).	1.75
target_height	float	Target height above the DEM surface that must be visible. Defaults to 0.0 (ground level).	0.0
max_distance	float	Maximum line-of-sight distance in CRS units. 0.0 (default) means unlimited.	0.0
mode	str	Cell-evaluation method — "max" (default), "min", "edge", or "diagonal" — mapping to the GDAL GVM_Max / GVM_Min / GVM_Edge / GVM_Diagonal viewshed modes.	'max'
visible_value	float	Output value written to visible cells. Defaults to 255.0.	255.0
invisible_value	float	Output value written to hidden cells. Defaults to 0.0.	0.0
out_of_range_value	float	Output value for cells beyond max_distance. Defaults to 0.0.	0.0
no_data_value	float	Output no-data value. Defaults to -1.0.	-1.0
curvature_coefficient	float	Earth-curvature / refraction coefficient. Defaults to 0.85714 (GDAL's standard atmospheric value); use 1.0 to ignore curvature.	0.85714
path	str | None	If given, write the result to this GeoTIFF path. Otherwise the raster is created in memory.	None
creation_options	list[str] | None	GDAL creation options. Defaults to ['COMPRESS=DEFLATE', 'PREDICTOR=2'].	None

Returns:

Name	Type	Description
Dataset	'Dataset'	Single-band raster encoding visibility (visible_value for visible cells, invisible_value otherwise).

Raises:

Type	Description
ValueError	If mode is not one of "max", "min", "edge", "diagonal".

Examples:

• Compute the viewshed from the top-left corner of a small DEM.

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from digitalrivers import Terrain
  >>> arr = np.array(
  ...     [[10, 11, 12, 40], [10, 9, 8, 7],
  ...      [5, 6, 30, 6], [4, 3, 2, 1]],
  ...     dtype=np.float32,
  ... )
  >>> ds = Dataset.create_from_array(
  ...     arr, top_left_corner=(0, 0), cell_size=1.0,
  ...     epsg=32636, no_data_value=-9999.0,
  ... )
  >>> vs = Terrain(ds.raster).viewshed(
  ...     observer_x=0.5, observer_y=-0.5,
  ... )
  >>> vs.read_array().shape
  (4, 4)
  

See Also

Terrain.hill_shade: Shaded-relief visualisation of the surface.

Source code in src/digitalrivers/terrain.py

def viewshed(
    self,
    observer_x: float,
    observer_y: float,
    band: int = 0,
    observer_height: float = 1.75,
    target_height: float = 0.0,
    max_distance: float = 0.0,
    mode: str = "max",
    visible_value: float = 255.0,
    invisible_value: float = 0.0,
    out_of_range_value: float = 0.0,
    no_data_value: float = -1.0,
    curvature_coefficient: float = 0.85714,
    path: str | None = None,
    creation_options: list[str] | None = None,
) -> "Dataset":
    """Compute the viewshed (line-of-sight visibility) from an observer point.

    Wraps GDAL `ViewshedGenerate` to flag, for every cell, whether it
    is visible from an observer standing at `(observer_x, observer_y)`,
    accounting for the intervening terrain. The observer and target
    heights are added above the DEM surface, and Earth curvature /
    atmospheric refraction can be modelled via `curvature_coefficient`.

    Args:
        observer_x: Observer X coordinate, in the DEM's CRS.
        observer_y: Observer Y coordinate, in the DEM's CRS.
        band: Zero-based band index of the elevation band. Defaults
            to 0.
        observer_height: Observer height above the DEM surface, in the
            DEM's vertical units. Defaults to 1.75 (roughly eye level).
        target_height: Target height above the DEM surface that must be
            visible. Defaults to 0.0 (ground level).
        max_distance: Maximum line-of-sight distance in CRS units.
            `0.0` (default) means unlimited.
        mode: Cell-evaluation method — `"max"` (default), `"min"`,
            `"edge"`, or `"diagonal"` — mapping to the GDAL
            `GVM_Max` / `GVM_Min` / `GVM_Edge` / `GVM_Diagonal`
            viewshed modes.
        visible_value: Output value written to visible cells. Defaults
            to 255.0.
        invisible_value: Output value written to hidden cells. Defaults
            to 0.0.
        out_of_range_value: Output value for cells beyond
            `max_distance`. Defaults to 0.0.
        no_data_value: Output no-data value. Defaults to -1.0.
        curvature_coefficient: Earth-curvature / refraction
            coefficient. Defaults to 0.85714 (GDAL's standard
            atmospheric value); use 1.0 to ignore curvature.
        path: If given, write the result to this GeoTIFF path.
            Otherwise the raster is created in memory.
        creation_options: GDAL creation options. Defaults to
            `['COMPRESS=DEFLATE', 'PREDICTOR=2']`.

    Returns:
        Dataset: Single-band raster encoding visibility (`visible_value`
            for visible cells, `invisible_value` otherwise).

    Raises:
        ValueError: If `mode` is not one of `"max"`, `"min"`,
            `"edge"`, `"diagonal"`.

    Examples:
        - Compute the viewshed from the top-left corner of a small DEM.
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from digitalrivers import Terrain
            >>> arr = np.array(
            ...     [[10, 11, 12, 40], [10, 9, 8, 7],
            ...      [5, 6, 30, 6], [4, 3, 2, 1]],
            ...     dtype=np.float32,
            ... )
            >>> ds = Dataset.create_from_array(
            ...     arr, top_left_corner=(0, 0), cell_size=1.0,
            ...     epsg=32636, no_data_value=-9999.0,
            ... )
            >>> vs = Terrain(ds.raster).viewshed(
            ...     observer_x=0.5, observer_y=-0.5,
            ... )
            >>> vs.read_array().shape
            (4, 4)

            ```

    See Also:
        Terrain.hill_shade: Shaded-relief visualisation of the surface.
    """
    modes = {
        "diagonal": gdal.GVM_Diagonal,
        "edge": gdal.GVM_Edge,
        "max": gdal.GVM_Max,
        "min": gdal.GVM_Min,
    }
    if mode not in modes:
        raise ValueError(f"mode must be one of {sorted(modes)}; got {mode!r}")

    if path is None:
        driver = "MEM"
        path = ""
    else:
        driver = "GTiff"

    if creation_options is None:
        creation_options = CREATION_OPTIONS.copy()

    src_band = self.raster.GetRasterBand(band + 1)
    # `gdal.ViewshedGenerate` + the `GVM_*` enums require a recent GDAL; that
    # floor is supplied transitively by the `pyramids-gis >=0.29.0` pin (see
    # the dependency note in pyproject.toml), not a direct GDAL dependency.
    dst = gdal.ViewshedGenerate(
        src_band,
        driver,
        path,
        creation_options,
        observer_x,
        observer_y,
        observer_height,
        target_height,
        visible_value,
        invisible_value,
        out_of_range_value,
        no_data_value,
        curvature_coefficient,
        modes[mode],
        max_distance,
        heightMode=gdal.GVOT_NORMAL,
    )
    return Dataset(dst, access="write")