Deriving Slope, Rugosity & Terrain Metrics #

Operational intent: this topic turns a finished bathymetric surface into the derived layers that downstream science and engineering actually consume — slope for stability and trawlability, rugosity for benthic habitat, bathymetric position index for identifying banks and channels. It sits inside Bathymetric Processing & Terrain Modeling, immediately after surface smoothing and before any habitat or engineering model. The computation is cheap; the interpretation is not, and almost every mistake in this stage is a mistake about scale rather than about arithmetic.

The governing fact is this: a terrain derivative is a property of the grid at least as much as of the seabed. Slope computed on a 1 m grid and slope computed on the same seabed at 50 m are different quantities, not different-quality estimates of one quantity. Neither is the “true” slope, because the seabed has structure at both scales and the question “how steep is it here” has no answer until a neighbourhood size is supplied. A slope raster published without its analysis scale is uninterpretable, and it is published without its scale most of the time.

Reference Configuration #

Two configuration values are not in that table and matter more than anything in it: the cell size and the vertical exaggeration between the horizontal and vertical units. Slope is a ratio; if the grid’s horizontal unit is degrees and its depths are metres, the ratio is dimensionally meaningless and the resulting raster still renders as a plausible-looking slope map. Assert projected units before computing anything.

Slope: the Kernel Choice Matters Less Than the Scale #

Horn’s method — the 3×3 weighted difference used by GDAL, ArcGIS and QGIS — is the de-facto standard and is a fine choice. Zevenbergen–Thorne fits a quadratic and gives slightly different values on rough surfaces. Fitting a plane by least squares over a wider radius gives a scale-controlled estimate. The differences between these methods at fixed scale are a few per cent; the difference between scales is a factor of twenty, as the diagram above shows.

That ordering has a practical consequence: effort spent choosing between kernels is misdirected if the scale has not been chosen deliberately. The scale should come from the physical question. Trawl gear responds to slope over its own footprint, tens of metres. A pipeline free-span analysis responds to slope over the span length. A sediment model responds to slope at the scale of the bedforms it resolves. Each is a different raster and they should be different files, named with their radius.

"""Scale-controlled seafloor slope by least-squares plane fit."""
from __future__ import annotations

import logging

import numpy as np
import xarray as xr

logger = logging.getLogger(__name__)


def slope_degrees(
    dem: xr.DataArray, *, radius_m: float, cell_size_m: float
) -> xr.DataArray:
    """Slope from a plane fitted over a circular neighbourhood of `radius_m`.

    The DEM must be in a projected CRS with the same linear unit vertically and
    horizontally, or the returned angle is meaningless.
    """
    if cell_size_m <= 0:
        raise ValueError("cell_size_m must be positive")
    k = int(round(radius_m / cell_size_m))
    if k < 1:
        raise ValueError(
            f"radius {radius_m} m is smaller than one {cell_size_m} m cell; "
            "either the analysis scale or the grid resolution is wrong"
        )

    # Offsets of every cell in the neighbourhood, in metres from the centre.
    ii, jj = np.mgrid[-k:k + 1, -k:k + 1]
    inside = (ii ** 2 + jj ** 2) <= k ** 2
    x = (jj * cell_size_m)[inside]
    y = (-ii * cell_size_m)[inside]

    # Least-squares plane z = ax + by + c reduces to fixed weights per offset,
    # so the fit is a pair of convolutions rather than a per-cell solve.
    denom_x = float((x ** 2).sum())
    denom_y = float((y ** 2).sum())
    wx = np.zeros(inside.shape); wx[inside] = x / denom_x
    wy = np.zeros(inside.shape); wy[inside] = y / denom_y

    dzdx = _convolve_nodata_aware(dem, wx)
    dzdy = _convolve_nodata_aware(dem, wy)
    slope = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))

    slope.attrs.update(
        long_name="seafloor slope",
        units="degree",
        analysis_radius_m=float(radius_m),
        cell_size_m=float(cell_size_m),
        method="least-squares plane fit over a circular neighbourhood",
    )
    logger.info("slope computed at radius %.1f m (%d-cell kernel)", radius_m, k)
    return slope

The attributes written at the end are not decoration. analysis_radius_m is the single field that makes the output interpretable, and it is the field that is lost every time a slope raster is produced by a menu click.

Rugosity Is Not a Slope Proxy #

The oldest rugosity measure — surface area divided by planar area — is monotonic in slope. On a uniformly steep plane it returns a large value, and on a flat but chaotic surface it returns a small one, which is the opposite of what a habitat model wants: a steep smooth wall is poor habitat, a flat boulder field is excellent. Any result computed with the surface-area ratio and interpreted as complexity is partly measuring slope.

The vector ruggedness measure (VRM) separates the two. It converts each cell’s slope and aspect into a unit normal vector, sums those vectors over the neighbourhood, and reports one minus the magnitude of the resultant divided by the count. A uniform plane, however steep, has parallel normals, a resultant of full magnitude, and a VRM of zero. A chaotic surface has dispersed normals, a short resultant, and a VRM approaching one. That is the quantity habitat models mean.

Bathymetric Position Index: an Annulus, Not a Disc #

BPI asks a different question from slope: not “how steep is it here” but “is this cell high or low relative to its surroundings”. It is computed as the depth at a cell minus the mean depth over an annulus around it, standardised by the annulus’s standard deviation, so the output is in units of local variability rather than metres.

The annulus geometry is the whole of the method, and it is where implementations diverge. The inner radius must exclude the cell’s immediate neighbourhood, or the comparison is against terrain the cell is part of and every value collapses toward zero. The outer radius sets what “surroundings” means, and therefore what features the index can see: an annulus running from 3 to 25 cells finds ridges and channels a few tens of metres across; one running from 25 to 250 cells finds banks and basins kilometres across. A feature much larger than the outer radius is invisible because the annulus sits entirely on it; a feature much smaller than the inner radius is invisible because the annulus never touches it.

That gives BPI a genuine band-pass character, and it is the reason the standard workflow computes two — a fine and a broad index — and classifies from the pair. A cell that is high on the fine index and neutral on the broad one is a ridge on a plain; high on both is a peak on a bank; low on the fine and high on the broad is a channel incised into a bank. The four-way classification from two indices is more informative than either alone and costs one extra convolution.

def bpi(dem: xr.DataArray, inner_cells: int, outer_cells: int) -> xr.DataArray:
    """Standardised bathymetric position index over an annulus."""
    if not 0 < inner_cells < outer_cells:
        raise ValueError(f"need 0 < inner < outer, got {inner_cells}, {outer_cells}")
    ii, jj = np.mgrid[-outer_cells:outer_cells + 1, -outer_cells:outer_cells + 1]
    r2 = ii ** 2 + jj ** 2
    ring = ((r2 > inner_cells ** 2) & (r2 <= outer_cells ** 2)).astype("float64")
    ring /= ring.sum()

    mean_ring = _convolve_nodata_aware(dem, ring)
    anomaly = dem - mean_ring
    sigma = float(anomaly.std(skipna=True))
    if sigma <= 0:
        raise ValueError("zero variance: the annulus is degenerate or the DEM is flat")

    out = anomaly / sigma
    out.attrs.update(
        long_name="bathymetric position index",
        units="1", inner_radius_cells=inner_cells, outer_radius_cells=outer_cells,
    )
    return out

Standardising by a single global sigma, as above, makes values comparable across the whole grid; standardising per-window instead makes them comparable only locally and turns a flat basin into a field of extreme values because its local variability is near zero. Global is almost always what is wanted, and the choice belongs in the metadata either way.

Nodata Is the Whole Difficulty #

Terrain derivatives are neighbourhood operations, and a bathymetric grid is full of holes: the coastline, unsurveyed gaps, cells masked by the gap-radius rule. Every one of those holes has a rim, and what the kernel does at the rim decides whether the product has a fabricated fringe.

Three behaviours are possible and only one is defensible by default. Filling nodata with a constant — zero, or the grid mean — invents a cliff at every hole edge and produces spectacular false slopes. Renormalising over the valid cells, as a normalised convolution does, produces a value from fewer samples; it is reasonable for smoothing and wrong for slope, because a gradient fitted to a half-neighbourhood is biased toward the side that had data. Propagating nodata — any kernel touching a hole emits nodata — shrinks the product by one kernel radius around every hole and is the only option that never states something the data cannot support.

Propagation is the right default, with one caveat worth building in: at a large analysis radius the erosion is large too, and a 250-cell BPI annulus will erase a strip 250 cells wide along the entire coast. That is honest but often useless, so the practical compromise is a minimum valid fraction: emit a value where at least, say, 70% of the neighbourhood is valid, renormalise over what is there, and record the valid fraction as a companion band so the consumer can filter further.

def _convolve_nodata_aware(
    dem: xr.DataArray, weights: np.ndarray, *, min_valid_fraction: float = 0.7
) -> xr.DataArray:
    """Weighted neighbourhood sum that neither invents nor over-erodes data."""
    from scipy.ndimage import convolve

    valid = np.isfinite(dem.values)
    filled = np.where(valid, dem.values, 0.0)

    numer = convolve(filled, weights, mode="nearest")
    # Fraction of the kernel's total weight that landed on real cells.
    frac = convolve(valid.astype("float64"), np.abs(weights), mode="nearest")
    frac /= np.abs(weights).sum()

    out = np.where(frac >= min_valid_fraction, numer, np.nan)
    return xr.DataArray(out, coords=dem.coords, dims=dem.dims)

Chunked Computation and the Halo #

At survey scale these grids do not fit in memory, and every metric here is a windowed operation, so the chunking rule from Dask vs Ray applies directly: use map_overlap with a depth of at least the kernel radius. Getting the depth wrong by one cell produces a seam grid — a faint rectangular lattice across the product, at exactly the chunk size, which is unmistakable once you know to look for it and easy to mistake for a real bedform field if you do not.

def slope_chunked(dem: xr.DataArray, radius_m: float, cell_size_m: float) -> xr.DataArray:
    k = int(round(radius_m / cell_size_m))
    return dem.map_overlap(
        lambda block: slope_degrees(block, radius_m=radius_m, cell_size_m=cell_size_m),
        depth=k, boundary=np.nan,
    )

boundary=np.nan rather than a reflection is deliberate: reflecting the surface at a chunk boundary manufactures a mirror-symmetric seabed that does not exist, and the fabricated symmetry shows up in curvature and BPI as a line of zeros.

Validation Gates #

Gate 1 — analytic recovery. Run the implementation over a synthetic plane of known dip and a cone of known half-angle. Slope must recover the plane’s dip to within the discretisation error, aspect must recover its azimuth, and curvature on the plane must be zero everywhere. This catches unit errors, sign errors and transposed axes in one test, and it runs in milliseconds.

Gate 2 — nodata containment. The count of non-nodata cells in the output must be less than or equal to the count in the input. A derivative with more valid cells than its source has invented terrain, which is the signature of a constant fill.

Gate 3 — seam detection. Compute the mean absolute difference across every chunk boundary in the output and compare it against the mean absolute difference across an equal number of interior columns. A ratio above about 1.2 means the halo is too small. This is the only gate that reliably catches an off-by-one in depth, and it costs one pass over the array.

Pipeline Integration and Downstream Handoff #

Each metric is written as its own band or its own COG, named with its analysis radius — slope_r40m.tif, not slope.tif — and carrying analysis_radius_m, cell_size_m, method and min_valid_fraction in its metadata. The valid-fraction band travels with them.

Downstream, habitat models consume rugosity and BPI, engineering studies consume slope and curvature, and both consume the uncertainty band produced upstream. Where a consumer asks for “the slope layer”, the correct response is to ask what scale their process operates at; where that question has no answer, publishing two rasters at bracketing scales is more useful than picking one and hoping.

Related #

• Computing Seafloor Slope From a DEM in Python — the Horn kernel, units, and the nodata rim in detail
• Calculating Benthic Rugosity and BPI — VRM and the annulus geometry that defines a position index
• Handling Data Gaps in Sparse Bathymetry — where the holes these kernels must propagate around come from
• Surface Smoothing Algorithms in Python — the stage immediately upstream, and how it changes every metric here
• Dask vs Ray for Bathymetric Processing at Scale — why windowed metrics belong in the chunked-array model

Up: Bathymetric Processing & Terrain Modeling