Computing Terrain Roughness and TPI

TL;DR: gdaldem TRI, gdaldem TPI and gdaldem roughness compute 3×3 ruggedness, topographic position and roughness in one command each. On a 1 m LiDAR DTM a 3×3 window measures micro-relief only, so for landform work compute TPI in Python as elevation minus the mean of an annulus of 20–200 m radius, standardize it, and classify ridges (high TPI), valleys (low TPI) and slopes or flats (near zero, split by slope).

# Context and Motivation

This guide is part of Hillshade, Slope and Aspect. Slope and aspect describe a surface point by point. Many questions need a neighbourhood view instead: is this cell higher or lower than its surroundings, how rugged is the terrain here, is this a ridge, a valley bottom or a mid-slope? Terrain indices answer that. They feed habitat models, landslide susceptibility, soil mapping, hydrological modelling and archaeological prospection. LiDAR DTMs make them very detailed — sometimes too detailed: indices designed for 30 m DEMs capture only centimetre-scale texture when applied with a 3×3 window at 1 m. The key decision is the neighbourhood scale.

# Prerequisites and Assumptions

• A void-filled bare-earth DTM in a projected CRS with metric units.
• GDAL command-line tools; NumPy, SciPy and rasterio for multi-scale indices.
• A target landform scale: field-scale features (20–50 m), hillslopes (100–300 m) or regional landforms (500 m+).

# Step-by-Step Implementation

# Step 1 — Quick 3×3 indices with gdaldem

TRI (mean absolute difference to the 8 neighbours), TPI (difference to the mean of neighbours) and roughness (max − min in the window). Useful for micro-relief and data-quality checks.

# Step 2 — Choose landform radii

Pick inner and outer radii for an annulus that matches the features of interest; excluding the centre avoids the cell’s own neighbourhood dominating.

# Step 3 — Compute annulus TPI

Convolve the DTM with a normalized annulus kernel (FFT convolution is fast for large radii) and subtract.

# Step 4 — Standardize

Divide by the TPI’s standard deviation so values are comparable across scales and areas.

# Step 5 — Classify

Standardized TPI above +1: ridge or upper slope; below −1: valley or lower slope; between: mid-slope if slope > 5°, flat otherwise.

# Complete Working Example

Command line, 3×3 indices:

gdaldem TRI mosaic/dtm_1m.tif out/tri_3x3.tif -compute_edges
gdaldem TPI mosaic/dtm_1m.tif out/tpi_3x3.tif -compute_edges
gdaldem roughness mosaic/dtm_1m.tif out/rough_3x3.tif -compute_edges

Python, multi-scale annulus TPI and a simple landform classification:

"""Multi-scale annulus TPI and a slope-position classification from a DTM."""
from __future__ import annotations

import numpy as np
import rasterio
from scipy.signal import fftconvolve


def annulus(inner_cells: int, outer_cells: int) -> np.ndarray:
    y, x = np.ogrid[-outer_cells:outer_cells + 1, -outer_cells:outer_cells + 1]
    r = np.hypot(x, y)
    k = ((r >= inner_cells) & (r <= outer_cells)).astype(float)
    return k / k.sum()


def tpi(z: np.ndarray, inner_m: float, outer_m: float, res: float) -> np.ndarray:
    k = annulus(max(1, int(inner_m / res)), int(outer_m / res))
    valid = ~np.isnan(z)
    zf = np.where(valid, z, 0.0)
    num = fftconvolve(zf, k, mode="same")
    den = fftconvolve(valid.astype(float), k, mode="same")
    mean = num / np.maximum(den, 1e-6)
    return np.where(valid, z - mean, np.nan)


with rasterio.open("mosaic/dtm_1m.tif") as ds:
    z = ds.read(1, masked=True).filled(np.nan).astype(float)
    profile, res = ds.profile, ds.res[0]

tpi_small = tpi(z, 5, 30, res)        # field-scale features
tpi_large = tpi(z, 50, 250, res)      # hillslope position
std = (tpi_large - np.nanmean(tpi_large)) / np.nanstd(tpi_large)

gy, gx = np.gradient(z, res)
slope = np.degrees(np.arctan(np.hypot(gx, gy)))
cls = np.full(z.shape, 0, dtype="uint8")
cls[std > 1] = 1                                   # ridge / upper slope
cls[std < -1] = 4                                  # valley / lower slope
mid = (std >= -1) & (std <= 1)
cls[mid & (slope > 5)] = 2                         # mid-slope
cls[mid & (slope <= 5)] = 3                        # flat
cls[np.isnan(z)] = 255

profile.update(dtype="uint8", nodata=255)
with rasterio.open("out/slope_position.tif", "w", **profile) as out:
    out.write(cls, 1)
print({name: f"{np.mean(cls == v):.1%}" for v, name in
       {1: "ridge", 2: "mid-slope", 3: "flat", 4: "valley"}.items()})

FFT convolution handles a 250-cell radius on a 5000 × 5000 raster in seconds; a direct convolution would take hours.

# Key Parameter Table

# Verification

• Known landforms. Ridge crests and valley bottoms visible in a hillshade should fall into the right classes at the large scale.
• Scale sanity. The small-scale TPI should highlight banks and ditches; the large-scale TPI should not respond to them.
• Edge effects. Within one outer radius of the raster edge, the annulus is truncated; the normalized convolution handles it, but classes there are less reliable. Compute on a mosaic larger than the area of interest.

# Gotchas and Edge Cases

3×3 on LiDAR measures noise. At 1 m, 3×3 TRI and roughness respond to ploughing, vegetation residue and interpolation texture. They are good data-quality indicators and poor landform indicators.

Buildings in the DTM. Unfiltered buildings create strong positive TPI spikes. Use a clean bare-earth DTM.

Non-comparable raw TPI. Raw TPI grows with relief and radius, so thresholds on raw values do not transfer. Standardize before classifying.

Memory. FFT convolution needs several copies of the raster in float64. For county-scale mosaics, process in overlapping blocks with a margin of one outer radius.

# Frequently Asked Questions

What is the topographic position index?

The difference between a cell’s elevation and the mean elevation of its neighbourhood. Positive values mean the cell is higher than its surroundings, as on ridges; negative values mean lower, as in valleys; values near zero indicate uniform slopes or flat ground.

Which neighbourhood size should I use for TPI?

One matched to the landforms of interest. Tens of metres pick out field-scale features such as banks and ditches; hundreds of metres identify hillslope positions. Using several scales and comparing them is common.

What is the difference between TRI and roughness?

TRI is the mean absolute elevation difference between a cell and its eight neighbours; roughness is the range, maximum minus minimum, within the 3 by 3 window. Both measure local ruggedness.

Why standardize TPI?

Raw TPI values depend on relief and neighbourhood size, so thresholds do not transfer between areas or scales. Dividing by the standard deviation makes values comparable and lets fixed thresholds classify landforms.

# Related

• Hillshade, Slope and Aspect — terrain derivatives overview
• Generating Slope and Aspect Rasters with gdaldem — the slope input to classification
• Multidirectional Hillshade for LiDAR DTMs — visual interpretation alongside indices
• Smoothing a LiDAR DTM Before Contouring — normalized convolution in another context
• Building a Seamless DTM Mosaic from Tiles — avoiding edge effects