Out-of-core flow accumulation & slope¶

The same engine= switch streams flow accumulation (Barnes 2017 perimeter-graph, D8/Rho8) and local slope tile-by-tile. Both match the in-memory result (accumulation to floating-point tolerance; slope bit-for-bit), at constant memory.

In [1]:

import matplotlib.pyplot as plt
# keep the executed notebook small
plt.rcParams["figure.dpi"] = 80
plt.rcParams["savefig.dpi"] = 80

import matplotlib.pyplot as plt # keep the executed notebook small plt.rcParams["figure.dpi"] = 80 plt.rcParams["savefig.dpi"] = 80

Setup¶

In [2]:

import gzip, os, shutil, tempfile, urllib.request
from pathlib import Path

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


def _data_dir() -> Path:
    """Cache downloads under <project>/examples/data/ (project root marked by pyproject.toml)."""
    cwd = Path.cwd().resolve()
    for parent in (cwd, *cwd.parents):
        if (parent / "pyproject.toml").exists():
            d = parent / "examples" / "data"
            d.mkdir(parents=True, exist_ok=True)
            return d
    d = cwd / "data"
    d.mkdir(exist_ok=True)
    return d


# Download a real 1-arc-second SRTM tile (AWS open data, no auth), cached locally.
DATA = _data_dir()
hgt, gz = DATA / "N45W122.hgt", DATA / "N45W122.hgt.gz"
URL = "https://elevation-tiles-prod.s3.amazonaws.com/skadi/N45/N45W122.hgt.gz"
if not hgt.exists():
    if not gz.exists():
        print(f"Downloading {URL} ...")
        urllib.request.urlretrieve(URL, gz)
    with gzip.open(gz, "rb") as fi, open(hgt, "wb") as fo:
        shutil.copyfileobj(fi, fo)

# Crop a ~0.1 degree (360x360) window around Mount Hood and save it as a file-backed GeoTIFF DEM.
raw = np.fromfile(hgt, dtype=">i2").reshape(3601, 3601).astype(np.float32)
crop = raw[940:1300, 1080:1440].copy()
WORK = Path(tempfile.mkdtemp(prefix="dr_ooc_"))
DEM_PATH = WORK / "mt_hood_dem.tif"
Dataset.create_from_array(
    crop, top_left_corner=(-121.9, 45.5), cell_size=1 / 3600, epsg=4326,
    no_data_value=-32768.0, driver_type="GTiff", path=str(DEM_PATH),
).close()
dem = DEM(Dataset.read_file(str(DEM_PATH)).raster)
print(f"DEM: {dem.rows} x {dem.columns} cells, elevation {crop.min():.0f}-{crop.max():.0f} m")

import gzip, os, shutil, tempfile, urllib.request from pathlib import Path import numpy as np from pyramids.dataset import Dataset from digitalrivers import DEM def _data_dir() -> Path: """Cache downloads under <project>/examples/data/ (project root marked by pyproject.toml).""" cwd = Path.cwd().resolve() for parent in (cwd, *cwd.parents): if (parent / "pyproject.toml").exists(): d = parent / "examples" / "data" d.mkdir(parents=True, exist_ok=True) return d d = cwd / "data" d.mkdir(exist_ok=True) return d # Download a real 1-arc-second SRTM tile (AWS open data, no auth), cached locally. DATA = _data_dir() hgt, gz = DATA / "N45W122.hgt", DATA / "N45W122.hgt.gz" URL = "https://elevation-tiles-prod.s3.amazonaws.com/skadi/N45/N45W122.hgt.gz" if not hgt.exists(): if not gz.exists(): print(f"Downloading {URL} ...") urllib.request.urlretrieve(URL, gz) with gzip.open(gz, "rb") as fi, open(hgt, "wb") as fo: shutil.copyfileobj(fi, fo) # Crop a ~0.1 degree (360x360) window around Mount Hood and save it as a file-backed GeoTIFF DEM. raw = np.fromfile(hgt, dtype=">i2").reshape(3601, 3601).astype(np.float32) crop = raw[940:1300, 1080:1440].copy() WORK = Path(tempfile.mkdtemp(prefix="dr_ooc_")) DEM_PATH = WORK / "mt_hood_dem.tif" Dataset.create_from_array( crop, top_left_corner=(-121.9, 45.5), cell_size=1 / 3600, epsg=4326, no_data_value=-32768.0, driver_type="GTiff", path=str(DEM_PATH), ).close() dem = DEM(Dataset.read_file(str(DEM_PATH)).raster) print(f"DEM: {dem.rows} x {dem.columns} cells, elevation {crop.min():.0f}-{crop.max():.0f} m")

2026-06-09 23:51:02 | INFO | pyramids.base.config | Logging is configured.

DEM: 360 x 360 cells, elevation 23-1195 m

Tiled flow accumulation¶

Fill, derive D8 flow direction, then accumulate both in-memory and tiled.

In [3]:

filled = dem.fill_depressions(engine="in_memory")
fdir = filled.flow_direction(method="d8")

acc_mem = fdir.accumulate()
acc_tiled = fdir.accumulate(
    engine="tiled", out_path=str(WORK / "acc_tiled.tif"), tile_size=128
)
am = np.asarray(acc_mem.read_array()).astype(float)
at = np.asarray(acc_tiled.read_array()).astype(float)
print("tiled accumulation == in-memory (allclose):", np.allclose(am, at))
print("max upstream cells:", int(am.max()))

filled = dem.fill_depressions(engine="in_memory") fdir = filled.flow_direction(method="d8") acc_mem = fdir.accumulate() acc_tiled = fdir.accumulate( engine="tiled", out_path=str(WORK / "acc_tiled.tif"), tile_size=128 ) am = np.asarray(acc_mem.read_array()).astype(float) at = np.asarray(acc_tiled.read_array()).astype(float) print("tiled accumulation == in-memory (allclose):", np.allclose(am, at)) print("max upstream cells:", int(am.max()))

tiled accumulation == in-memory (allclose): True
max upstream cells: 4000

In [4]:

# Render the log-scaled accumulation through pyramids (wrap the derived array in a Dataset).
log_acc = np.log1p(at).astype(np.float32)
g = Dataset.dataset_like(acc_tiled, log_acc).plot(
    band=0, title="log(flow accumulation) — tiled", cmap="cubehelix", figsize=(5, 5)
)
g.ax.set_xticks([])
g.ax.set_yticks([])

# Render the log-scaled accumulation through pyramids (wrap the derived array in a Dataset). log_acc = np.log1p(at).astype(np.float32) g = Dataset.dataset_like(acc_tiled, log_acc).plot( band=0, title="log(flow accumulation) — tiled", cmap="cubehelix", figsize=(5, 5) ) g.ax.set_xticks([]) g.ax.set_yticks([])

Out[4]:

[]

DEM.flow_accumulation forwards engine= too¶

The convenience dispatcher streams the float32 Accumulation straight to disk for the tiled engine (no whole-array int32 cast).

In [5]:

acc2 = filled.flow_accumulation(
    fdir, engine="tiled", out_path=str(WORK / "acc2.tif"), tile_size=128
)
print("matches:", np.allclose(np.asarray(acc2.read_array()).astype(float), at))

acc2 = filled.flow_accumulation( fdir, engine="tiled", out_path=str(WORK / "acc2.tif"), tile_size=128 ) print("matches:", np.allclose(np.asarray(acc2.read_array()).astype(float), at))

matches: True

Tiled slope (a local stencil)¶

Local terrain ops are stencils — they tile trivially with a one-cell halo. DEM.slope streams tile-by-tile and is bit-for-bit identical to the whole-array slope.

In [6]:

slope_mem = dem.slope(engine="in_memory")
slope_tiled = dem.slope(
    engine="tiled", out_path=str(WORK / "slope_tiled.tif"), tile_size=128
)
sm = np.asarray(slope_mem.read_array())
st = np.asarray(slope_tiled.read_array())
print("tiled slope == in-memory slope:", np.allclose(np.nan_to_num(sm), np.nan_to_num(st)))

# slope_tiled is itself a Dataset — render it directly through pyramids.
g = slope_tiled.plot(band=0, title="max downhill slope — tiled", cmap="magma", figsize=(5, 5))
g.ax.set_xticks([])
g.ax.set_yticks([])

slope_mem = dem.slope(engine="in_memory") slope_tiled = dem.slope( engine="tiled", out_path=str(WORK / "slope_tiled.tif"), tile_size=128 ) sm = np.asarray(slope_mem.read_array()) st = np.asarray(slope_tiled.read_array()) print("tiled slope == in-memory slope:", np.allclose(np.nan_to_num(sm), np.nan_to_num(st))) # slope_tiled is itself a Dataset — render it directly through pyramids. g = slope_tiled.plot(band=0, title="max downhill slope — tiled", cmap="magma", figsize=(5, 5)) g.ax.set_xticks([]) g.ax.set_yticks([])

tiled slope == in-memory slope: True

Out[6]:

[]

In [7]:

for d in (filled, acc_mem, acc_tiled, acc2, slope_mem, slope_tiled):
    d.close()

for d in (filled, acc_mem, acc_tiled, acc2, slope_mem, slope_tiled): d.close()