Stochastic depression analysis¶

Monte-Carlo uncertainty quantification on depression detection. The method adds Gaussian noise (σ = supplied vertical-error estimate) to the DEM, runs fill_depressions on each noisy realisation, and aggregates per-cell occurrence probability across N realisations.

Reference: Lindsay & Creed (2006). "Distinguishing actual and artefact depressions in digital elevation data." Computers & Geosciences 32(8): 1192–1204.

Use: identify which depressions are robust to DEM vertical error vs noise artefacts.

In [1]:

import matplotlib.pyplot as plt
plt.rcParams["figure.dpi"] = 80
plt.rcParams["savefig.dpi"] = 80

import matplotlib.pyplot as plt plt.rcParams["figure.dpi"] = 80 plt.rcParams["savefig.dpi"] = 80

In [2]:

import numpy as np
from pyramids.dataset import Dataset
from digitalrivers import DEM

# Build a DEM with one clear deep pit + one ambiguous shallow pit.
rows, cols = 20, 20
z = np.full((rows, cols), 100.0, dtype=np.float32)
z[5, 5] = 80.0    # deep pit — 20 m below surroundings
z[15, 15] = 99.0  # shallow pit — only 1 m below surroundings

ds = Dataset.create_from_array(
    z, top_left_corner=(0.0, 0.0), cell_size=30.0,
    epsg=32618, no_data_value=-9999.0,
)
dem = DEM(ds.raster)
print(f"DEM with deep pit at (5, 5)={z[5, 5]}, shallow pit at (15, 15)={z[15, 15]}")

import numpy as np from pyramids.dataset import Dataset from digitalrivers import DEM # Build a DEM with one clear deep pit + one ambiguous shallow pit. rows, cols = 20, 20 z = np.full((rows, cols), 100.0, dtype=np.float32) z[5, 5] = 80.0 # deep pit — 20 m below surroundings z[15, 15] = 99.0 # shallow pit — only 1 m below surroundings ds = Dataset.create_from_array( z, top_left_corner=(0.0, 0.0), cell_size=30.0, epsg=32618, no_data_value=-9999.0, ) dem = DEM(ds.raster) print(f"DEM with deep pit at (5, 5)={z[5, 5]}, shallow pit at (15, 15)={z[15, 15]}")

2026-05-17 22:36:34 | INFO | pyramids.base.config | Logging is configured.

DEM with deep pit at (5, 5)=80.0, shallow pit at (15, 15)=99.0

Visualise: synthetic DEM with the two pits marked¶

In [3]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(z, cmap="terrain", origin="upper")
ax.plot(5, 5, "ko", markersize=14, markerfacecolor="red", label="deep pit (z=80)")
ax.plot(15, 15, "ko", markersize=10, markerfacecolor="orange",
        label="shallow pit (z=99)")
ax.set_title("Synthetic DEM"); ax.legend()
fig.colorbar(im, ax=ax, shrink=0.7, label="elevation (m)")
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(6, 5)) im = ax.imshow(z, cmap="terrain", origin="upper") ax.plot(5, 5, "ko", markersize=14, markerfacecolor="red", label="deep pit (z=80)") ax.plot(15, 15, "ko", markersize=10, markerfacecolor="orange", label="shallow pit (z=99)") ax.set_title("Synthetic DEM"); ax.legend() fig.colorbar(im, ax=ax, shrink=0.7, label="elevation (m)") fig.tight_layout() plt.show()

Low-noise scenario¶

σ = 0.5 m vertical error (typical for high-quality LiDAR-derived DEMs). The deep pit is far below any noise realisation; the shallow pit is right at the noise scale.

In [4]:

prob_low = dem.stochastic_depressions(sigma=0.5, n_runs=50, seed=42).read_array()
print(f"Probability at deep pit    (5, 5):  {prob_low[5, 5]:.3f}")
print(f"Probability at shallow pit (15, 15): {prob_low[15, 15]:.3f}")
print(f"Probability at flat cell   (0, 0):   {prob_low[0, 0]:.3f}")

prob_low = dem.stochastic_depressions(sigma=0.5, n_runs=50, seed=42).read_array() print(f"Probability at deep pit (5, 5): {prob_low[5, 5]:.3f}") print(f"Probability at shallow pit (15, 15): {prob_low[15, 15]:.3f}") print(f"Probability at flat cell (0, 0): {prob_low[0, 0]:.3f}")

Probability at deep pit    (5, 5):  1.000
Probability at shallow pit (15, 15): 0.880
Probability at flat cell   (0, 0):   0.000

High-noise scenario¶

σ = 5 m vertical error. Now the deep pit's signal-to-noise ratio degrades — it should still register most of the time, but no longer ~100%.

In [5]:

prob_high = dem.stochastic_depressions(sigma=5.0, n_runs=50, seed=42).read_array()
print(f"Probability at deep pit    (5, 5):  {prob_high[5, 5]:.3f}")
print(f"Probability at shallow pit (15, 15): {prob_high[15, 15]:.3f}")

prob_high = dem.stochastic_depressions(sigma=5.0, n_runs=50, seed=42).read_array() print(f"Probability at deep pit (5, 5): {prob_high[5, 5]:.3f}") print(f"Probability at shallow pit (15, 15): {prob_high[15, 15]:.3f}")

Probability at deep pit    (5, 5):  1.000
Probability at shallow pit (15, 15): 0.380

Visualise: depression-occurrence probability at two noise levels¶

In [6]:

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for ax, prob, sigma in zip(axes, (prob_low, prob_high), (0.5, 5.0)):
    show = np.where(prob <= 0, np.nan, prob)
    ax.imshow(z, cmap="gray", origin="upper", alpha=0.5)
    im = ax.imshow(show, cmap="hot_r", vmin=0, vmax=1, origin="upper")
    ax.set_title(f"σ = {sigma} m"); ax.axis("off")
    fig.colorbar(im, ax=ax, shrink=0.7, label="probability")
fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt fig, axes = plt.subplots(1, 2, figsize=(10, 4)) for ax, prob, sigma in zip(axes, (prob_low, prob_high), (0.5, 5.0)): show = np.where(prob <= 0, np.nan, prob) ax.imshow(z, cmap="gray", origin="upper", alpha=0.5) im = ax.imshow(show, cmap="hot_r", vmin=0, vmax=1, origin="upper") ax.set_title(f"σ = {sigma} m"); ax.axis("off") fig.colorbar(im, ax=ax, shrink=0.7, label="probability") fig.tight_layout() plt.show()

Reproducibility¶

Same seed → bit-identical output. Useful for reproducible reports.

In [7]:

p1 = dem.stochastic_depressions(sigma=1.0, n_runs=10, seed=99).read_array()
p2 = dem.stochastic_depressions(sigma=1.0, n_runs=10, seed=99).read_array()
np.testing.assert_array_equal(p1, p2)
print("Same seed produces bit-identical output across calls.")

p1 = dem.stochastic_depressions(sigma=1.0, n_runs=10, seed=99).read_array() p2 = dem.stochastic_depressions(sigma=1.0, n_runs=10, seed=99).read_array() np.testing.assert_array_equal(p1, p2) print("Same seed produces bit-identical output across calls.")

Same seed produces bit-identical output across calls.

Summary¶

• High-confidence depressions (probability ≈ 1) are robust to DEM noise.
• Marginal cells (probability ~ 0.3–0.7) likely include artefacts.
• Cells with probability ≈ 0 never registered — clearly not pits at this noise scale.

Use the probability raster to threshold a confidence-aware depression mask before downstream fill / breach operations.