Coastal & Marine Spatial Analysis Pipelines

Computing Seafloor Slope From a DEM in Python

Slope is the first derived layer anyone asks for from a bathymetric grid, it is four lines of NumPy, and the four lines are wrong in a specific way on roughly half the marine DEMs they are first tried on. This page is the correct implementation, the test that proves it, and the two traps that make an incorrect one look plausible. It is a step within deriving slope, rugosity and terrain metrics, which covers why the scale of the calculation matters more than the kernel used for it.

Root cause: slope is a ratio of two units, and marine grids mix them

Slope is rise over run. The rise comes from the DEM’s values, in metres of depth; the run comes from its coordinates. On a projected grid both are metres and the ratio is dimensionless, as it must be. On a geographic grid — EPSG:4326, which is what a great many delivered bathymetric products are in — the run is in degrees, and a degree of longitude is about 111,320 m at the equator and about 78,800 m at 45° N.

The result is a slope raster that is too small by five orders of magnitude, varies with latitude, and renders as an entirely plausible near-flat seabed:

import numpy as np
import rioxarray as rxr

dem = rxr.open_rasterio("bathy_wgs84.tif", masked=True).squeeze()
print(dem.rio.crs)                       # EPSG:4326
gy, gx = np.gradient(dem.values)         # per pixel, and a pixel is 0.0001 deg
slope = np.degrees(np.arctan(np.hypot(gx, gy)))
print(float(np.nanmax(slope)))           # 0.0000042  — a "flat" continental slope

Nothing raises. The array is finite, the values are in a valid range for an angle, and the map is a uniform pale wash that reads as calm seabed. The only symptom is that the number is absurd, and absurdity is only visible to someone who knows what a continental slope should measure.

The Same Seabed, Two Grids, Two Incompatible Answers A cross-section of a seabed dipping at 30 degrees is drawn once. Beneath it, two calculations are shown. On a UTM grid with 10 metre cells, the rise of 5.77 metres over a run of 10 metres gives a gradient of 0.577 and a slope of 30.0 degrees. On a geographic grid with 0.0001 degree cells at latitude 45 north, the same rise of 5.77 metres is divided by a run recorded as 0.0001, giving a gradient of 57,700 if interpreted directly, or, as the code more commonly does, dividing metres by degrees to yield a number five orders of magnitude away from the truth. A note adds that because a degree of longitude shortens with latitude, the error is not even constant across one raster. run rise 5.77 m one seabed, dipping 30° projected grid, 10 m cells 5.77 m / 10 m = 0.577 → arctan → 30.0° both terms in metres; the ratio is dimensionless geographic grid, 0.0001° cells 5.77 m / 0.0001 deg → 0.0000042° reported metres divided by degrees — not a ratio at all and the run shortens with latitude, so it drifts The wrong answer is finite, smooth, in a valid angular range, and renders as an unremarkable map. Only its magnitude gives it away. Which is why the fix begins with an assertion rather than with arithmetic.

Step-by-step fix with production code

1. Refuse to compute on a geographic grid 

import logging

import numpy as np
import xarray as xr

logger = logging.getLogger(__name__)


def _assert_projected_metres(dem: xr.DataArray) -> tuple[float, float]:
    """Return (cell_x, cell_y) in metres, or raise with an actionable message."""
    crs = dem.rio.crs
    if crs is None:
        raise ValueError("DEM has no CRS; slope is undefined without one")
    if crs.is_geographic:
        raise ValueError(
            f"DEM is in {crs.to_string()} (degrees). Slope is a ratio of vertical "
            "to horizontal distance and both must share a unit. Reproject to a "
            "local projected CRS first — do not scale degrees to metres, because "
            "the scale factor varies with latitude across the grid."
        )
    unit = crs.linear_units.lower()
    if unit not in {"metre", "meter", "m"}:
        raise ValueError(f"expected metres, got {crs.linear_units!r}; convert first")
    t = dem.rio.transform()
    return abs(t.a), abs(t.e)

The error message matters more than the check. The tempting “fix” for a geographic grid is to multiply by 111,320, and the message says why that is also wrong: the factor varies across the raster, so the corrected slope is right in the middle and increasingly wrong toward the edges.

2. Apply the Horn kernel

Horn’s 3×3 weighted difference is the standard used by GDAL, QGIS and ArcGIS, which means a slope raster computed this way is comparable with one produced by any of them. The weights give the orthogonal neighbours twice the influence of the diagonals.

from scipy.ndimage import convolve

HORN_X = np.array([[-1.0, 0.0, 1.0], [-2.0, 0.0, 2.0], [-1.0, 0.0, 1.0]]) / 8.0
HORN_Y = np.array([[1.0, 2.0, 1.0], [0.0, 0.0, 0.0], [-1.0, -2.0, -1.0]]) / 8.0


def slope_degrees(dem: xr.DataArray) -> xr.DataArray:
    """Horn slope in degrees, with nodata propagated rather than filled."""
    cell_x, cell_y = _assert_projected_metres(dem)

    z = dem.values.astype("float64")
    valid = np.isfinite(z)
    filled = np.where(valid, z, 0.0)

    dzdx = convolve(filled, HORN_X, mode="nearest") / cell_x
    dzdy = convolve(filled, HORN_Y, mode="nearest") / cell_y

    # Any output cell whose 3x3 window touched a hole is itself unknown.
    touched_hole = convolve((~valid).astype("float64"), np.ones((3, 3)), mode="nearest") > 0
    grad = np.hypot(dzdx, dzdy)
    out = np.where(touched_hole, np.nan, np.degrees(np.arctan(grad)))

    result = xr.DataArray(out, coords=dem.coords, dims=dem.dims)
    result.attrs.update(
        long_name="seafloor slope", units="degree", method="Horn 3x3",
        cell_size_x_m=cell_x, cell_size_y_m=cell_y,
    )
    logger.info(
        "slope: %d/%d cells valid, max %.2f deg",
        int(np.isfinite(out).sum()), out.size, float(np.nanmax(out)),
    )
    return result

Note the sign convention on HORN_Y. Raster rows increase downward while northing increases upward, so the y kernel is flipped relative to the x kernel. Getting this wrong leaves slope magnitude correct and aspect reflected about the east-west axis — a mistake that survives every check except an aspect test, which is why the analytic test below checks aspect too.

3. Propagate nodata, do not fill it

The touched_hole line is the difference between a product that stops at the coastline and one that grows a fringe of fabricated cliffs along it. Filling holes with zero before convolving turns every land cell into a depth of 0 m, so the kernel at the shoreline sees a 30 m step over one cell and reports a near-vertical slope for the entire coast.

What the Kernel Does at the Coastline A profile crossing from water into land, with the land side marked as nodata. Three slope traces are drawn beneath. The zero-fill trace spikes to near 90 degrees at the boundary, because the kernel sees a step from 30 metres depth to zero across one cell. The renormalised trace produces a moderate but biased value, because the gradient is fitted only to the water side of the window. The propagate trace stops one cell short of the boundary and emits nodata beyond, showing a small gap and no fabricated value. A caption notes the zero-fill artefact is a continuous ribbon along every coastline and every survey hole, and that it reads convincingly as a scarp. land — nodata seabed, ~30 m depth zero fill → 88° cliff, along the whole coast renormalise → a value, fitted to half a window propagate → stops one cell short, says nothing false The zero-fill artefact is a continuous ribbon along every coastline and every survey hole, and it looks exactly like a scarp.

Verification and acceptance test

A synthetic surface with a known answer tests units, kernel weights, sign conventions and cell-size handling in one run, and it takes milliseconds.

import pytest
from affine import Affine


def _synthetic_plane(dip_deg: float, azimuth_deg: float, cell: float = 10.0, n: int = 64):
    """A plane of exact dip and aspect, on a projected grid."""
    g = np.tan(np.radians(dip_deg))
    a = np.radians(azimuth_deg)
    jj, ii = np.meshgrid(np.arange(n), np.arange(n))
    x, y = jj * cell, -ii * cell
    z = -(x * np.sin(a) + y * np.cos(a)) * g
    da = xr.DataArray(z, dims=("y", "x"),
                      coords={"y": -np.arange(n) * cell, "x": np.arange(n) * cell})
    da.rio.write_crs("EPSG:32618", inplace=True)
    da.rio.write_transform(Affine(cell, 0, 0, 0, -cell, 0), inplace=True)
    return da


@pytest.mark.parametrize("dip", [0.5, 5.0, 30.0, 60.0])
def test_slope_recovers_a_known_plane(dip):
    dem = _synthetic_plane(dip, azimuth_deg=135.0)
    s = slope_degrees(dem)
    interior = s.values[2:-2, 2:-2]          # avoid the propagated edge
    assert np.allclose(interior, dip, atol=1e-6), f"recovered {np.nanmean(interior)}"


def test_slope_refuses_a_geographic_grid():
    dem = _synthetic_plane(10.0)
    dem.rio.write_crs("EPSG:4326", inplace=True)
    with pytest.raises(ValueError, match="degrees"):
        slope_degrees(dem)

A plane is exactly recoverable by Horn’s kernel, so atol=1e-6 is a real assertion rather than a loose one — any unit, weight or cell-size error moves the result far outside it.

4. Decide what “slope” means at a chunk boundary

At survey scale the DEM is chunked, and Horn’s kernel needs one cell of context on every side. Running it per chunk without a halo produces a one-cell band of wrong values along every internal boundary — a faint rectangular lattice across the whole product at exactly the chunk size.

def slope_chunked(dem: xr.DataArray) -> xr.DataArray:
    """Horn slope over a Dask-backed DEM, with the halo the kernel requires."""
    if dem.chunks is None:
        return slope_degrees(dem)
    return dem.map_overlap(slope_degrees, depth=1, boundary=np.nan)

depth=1 is the kernel radius, not a safety margin, and boundary=np.nan rather than "reflect" is deliberate: reflecting the surface at the outer edge of the grid manufactures a mirror-symmetric seabed, and slope computed across the mirror line is exactly zero — a thin line of artificially flat seabed around the whole product.

The seam is easy to test for and almost impossible to see by eye once the raster is stretched for display:

def assert_no_chunk_seams(slope: xr.DataArray, chunk_cols: int, tol: float = 1.2) -> None:
    """Mean |difference| across chunk boundaries vs across interior columns."""
    v = slope.values
    edges = np.arange(chunk_cols, v.shape[1] - 1, chunk_cols)
    if edges.size == 0:
        return
    at_edge = np.nanmean(np.abs(v[:, edges] - v[:, edges + 1]))
    interior = np.arange(chunk_cols // 2, v.shape[1] - 1, chunk_cols)
    inside = np.nanmean(np.abs(v[:, interior] - v[:, interior + 1]))
    if inside > 0 and at_edge / inside > tol:
        raise ValueError(
            f"chunk seams detected: boundary discontinuity is {at_edge / inside:.2f}× "
            "the interior, so the overlap depth is too small"
        )
The Lattice That Is Not a Bedform Field Two renderings of the same slope raster. On the left, computed chunk by chunk with no overlap, a faint regular grid of lines runs across the whole image at the chunk spacing; each line is one cell wide and holds slope values computed from a truncated kernel. On the right, computed with an overlap depth of one cell, the lattice is absent and only the seabed's own structure remains. A note observes that the lattice is regular, low-amplitude and continuous, which is what a real field of sand waves also looks like, and that the distinguishing feature is that it aligns exactly with the chunk grid rather than with any current direction. depth=0 — a lattice at the chunk size regular, low-amplitude, continuous — exactly like sand waves depth=1 — only the seabed the same data, one cell of overlap Tell them apart by orientation: the lattice aligns with the chunk grid, not with any current that ever ran there.

Edge cases and gotchas

  • Anisotropic cells. A grid with 2 m easting cells and 3 m northing cells needs the two partials divided by different lengths, as the code above does. A single cell_size parameter silently skews aspect and inflates slope along one axis; the symptom is a slope map with a directional grain that follows the grid rather than the seabed.
  • Vertical exaggeration in the source. Some legacy products store depth in decimetres or feet to fit an integer type. The CRS says nothing about the vertical unit, so the check above cannot catch it. Assert that the DEM’s value range matches the survey’s expected depth range before computing anything.
  • Slope on a smoothed DEM is a different quantity. If the grid has been through the smoothing stage, the derived slope reflects the smoothing kernel as much as the seabed. Record which surface the slope came from, and prefer computing derivatives from the unsmoothed grid with an explicit analysis radius instead.

Up: Deriving Slope, Rugosity & Terrain Metrics