feat: out-of-core (tiled, larger-than-RAM) DEM hydrology — fill, accumulation, local ops (B0–B7)

[MAfarrag]

Summary

Make digital-rivers able to process DEMs larger than RAM (e.g. a 10 GB / ~2.7 B-cell DEM, which the current
whole-array engine needs ~100–150 GB/stage to handle). The approach is the Barnes 2016/2017 tiled master-graph
method — process one tile at a time, reconcile only tile perimeters through a small global graph, write results
tile-by-tile — built on pyramids' chunked tile I/O and reusing our existing Numba kernels verbatim.

This issue collects the digital-rivers tasks (Workstream B: B0–B7). The upstream pyramids primitives these
depend on (Workstream A: A1–A4 — create_empty, map_overlap, store_windows, tiled_reduce) are tracked
separately in serapeum-org/pyramids#470. Maps to roadmap Phase 4 P30 (Dask / chunked-tile backend).

⚠️ Compiled against pyramids v0.20.0 and a DR snapshot; file:line anchors and primitive signatures have
since drifted — re-confirm with grep before implementing. Open questions are listed at the bottom.

Tasks

• B0 — wire local terrain ops (slope/aspect/hillshade/curvature/openness/TPI/TRI/color-relief) to blockwise tiling — quick win, independent — deps: A2 (soft)
• B1 — _outofcore scaffolding (TileSpec/plan_tiles/read_tile/TileStore/gid) — deps: A1
• B2 — labelled priority-flood Numba kernel (filled + labels + spillover graph) — deps: none
• B3 — tiled fill_depressions (Barnes 2016 master-graph, ε=0) — deps: A1, B1, B2
• B4 — tiled flow_accumulation (Barnes 2017 perimeter-graph, D8/ρ8) — deps: A1, B1
• B5 — public engine switch (auto/in_memory/tiled) + auto-RAM heuristic + guardrails — deps: B3, B4
• B6 — cross-seam flat / ε-fill resolution (lift the ε=0 restriction) — research — deps: B3
• B7 — dask-distributed execution under pyramids.configure(client=) — deps: B3, B4, A3 (opt)

---

0. Scope & non-goals

In scope. Out-of-core, single-machine-first (dask-distributable later) versions of the two global
DEM-hydrology operations that do not fit in RAM at scale:

• Tiled depression fill (fill_depressions) — Barnes 2016.
• Tiled flow accumulation (flow_accumulation / FlowDirection.accumulate) — Barnes 2017, D8/ρ8 only.
• Blockwise local ops (slope, aspect, hillshade, curvature, openness, TPI/TRI, color-relief) — trivial win
  via halo tiling; the cheap first deliverable.

Out of scope (flagged, deferred).

• Monotonic ε-fill (epsilon > 0) and flat-direction resolution across tile seams — these do not compose
  blockwise (see §2.3); tiled fill is correct for epsilon = 0 only, until a cross-seam flat solver is built (B6).
• Tiled D∞ / MFD accumulation — divergent flow has no fixed-halo closure; unsolved by Barnes 2017. Divergent
  routing stays in-memory-only (§2.4).
• Watershed / basins / stream-order tiling — same machinery applies later; not in this issue.

1. The decision (recap)

digital-rivers is whole-array in-RAM; a 10 GB DEM needs ~100–150 GB per stage. Local ops tile trivially;
global flow routing does not (flow crosses every tile edge). Of the five out-of-core strategies surveyed, the
tiled + global boundary-graph method (Barnes 2016/2017) is the right fit: single-machine, trillion-cell
proven, reuses our existing Numba kernels verbatim, and maps onto pyramids' tile I/O. Chosen over the
iterative-sweep alternative (§2.1).

2. Key design decisions

2.1 Master-graph, not iterative-sweep. Barnes's two-MapReduce design (map: per-tile solve → reduce: one
global graph solve → map: per-tile finalize) is single-pass, exact, and bit-identical to the single-tile
result, with fixed perimeter-sized per-tile communication. Iterative sweeps converge only with unbounded passes
and give no fixed bound.

2.2 Reuse the in-memory kernels as both the per-tile AND per-graph engine. Fill: the spillover graph "is
itself a DEM" → run the same priority-flood on the graph in the reduce step (per-tile step adds labels + spill
graph: B2). Accumulation: the global flow graph is summed with the same Kahn sweep (per-tile step adds a
perimeter FOLLOWPATH trace: B4). No new graph math.

2.3 epsilon = 0 scope; flats are the #1 correctness trap. epsilon = 0 plateau-fill composes correctly
across seams (one spill elevation per straddling depression → both halves fill to the same level → bit-identical).
Monotonic ε-fill and Garbrecht–Martz flat-direction resolution do NOT localize and silently diverge from the
single-tile result under naive tiling. Tiled fill therefore guarantees correctness for epsilon = 0 only; ε>0 /
flat routing must run on the merged surface (B6) or stay on the in-memory engine. Enforced: raise if
engine="tiled" + epsilon>0 until B6 lands.

2.4 Tiled accumulation is D8/ρ8 only. Barnes 2017 is restricted to non-divergent flow (each cell has one
downstream neighbour, so the global graph stays a forest). D∞/MFD fan flow to multiple neighbours — the perimeter
graph becomes a weighted DAG with no fixed-halo closure. Keep D∞/MFD in-memory-only, with a clear error on
engine="tiled" + divergent routing.

2.5 Memory modes: EVICT default, RETAIN fast path, CACHE later. EVICT (2R/1W): discard per-tile
intermediates, recompute them in finalize — smallest footprint, the safe default. RETAIN (1R/1W): keep per-tile
outputs in a dict (only when the dataset fits RAM). CACHE (3R/3W): spill intermediates to a scratch Zarr/GeoTIFF;
add after profiling.

2.6 Single-machine first, dask later. Stages 1 and 3 are embarrassingly parallel; stage 2 is a cheap serial
graph solve on the producer. Implement serial/threaded first; swap the consumer loop for dask tasks under
pyramids.configure(client=...) later (B7).

3. Dependency on pyramids (Workstream A — serapeum-org/pyramids#470)

A-task	What it gives B	Needed by	Status
A1 create_empty / empty_like	disk-backed empty raster sink for write_array(window=)	B1, B3, B4, B7	hard blocker
A2 public map_overlap	clean halo-aware local-op streaming	B0	recommended
A3 store_windows	chunk-parallel GeoTIFF dask sink	B7	optional/later
A4 tiled_reduce	extracted map→reduce→map scaffold	(none yet)	optional/later

A1 is the only hard blocker. Until it lands, B1/B3/B4 cannot write out-of-core. See pyramids#470 for the full
A-task specs and acceptance criteria.

Per-task files & touchpoints (what to create / edit)

Task	Create / edit	Reuse	Test
B0	dem.py (slope/aspect/hillshade/curvature/openness/TPI/TRI), terrain.py (color-relief)	A2 map_overlap; existing per-block kernels unchanged	tests/dem/, tests/terrain/ (tiled == whole-array)
B1	NEW src/digitalrivers/_outofcore/{__init__,tiling,cache}.py	mirror cloud_io.tile_windows (cloud_io.py:23); pyramids read_array(window=) / to_zarr	tests/_outofcore/test_tiling.py
B2	_numba.py (new priority_flood_labels_numba)	extend priority_flood_numba + _heap_push/_heap_pop	tests/_outofcore/test_labelled_flood.py
B3	NEW _outofcore/{spillgraph,fill}.py	B2 kernel; A1 sink; B1 tiling	tests/_outofcore/test_fill_tiled.py + tests/e2e/
B4	NEW _outofcore/{spillgraph,accumulate}.py (+ follow_path trace)	kahn_accumulate_d8_numba; A1 sink; B1 tiling	tests/_outofcore/test_accum_tiled.py + tests/e2e/
B5	dem.py, flow_direction.py (add engine= / out_path=); _resolve_engine helper	psutil (add as optional dep)	tests/dem/, tests/flow_direction/
B6	_conditioning/flats.py (seam-aware), _outofcore/fill.py	resolve_flats	tests/_outofcore/ (ε>0 seam equivalence)
B7	_outofcore/{fill,accumulate}.py (swap serial loop for dask)	pyramids.configure(client=), default_lock(); A3 sink	tests/_outofcore/ (distributed == serial)

---

4. Tasks

New internal subpackage src/digitalrivers/_outofcore/ (underscore = internal, per CLAUDE.md).

B0 — wire local ops to blockwise (QUICK WIN, independent)

Context. Local terrain derivatives are stencils — each output cell depends only on a fixed neighbourhood —
so they tile trivially with a halo and need no boundary reconciliation. Wiring them to map_overlap is the
cheapest, lowest-risk memory win, fully independent of the Barnes work. Same model as GDAL's gdaldem
(line/block streaming over a 3×3 window).

Per-op halo (depth). The halo equals the kernel radius:

Op	DR source	depth	boundary
slope, aspect, hillshade	dem.py (3×3 Horn/Zevenbergen)	1	"nearest"
curvature, TPI, TRI, roughness	dem.py	1 (or window radius)	"nearest"
openness, sky-view-factor	dem.py → _numba.horizon_walk_kernel	search_radius	"nearest"
color-relief	terrain.py (elementwise)	0 (map_blocks, no halo)	n/a

How. Keep the existing Numba/scipy per-block kernel; wrap it with A2 map_overlap (or
read_array(chunks=).map_overlap until A2 lands). Handle nodata by padding with the nodata sentinel / "nearest"
so edge slopes stay finite. Add an engine="tiled"/out_path= path that streams to a create_empty (A1) sink,
or return a lazy dask-backed Dataset.

# dem.py — slope, tiled path (illustrative)
def _slope_block(z):                       # the EXISTING per-array slope kernel
    ...                                    # unchanged
def slope(self, *, engine="auto", out_path=None, tile_size=512):
    if engine == "tiled":
        lazy = self.map_overlap(_slope_block, depth=1, boundary="nearest",
                                chunks=(tile_size, tile_size))
        return _write_lazy(lazy, out_path, like=self)   # A1 create_empty + A3/loop store
    return Terrain(self._slope_in_memory())             # current behaviour

Acceptance. Each op's tiled output == whole-array output bit-for-bit on fixtures, invariant to
tile_size; peak RAM ~ O(tile). Effort: S. Deps: A2 (soft — can inline map_overlap meanwhile).
References. dask array-overlap https://docs.dask.org/en/stable/array-overlap.html; gdaldem
https://gdal.org/en/stable/programs/gdaldem.html; xarray-spatial blockwise terrain kernels.

B1 — _outofcore scaffolding (tiling.py, cache.py, gid helpers)

Context. The plumbing every tiled algorithm shares: a tile plan with halos, halo-aware read/write over
pyramids, a global cell-id scheme that survives >2³¹ cells, and the Evict/Retain/Cache intermediate store (Barnes
2016 §4.2). Mostly arithmetic + thin wrappers — the algorithmic content is in B2–B4.

Notes.

• plan_tiles reuses the cloud_io.tile_windows arithmetic but separates the core window (cells a tile owns,
  written back) from the halo-expanded window (core + 1-cell ring, clipped at the domain edge). For Barnes D8
  the halo is the 1-cell perimeter; for B0 stencils it is the kernel radius.
• gid = r * full_cols + c as int64 — full_cols * full_rows ≈ 2.7e9 > 2³¹, so int32 overflows; int64 is
  what perimeter graphs key on.
• TileStore realises the three memory modes: retain = in-RAM dict; cache = spill to scratch Zarr /
  sparse GTiff (A1); evict = store nothing, recompute the per-tile pass in stage 3.

# _outofcore/tiling.py
from __future__ import annotations
from dataclasses import dataclass

@dataclass(frozen=True)
class TileSpec:
    tid: int; row: int; col: int                 # tile id + grid position
    row_off: int; col_off: int                   # CORE origin in the full raster
    n_rows: int; n_cols: int                     # CORE size (edge tiles clipped)
    halo: int = 1

    def halo_window(self, full_rows, full_cols):
        r0 = max(0, self.row_off - self.halo); c0 = max(0, self.col_off - self.halo)
        r1 = min(full_rows, self.row_off + self.n_rows + self.halo)
        c1 = min(full_cols, self.col_off + self.n_cols + self.halo)
        return (r0, c0, r1 - r0, c1 - c0)        # (row_off, col_off, n_rows, n_cols)

def plan_tiles(rows, cols, tile_rows, tile_cols, halo=1) -> list[TileSpec]: ...

def gid(r: int, c: int, full_cols: int) -> int:  # int64, overflow-safe
    return int(r) * int(full_cols) + int(c)

Acceptance. plan_tiles matches cloud_io.tile_windows for halo=0; cores tile the domain exactly (no
overlap, edge clipping); halo windows clip at the border; gid round-trips at 51810² cells. Effort: S.
Deps: A1. References. Barnes 2016 §4.2 (Evict/Cache/Retain R/W counts); cloud_io.tile_windows.

B2 — labelled priority-flood kernel (NEW Numba kernel)

Context. The one genuinely new kernel. Barnes 2016's per-tile step is the watershed variant of
Priority-Flood (Algorithm 1): it fills internal depressions and paints every cell with a watershed label,
recording for each pair of adjacent labels the lowest elevation at which they meet (the spillover graph). The
reduce step then treats that graph as a tiny DEM. Our priority_flood_numba already does the flood; it just
throws away the labels and meeting points.

Algorithm (Barnes 2016, Algorithm 1).

1. Push all edge cells (and cells adjacent to nodata) onto the min-PQ at their own elevation; all belong to
   the outlet label 1 ("ocean", drains to −∞).
2. Pop the lowest cell c. For each unvisited neighbour n: mark visited, inherit c's label, raise
   z[n] = max(z[n], z[c]), push n.
3. When a popped c touches an already-labelled neighbour n of a different label, record a spill candidate
   for (label[c], label[n]) at max(z[c], z[n]), keeping the minimum over all touching pairs.
4. A new label is minted when a popped cell is a fresh local minimum (an interior depression).

Fill values are identical to a plain Priority-Flood (epsilon=0); labels + spill graph are the additions.
Keep the existing binary-heap PQ; add a labels int32[:, :] and two growable buffers for (label_a, label_b)
and spill_z.

# _numba.py  — lazy-imported (never at module top: keeps `import digitalrivers` numba-free)
@njit(cache=True)
def priority_flood_labels_numba(z, nodata_mask, dir_dr, dir_dc):
    """Watershed Priority-Flood (Barnes 2016 Algorithm 1).
    Returns: filled_z float64[:,:] (== priority_flood_numba(z, eps=0));
             labels int32[:,:] (1 = edge/nodata outlet; >=2 = interior watersheds);
             spill_pairs int64[:,2] (label_a < label_b); spill_z float64[:] (min meet elevation per pair)."""
    ...

Acceptance. filled_z bit-for-bit == priority_flood_numba(eps=0) on every fixture; labels partition the
tile (every data cell ≥ 1); every interior depression yields ≥1 spill entry; spill_pairs symmetric-deduplicated
(a < b, min elevation kept); numba stays lazily imported. Effort: M. Deps: none (pure kernel).
References. Barnes 2016 Algorithm 1 https://arxiv.org/abs/1606.06204; reference C++
https://github.com/r-barnes/Barnes2016-ParallelPriorityFlood.

B3 — tiled fill_depressions (Barnes 2016)

Context. The orchestrator turning B2's per-tile kernel into a whole-DEM out-of-core fill: two MapReduce passes
(map-per-tile → reduce-the-edges → map-per-tile). The reduce step holds the correctness: the HANDLEEDGE join
(Algorithm 2) and the graph Priority-Flood.

The join rule (the crux). Where tile A's perimeter cell c abuts tile B's cell n (orthogonal and
diagonal), if label[c] ≠ label[n] the two watersheds spill at e = max(z[c], z[n]) — the saddle is the higher
of the two cells. Across all touching cell-pairs for a label pair, keep the minimum such e. This
max-then-min is the entire reconciliation. Worked example: min(max(4,5), max(4,1)) = min(5,4) = 4.

The global solve. The spillover graph "is itself a DEM": run an ordinary Priority-Flood on the graph (nodes =
global labels, edge weights = spill elevations), seeded from the outlet node (label 1, −∞). Result is
drain[label] = the elevation each watershed must be raised to. Reuse the same heap as B2. Broadcast drain to
the tiles and raise.

Memory mode. EVICT (default): stage 1 keeps only edge strips + the local graph; stage 3 re-reads the source
tile and recomputes the labelled flood, then raises by drain[label]. RETAIN: keep labels and skip the
recompute.

# _outofcore/fill.py
def fill_depressions_tiled(dem, out_path, *, tile_rows=2048, tile_cols=2048,
                           epsilon=0.0, cache="evict", workers=1, scratch_dir=None):
    if epsilon != 0.0:
        raise NotImplementedError("tiled fill is seam-correct for epsilon=0 only (see B6)")
    R, C = dem.rows, dem.columns
    specs = plan_tiles(R, C, tile_rows, tile_cols, halo=1)
    out = Dataset.create_empty(R, C, dtype="float32", geo=dem.geotransform,
                               epsg=dem.epsg, no_data_value=dem.band_no_data(), path=out_path)
    graph, store, edges = GlobalSpillGraph(), TileStore(cache, scratch_dir), {}
    # stage 1: map (per-tile labelled flood)
    for s in specs:
        _, core, z = read_tile(dem, s, R, C)
        fz, lab, pairs, sz = priority_flood_labels_numba(z, _nodata(z), _DIR_DR_I32, _DIR_DC_I32)
        lab = relabel_global(lab, s)
        graph.add_tile_edges(pairs, sz, s)
        edges[s.tid] = extract_perimeter(lab, fz, core)
        store.put(s.tid, labels=lab, filled=fz) if cache != "evict" else None
    # stage 2: reduce (stitch + global solve)
    for a, b in adjacent_pairs(specs):
        graph.join_border(edges[a.tid], edges[b.tid])     # max-of-two, keep-min
    drain = graph.solve()                                  # priority-flood on the label graph
    # stage 3: map (raise + write)
    for s in specs:
        if cache == "evict":
            _, core, z = read_tile(dem, s, R, C)
            fz, lab, *_ = priority_flood_labels_numba(z, _nodata(z), _DIR_DR_I32, _DIR_DC_I32)
            lab = relabel_global(lab, s)
        else:
            fz, lab = store.get(s.tid)["filled"], store.get(s.tid)["labels"]
        raised = np.maximum(fz, apply_drain(lab, drain))[core]
        out.write_array(raised, window=(s.row_off, s.col_off, s.n_rows, s.n_cols))
    return DEM(out.raster)

Acceptance (headline). For epsilon=0, tiled output bit-for-bit == in-memory fill_depressions on every
fixture and at several tile_size values (incl. a tile size that splits a depression); a depression straddling a
seam fills to one consistent level; peak RAM ~ O(tile + perimeter graph); epsilon>0 raises. Effort: L.
Deps: A1, B1, B2. References. Barnes 2016 §3 + Algorithm 1/2 + Fig. 4/5
https://arxiv.org/abs/1606.06204; C++/MPI https://github.com/r-barnes/Barnes2016-ParallelPriorityFlood.

B4 — tiled flow_accumulation (Barnes 2017, D8/ρ8)

Context. Same producer/consumer skeleton as B3, but the boundary object is a perimeter flow graph (Barnes
2017 §3). Flow directions must already exist (this op does not fill or route flats — that's B3's job). The
per-tile engine is our existing Kahn accumulation, run with inflow = 0; the new piece is the perimeter trace
and offset propagation. D8/ρ8 only — divergent flow has no fixed-halo closure (§2.4); raise on MFD/D∞.

Per-tile trace (FOLLOWPATH, Barnes 2017 Algorithm 2). After local Kahn accumulation, for each perimeter cell
c follow its flow path inside the tile until it exits an edge or dies. Store the link L(c): exit_gid
(through-router), FLOWTERMINATES (interior sink/leaf), or FLOWEXTERNAL (F(c) immediately leaves the tile —
the origin of cross-tile flow). Send (gid, flow_dir, partial_acc, L) per perimeter cell — strictly
perimeter-proportional.

Global solve. Stitch all perimeters into one flow graph via the L links + edge/corner adjacency. Reset
non-FLOWEXTERNAL cells' accumulation to 0 (interior contribution already counted locally — avoids double
counting), then run the same Kahn accumulation on the graph into an offset array A'. A'[gid] is the flow
arriving at each tile inlet from outside.

Finalize. Re-stream each tile; for each perimeter inlet, add its A' offset to every downstream cell of the
path it feeds (a second cheap Kahn pass seeded with non-zero boundary inflow), then write_array(window=).

# _outofcore/accumulate.py
def flow_accumulation_tiled(fdir, out_path, *, weights=None, tile_rows=2048,
                            tile_cols=2048, cache="evict", workers=1, scratch_dir=None):
    if fdir.routing not in ("d8", "rho8"):
        raise NotImplementedError("tiled accumulation is D8/rho8 only (see §2.4)")
    R, C = fdir.rows, fdir.columns
    specs = plan_tiles(R, C, tile_rows, tile_cols, halo=1)
    out = Dataset.create_empty(R, C, dtype="float32", geo=fdir.geotransform,
                               epsg=fdir.epsg, no_data_value=-1.0, path=out_path)
    pgraph, store = PerimeterFlowGraph(), TileStore(cache, scratch_dir)
    # stage 1: local accumulation (inflow=0) + perimeter follow-path
    for s in specs:
        _, core, fd = read_tile(fdir, s, R, C)
        w = read_tile(weights, s, R, C)[2] if weights else None
        acc = kahn_accumulate_d8_numba(fd, _w(w, fd), _DIR_DR_I32, _DIR_DC_I32)
        recs = follow_path(fd, acc, s, _DIR_DR_I32, _DIR_DC_I32)
        pgraph.add_tile_perimeter(recs, s)
        store.put(s.tid, acc=acc) if cache != "evict" else None
    # stage 2: perimeter-only global solve -> per-inlet offsets
    offsets = pgraph.solve_offsets()
    # stage 3: add offsets down each inlet's path, write
    for s in specs:
        acc = store.get(s.tid)["acc"] if cache != "evict" else _recompute(fdir, s, weights, R, C)
        acc = propagate_offsets(acc, offsets, s, _DIR_DR_I32, _DIR_DC_I32)
        out.write_array(acc[core_slice(s)], window=(s.row_off, s.col_off, s.n_rows, s.n_cols))
    return Accumulation.from_dataset(out, routing=fdir.routing)

Acceptance. D8/ρ8 tiled accumulation bit-for-bit == in-memory accumulate on fixtures and across
tile_size; weighted accumulation matches; a river straddling a seam accumulates the correct downstream total;
MFD/D∞ raises. Effort: L. Deps: A1, B1. References. Barnes 2017 §3 + Algorithms 1/2 + Fig. 1
https://arxiv.org/abs/1608.04431; C++/MPI https://github.com/r-barnes/Barnes2016-ParallelFlowAccum.

B5 — public API integration

Context. Surface the tiled engines without breaking the in-memory API. A single engine= keyword plus an
auto heuristic that flips to tiled only when the DEM would blow RAM. engine="tiled" requires out_path
(results stream to disk). All the unsupported-scope guardrails (§2.3 ε>0, §2.4 divergent) live here so misuse
fails loudly, never silently.

Heuristic. Estimate per-stage peak as ≈ k · rows · cols · 8 bytes (the surveyed ~5–10× multiplier; k≈8
fill, k≈6 accumulation), compare to available RAM (psutil.virtual_memory().available, or a configurable
absolute cell threshold when psutil is absent). Pick tiled past a safety fraction (e.g. 0.5) of available RAM.

class DEM(Dataset):
    def fill_depressions(self, *, method="priority_flood", epsilon=0.0,
                         engine="auto", out_path=None, tile_size=2048,
                         cache="evict", workers=1, **kw):
        engine = _resolve_engine(engine, self.rows, self.columns, k=8)
        if engine == "tiled":
            if out_path is None:
                raise ValueError("engine='tiled' requires out_path (streams to disk)")
            from digitalrivers._outofcore.fill import fill_depressions_tiled
            return fill_depressions_tiled(self, out_path, epsilon=epsilon,
                                          tile_rows=_t(tile_size), tile_cols=_t(tile_size),
                                          cache=cache, workers=workers)
        return self._fill_in_memory(method=method, epsilon=epsilon, **kw)

FlowDirection.accumulate / DEM.flow_accumulation get the same switch. Acceptance. auto-heuristic picks
tiled past the threshold and in-memory below it; small DEMs are byte-identical to today; engine="tiled" without
out_path, with epsilon>0, or with divergent routing each raise a clear error. Effort: S. Deps: B3, B4.
References. psutil.virtual_memory https://psutil.readthedocs.io/.

B6 — cross-seam flat & ε-fill resolution (STRETCH / research)

Context. §2.3's restriction: epsilon=0 plateau-fill composes across seams, but monotonic ε-fill and
flat-direction resolution (Garbrecht–Martz / Barnes 2014a) do not — their gradients are computed across the
whole flat, which can span tiles, so naive tiling silently diverges. This is the single biggest correctness trap
and why B3 hard-errors on epsilon>0.

Options (rough order of effort).

1. Two-phase flats: run B3 (ε=0) for globally-consistent fill, then resolve flat directions in a separate
   tiled pass over the merged surface whose halo covers each flat — bounded only when flats are small.
2. Flat-graph reduction: treat each connected flat as a node (like the spillover graph) and assign
   away-from-higher / toward-lower gradients on the merged flat graph — a second MapReduce.
3. Defer ε>0 to in-memory only indefinitely if real workloads don't need monotonic fill at >RAM scale.

Effort: XL. Deps: B3. References. Barnes, Lehman & Mulla (2014) flat-drainage; Garbrecht & Martz
(1997). Barnes & Callaghan (2021) Fill-Spill-Merge https://esurf.copernicus.org/articles/9/105/2021/ (only if
partial-fill/lake routing is added later).

B7 — dask-distributed execution (LATER)

Context. B3/B4 are written as a serial consumer loop precisely so the only change here is execution: the
stage-1 and stage-3 per-tile maps are embarrassingly parallel; the stage-2 graph reduce stays a cheap serial step
on the producer. This is the dask blockwise → tree-reduce → blockwise shape.

How. Wrap each per-tile call as a dask.delayed/future submitted to the client; gather edge data on the
producer for the reduce; broadcast the solved drain/offsets (small) back and submit the finalize tasks. Use
pyramids.configure(client=...) to replay GDAL/cloud env on every worker and default_lock() for the GeoTIFF
write (or write per-worker tiles to a Zarr/sparse-GTiff sink, A3, to avoid the single-writer lock). Workers
re-open the source by path via pyramids' CachingFileManager (handles are never serialized).

import dask
futs = [dask.delayed(_consume_tile)(path, s) for s in specs]   # stage 1, parallel
edge_data = dask.compute(*futs)
drain = graph_from(edge_data).solve()                           # stage 2, producer
dask.compute(*[dask.delayed(_finalize_tile)(path, out_path, s, drain) for s in specs])

Acceptance. Distributed results bit-for-bit == the serial path; near-linear speedup to a modest worker count;
no GDAL handle is pickled. Effort: M. Deps: B3, B4, A3 (nice). References. dask delayed/futures
https://docs.dask.org/en/stable/delayed.html; pyramids configure(client=) / default_lock().

---

5. Correctness strategy (non-negotiable)

The headline guarantee is tiled == in-memory, bit-for-bit (for the supported scope: ε=0 fill, D8/ρ8 accum):

1. Equivalence harness: every existing fixture run through both engines; assert exact array equality. Primary
   regression gate; pins seam reconciliation.
2. Seam tests: synthetic DEMs with a depression / stream straddling a tile boundary, swept over several
   tile_size values — result invariant to tiling.
3. Outlet/nodata tests: all-nodata tile contributes nothing; domain-edge + nodata-adjacent cells wire to the
   outlet node; islands handled.
4. Scale smoke test: a multi-GB synthetic DEM under a capped RAM budget (container/ulimit) proving peak
   stays O(tile).
5. Guardrails: engine="tiled" raises on epsilon>0 (pre-B6) and on divergent routing — never silently
   produce a wrong-but-plausible result.

6. Phasing & sequencing

Phase 0 (quick win):          B0  (local ops blockwise)       ── depends soft on A2
Phase 1 (pyramids primitive): A1  (create_empty)              ── unblocks all B output  [pyramids#470]
Phase 2 (fill):               B1 → B2 → B3                    ── needs A1
Phase 3 (accumulation):       B4                              ── needs A1, B1
Phase 4 (API + ergonomics):   B5                              ── needs B3, B4
Phase 5 (scale-out):          B7  (dask)                      ── needs B3/B4; A3 nice
Phase 6 (research):           B6  (seam-correct flats/ε)

Critical path: A1 → B1 → B2 → B3 → B5. B0 and A2 proceed in parallel immediately. B4 parallels B3 once B1
exists.

7. Risks & open questions

Risk	Impact	Mitigation
ε>0 / flats not seam-correct (§2.3)	tiled silently ≠ in-memory	Restrict to ε=0; hard error otherwise; B6 later. Biggest trap.
Tiled D∞/MFD accumulation unsolved (§2.4)	no divergent out-of-core	D8/ρ8 only; in-memory for divergent; error on misuse
eager map_blocks output is MEM (in-RAM)	false sense of out-of-core	use A1 + write_array(window=), not map_blocks, for the sink
New labelled kernel correctness (B2)	wrong fills	bit-for-bit equivalence harness vs existing kernel
Numba stays lazy-imported	import digitalrivers pulls numba	keep all _numba imports inline (existing rule)
Scratch-disk volume (CACHE mode)	disk pressure	EVICT default; document CACHE disk needs (~2× intermediate)
GeoTIFF write serialization (no A3 yet)	slow finalize under dask	windowed writes are independent regions; or stage to Zarr then translate

Open questions. (1) does the current priority_flood_numba already track any region labels we can extend, or
is B2 fully new? (2) exact read_array(window=) kwarg name / return shape on the current pyramids engine. (3)
threshold policy for engine="auto" (fraction of free RAM vs absolute cell count).

---

Summary table

#	Priority	Title	Depends on	Effort
B0	quick win	Local terrain ops → blockwise tiling	A2 (soft)	S
B1	required	_outofcore scaffolding (tiling/cache/gid)	A1	S
B2	required	Labelled priority-flood Numba kernel	—	M
B3	required	Tiled fill_depressions (Barnes 2016, ε=0)	A1, B1, B2	L
B4	required	Tiled flow_accumulation (Barnes 2017, D8/ρ8)	A1, B1	L
B5	required	Public engine switch + auto-heuristic + guardrails	B3, B4	S
B6	research	Cross-seam flat / ε-fill resolution	B3	XL
B7	later	Dask-distributed execution	B3, B4, A3 (opt)	M

Upstream dependency: the A1–A4 pyramids primitives are tracked in serapeum-org/pyramids#470.