Coastal & Marine Spatial Analysis Pipelines

Calculating Benthic Rugosity and BPI

Habitat models want two things from a bathymetric surface that slope alone cannot provide: how complex the seabed is, and whether a given place is high or low relative to its surroundings. Rugosity and the bathymetric position index answer those, and both have a widely-used implementation that answers a subtly different question from the one the model needs. This page gives working implementations of both, within deriving slope, rugosity and terrain metrics.

Root cause: the legacy rugosity measure is a slope measure

The oldest rugosity metric is surface area divided by planar area. On a plane tilted at angle θ that ratio is exactly 1/cos θ, so it is a monotonic function of slope and nothing else. A steep smooth wall scores high; a flat boulder field scores low. For habitat work that ranking is backwards — the boulder field provides shelter and the wall does not.

The vector ruggedness measure (VRM) fixes it by measuring the dispersion of surface orientations rather than their magnitude. Each cell’s slope and aspect define a unit normal; over a neighbourhood those normals are summed as vectors and the length of the resultant is compared against the count. Parallel normals — any plane, at any tilt — give a resultant of full length and a VRM of zero. Scattered normals give a short resultant and a VRM approaching one.

import numpy as np

def unit_normals(slope_deg: np.ndarray, aspect_deg: np.ndarray):
    """Per-cell surface normal components from slope and aspect."""
    s, a = np.radians(slope_deg), np.radians(aspect_deg)
    return (np.sin(s) * np.sin(a),      # x
            np.sin(s) * np.cos(a),      # y
            np.cos(s))                  # z

The z component is cos(slope), so a flat cell has a normal pointing straight up and a vertical cell has one lying in the plane. That is the whole of the geometry, and everything below is a neighbourhood sum over it.

VRM Measures the Spread of Normals, Not Their Tilt Three panels each showing a seabed profile with its per-cell unit normals drawn as short arrows, and beside each the vector sum of those normals. The flat plain has all normals pointing straight up; their sum is nearly the full count in length, so VRM is near zero. The steep smooth slope has all normals tilted by the same thirty degrees; they are still parallel, the sum is still nearly full length, and VRM is again near zero despite the steepness. The boulder field is nearly level overall but its normals point in many directions; the vector sum is short, and VRM is high. A caption states that the surface-area ratio would rank the middle panel highest and the right panel lowest, which is the opposite of what a habitat model wants. flat plain resultant: full VRM ≈ 0.00 steep, smooth resultant: still full VRM ≈ 0.00 boulder field resultant: short VRM ≈ 0.34 The surface-area ratio ranks the middle panel highest and the right one lowest — the reverse of what habitat wants. VRM ranks them 0, 0, 0.34, which is the ordering the shelter argument implies.

Step-by-step implementation

1. Vector ruggedness, with nodata propagated 

import logging

import numpy as np
from scipy.ndimage import uniform_filter

logger = logging.getLogger(__name__)


def vector_ruggedness(
    slope_deg: np.ndarray, aspect_deg: np.ndarray, *, size: int = 3,
    min_valid_fraction: float = 0.7,
) -> np.ndarray:
    """Vector ruggedness measure over a size x size neighbourhood, 0 to 1."""
    if size < 3 or size % 2 == 0:
        raise ValueError("size must be an odd number >= 3")
    valid = np.isfinite(slope_deg) & np.isfinite(aspect_deg)
    nx, ny, nz = unit_normals(np.where(valid, slope_deg, 0.0),
                              np.where(valid, aspect_deg, 0.0))
    for c in (nx, ny, nz):
        c[~valid] = 0.0

    n = float(size * size)
    sx, sy, sz = (uniform_filter(c, size=size, mode="nearest") * n for c in (nx, ny, nz))
    count = uniform_filter(valid.astype("float64"), size=size, mode="nearest") * n

    with np.errstate(invalid="ignore", divide="ignore"):
        vrm = 1.0 - np.sqrt(sx ** 2 + sy ** 2 + sz ** 2) / count
    vrm[count / n < min_valid_fraction] = np.nan
    vrm = np.clip(vrm, 0.0, 1.0)
    logger.info("VRM %d×%d: median %.4f over %d valid cells",
                size, size, float(np.nanmedian(vrm)), int(np.isfinite(vrm).sum()))
    return vrm

Dividing by the valid count rather than by the window size is what makes this correct at a coastline or a data hole. Divide by the window size instead and every cell near an edge gets an artificially low resultant, which reads as high ruggedness — a rim of false complexity around every gap, which a habitat model will happily interpret as a reef.

2. BPI: two annuli, then a classification 

def annulus_kernel(inner: int, outer: int) -> np.ndarray:
    """Normalised ring kernel for a bathymetric position index."""
    if not 0 < inner < outer:
        raise ValueError(f"need 0 < inner < outer, got {inner}, {outer}")
    ii, jj = np.mgrid[-outer:outer + 1, -outer:outer + 1]
    r2 = ii ** 2 + jj ** 2
    ring = ((r2 > inner ** 2) & (r2 <= outer ** 2)).astype("float64")
    if ring.sum() == 0:
        raise ValueError("annulus contains no cells — inner and outer are too close")
    return ring / ring.sum()


def bpi_standardised(dem: np.ndarray, inner: int, outer: int) -> np.ndarray:
    """Depth minus the annulus mean, standardised by the global anomaly σ."""
    from scipy.ndimage import convolve
    valid = np.isfinite(dem)
    filled = np.where(valid, dem, 0.0)
    k = annulus_kernel(inner, outer)
    ring_mean = convolve(filled, k, mode="nearest")
    frac = convolve(valid.astype("float64"), k, mode="nearest")
    with np.errstate(invalid="ignore", divide="ignore"):
        anomaly = dem - ring_mean / np.maximum(frac, 1e-9)
    anomaly[frac < 0.7] = np.nan
    sigma = float(np.nanstd(anomaly))
    if sigma <= 0:
        raise ValueError("zero anomaly variance — the DEM is flat or the annulus is degenerate")
    return anomaly / sigma

Standardising by a single global sigma keeps values comparable across the whole grid. A per-window standardisation makes a flat basin — where local variability is near zero — into a field of extreme values, which is the most common way a BPI product ends up unusable.

3. Classify from the pair 

def benthic_class(fine: np.ndarray, broad: np.ndarray, t: float = 1.0) -> np.ndarray:
    """Four-way benthic position from a fine and a broad BPI."""
    out = np.zeros(fine.shape, dtype="uint8")          # 0 = flat / slope
    out[(fine >= t) & (np.abs(broad) < t)] = 1         # crest on a plain
    out[(fine >= t) & (broad >= t)] = 2                # crest on a high
    out[(fine <= -t) & (broad >= t)] = 3               # channel incised into a high
    out[(fine <= -t) & (np.abs(broad) < t)] = 4        # depression on a plain
    out[~np.isfinite(fine) | ~np.isfinite(broad)] = 255
    return out

The classes are the product. Raw standardised BPI values are hard for a non-specialist to act on; “channel incised into a bank” is a thing an ecologist can map against a species distribution.

Two Indices, One Map of Benthic Position A scatter space with the fine BPI on the horizontal axis and the broad BPI on the vertical, both in standard deviations. Threshold lines at plus and minus one standard deviation divide it into regions. The upper right region is labelled crest on a high, the right centre region crest on a plain, the upper left region channel incised into a high, the lower centre region depression on a plain, and the centre region flat or uniform slope. A profile beneath runs across a bank with a channel cut into it, and callouts map each part of the profile to its class: the bank flanks fall in the neutral centre, the bank crest in crest on a high, the channel floor in channel incised into a high, and the surrounding plain in flat. fine BPI (σ) → broad BPI (σ) → crest on a high crest on a plain channel in a high flat / slope the same classes on a profile flat slope crest on a high channel flat Raw standardised values are hard to act on. “Channel incised into a bank” is something an ecologist can map against a species distribution, and it costs one extra convolution to produce. Publish the classes and the two indices; the thresholds are a choice and someone will want to move them.

Verification and acceptance test 

import pytest


def test_vrm_is_zero_on_any_plane():
    """A plane at any tilt has parallel normals and therefore zero ruggedness."""
    for dip in (0.0, 5.0, 30.0, 60.0):
        slope = np.full((64, 64), dip)
        aspect = np.full((64, 64), 135.0)
        vrm = vector_ruggedness(slope, aspect, size=5)
        assert np.nanmax(vrm[3:-3, 3:-3]) < 1e-9, f"VRM non-zero on a {dip}° plane"


def test_bpi_is_zero_on_a_plane():
    yy, xx = np.mgrid[0:128, 0:128]
    dem = -20.0 - 0.01 * xx                      # uniform slope, no position structure
    b = bpi_standardised(dem, inner=3, outer=25)
    assert np.nanmax(np.abs(b[30:-30, 30:-30])) < 0.05

Both tests assert the negative — no complexity where there is none, no position anomaly on a uniform slope — and both fail loudly if a sign, an axis or a normalisation is wrong.

What these layers are used for, and what that demands of them

Rugosity and BPI are rarely the deliverable. They are predictor variables in a habitat model, and that use imposes three requirements that a standalone raster does not.

They must be comparable between surveys. A species-distribution model fitted on one survey and applied to another needs both to have been computed at the same analysis scale on grids of the same resolution. A VRM computed at 3×3 on a 1 m grid and one computed at 3×3 on a 5 m grid are different variables with the same name, and a model that mixes them is fitting the survey rather than the seabed.

Their nodata must be honest. A habitat model given a fabricated rim of high rugosity around every data hole will learn that species prefer the edges of survey coverage. Because the rim follows the coastline, and many species genuinely do prefer the coast, the artefact is confounded with the signal and will not show up as an obvious error.

Their correlation structure matters. VRM and BPI are both derived from the same DEM and are not independent. Fine-scale BPI magnitude correlates strongly with VRM, because both respond to short-wavelength structure. Feeding both into a model that assumes independent predictors inflates the apparent importance of terrain, and the standard remedy — check the variance inflation factor before fitting — belongs in the modelling step, not here, but the layers should be shipped with a note that it is needed.

The practical consequence for this stage is small and specific: publish the analysis scale, the grid resolution and the nodata rule in the layer metadata, and publish the pair of BPI annuli rather than a single derived class alone. Everything a downstream modeller needs to check comparability is then in the file.

The Fabricated Rim, and Why a Model Cannot See Through It A plan view of a survey with a coastline along the top and an unsurveyed hole in the middle. In the incorrect version, a bright rim of high rugosity follows both the coastline and the hole edge, produced by dividing the summed normals by the full window size instead of by the count of valid cells. In the correct version the rim is absent and only real seabed structure remains. Beneath, a note explains that because the rim follows the coast, and many species genuinely prefer coastal habitat, a distribution model fitted on the incorrect layer attributes real ecological signal to an artefact and cannot distinguish the two from the data alone. divide by window size hole bright rim on the coast and the hole edge divide by the valid count hole only real seabed structure remains The rim follows the coast, and many species genuinely prefer the coast. A model fitted on the left-hand layer attributes real ecology to an artefact, and cannot separate them from the data alone. One denominator. That is the whole difference.

Edge cases and gotchas

  • Aspect is undefined where slope is zero, and libraries return anything from 0 to −9999 there. Since the normal is straight up regardless of aspect when slope is zero, the formulation above is immune — but a VRM implementation that trigonometrically combines aspect before checking slope is not.
  • VRM is scale-dependent in the same way slope is. A 3×3 VRM measures cell-scale roughness and an 11×11 measures patch-scale complexity. Publish the window size in the layer name, not only in the metadata.
  • BPI annuli in cells, products in metres. The same annulus on a 1 m and a 5 m grid finds features five times apart in size. Specify the radii in metres and convert at run time, as the parent topic does for slope.

Up: Deriving Slope, Rugosity & Terrain Metrics