Tiled and Parallel Watershed Processing

A statewide 1 m DEM is a few terabytes, and no amount of patience makes it fit in memory. As part of the watershed delineation and catchment synchronization workflow, tiling is what makes that dataset tractable — and it is also the fastest way to produce a delineation that is subtly wrong everywhere, because the operations in a hydrologic pipeline do not all tile the same way. Some are local and parallelise trivially. One of them, flow accumulation, has a dependency structure that no fixed overlap can contain, and treating it like the others is the defining error of large-scale hydrologic processing.

This page is about that distinction and the machinery it forces: how to size tiles, which stages take a halo, why accumulation needs an ordered two-pass exchange instead, and how to prove afterwards that no seam survived.

Prerequisites and Environment Setup

conda create -n tiledhydro python=3.11
conda activate tiledhydro
conda install -c conda-forge rioxarray=0.15 dask=2024.2 distributed \
    rasterio=1.3 richdem=0.3 geopandas=0.14 zarr numcodecs

The most common environment mistake is running a Dask cluster with more workers than the memory budget supports. Each worker holds a tile plus its halo at several dtypes simultaneously — elevation as float32, the filled surface as float32, directions as uint8, accumulation as float64 — and the sum, not the input size, is what has to fit.

Mechanics: Which Stages Tile, and Which Do Not

The stages of a delineation pipeline fall into three groups with very different parallelisation properties.

The first group needs nothing. The second group needs a halo wide enough to contain the operation’s reach: for depression filling that means wide enough to contain the largest depression that straddles the boundary, which is usually a few dozen cells but should be measured rather than assumed. The third group is different in kind.

Why accumulation cannot use a halo

Flow accumulation at a cell is the count of every cell upstream of it. On a tiled grid, “upstream” can cross an arbitrary number of tiles: a cell on the downstream edge of a tile may drain a valley that begins four tiles away. There is no halo width that contains that, because the dependency is a property of the terrain rather than of the geometry.

The correct treatment is an ordered exchange. Each tile computes a local accumulation as though nothing drains into it. Each tile then reports, for every boundary cell where flow leaves, how much accumulation crosses. Those reports define a dependency graph over tiles. Tiles are then reprocessed in topological order, each one seeded with the totals arriving at its inbound edges.

Sizing the tile

Tile size is a memory calculation, not a preference. For a tile of n × n cells with a halo of h:

A 4096 × 4096 tile with a 128-cell halo is (4096 + 256)² ≈ 18.9 million cells, so about 400 MB at peak. Four workers on a 4 GB budget is comfortable; sixteen workers is not. Tiles smaller than about 1024 cells rarely pay off, because the halo starts to dominate: at n = 1024 and h = 128, the halo is 28 % of the processed area and every cell in it is computed twice.

Step-by-Step Workflow

1. Choose the tile size from the memory budget above, then round to a multiple of the file’s internal block size so reads are aligned.
2. Condition each tile with its halo. Read the tile plus margin, fill depressions, discard the margin. This is a map_overlap and parallelises cleanly.
3. Compute flow direction per tile, also with a halo, also discarding the margin.
4. Compute a local accumulation per tile with no inbound seeds. Record, for every boundary cell whose flow leaves the tile, the accumulation value and the neighbour it enters.
5. Build the tile dependency graph from those records and sort it topologically. A tile is ready when every tile that drains into it has been finalised.
6. Re-run accumulation in dependency order, seeding each tile’s inbound boundary cells with the totals reported by its upstream neighbours.
7. Delineate and label on the finalised accumulation, then validate the seams.

Production-Ready Code

The driver below runs the local pass across tiles in parallel, builds the dependency graph from the edge reports, and re-runs the tiles in topological order. It uses rioxarray for windowed I/O and richdem for the per-tile hydrology.

import logging
from collections import defaultdict
from concurrent.futures import ProcessPoolExecutor

import numpy as np
import rasterio
import richdem as rd
from rasterio.windows import Window

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


def tile_windows(width: int, height: int, tile: int, halo: int):
    """Yield (core_window, read_window, halo_offsets) for every tile."""
    for row in range(0, height, tile):
        for col in range(0, width, tile):
            core = Window(col, row, min(tile, width - col), min(tile, height - row))
            r0 = max(0, row - halo)
            c0 = max(0, col - halo)
            r1 = min(height, row + core.height + halo)
            c1 = min(width, col + core.width + halo)
            read = Window(c0, r0, c1 - c0, r1 - r0)
            yield core, read, (row - r0, col - c0)


def process_tile_local(args) -> dict:
    """Condition and route one tile with its halo; return edge exports.

    Runs in a worker process, so it takes and returns only picklable data.
    """
    dem_path, core, read, offsets, nodata = args
    with rasterio.open(dem_path) as src:
        patch = src.read(1, window=read).astype(np.float32)

    dem = rd.rdarray(patch, no_data=nodata)
    rd.FillDepressions(dem, epsilon=True, in_place=True)
    fdir = rd.FlowProportions(dem, method="D8")
    acc = rd.FlowAccumulation(dem, method="D8")

    r_off, c_off = offsets
    h, w = int(core.height), int(core.width)
    core_acc = np.asarray(acc)[r_off:r_off + h, c_off:c_off + w]
    core_dir = np.asarray(fdir)[r_off:r_off + h, c_off:c_off + w]

    # Report accumulation crossing each edge, keyed by the neighbour it enters.
    exports = defaultdict(float)
    exports["N"] = float(core_acc[0, :].sum())
    exports["S"] = float(core_acc[-1, :].sum())
    exports["W"] = float(core_acc[:, 0].sum())
    exports["E"] = float(core_acc[:, -1].sum())

    return {
        "key": (int(core.row_off), int(core.col_off)),
        "local_total": float(core_acc.max()),
        "exports": dict(exports),
        "acc": core_acc,
        "fdir": core_dir,
    }


def run_tiled_accumulation(
    dem_path: str,
    out_path: str,
    tile: int = 4096,
    halo: int = 128,
    workers: int = 4,
) -> dict:
    """
    Two-pass tiled flow accumulation: local pass in parallel, then an ordered
    pass that seeds each tile with the totals arriving at its edges.
    """
    with rasterio.open(dem_path) as src:
        width, height = src.width, src.height
        profile = src.profile.copy()
        nodata = src.nodata if src.nodata is not None else -9999.0

    jobs = [(dem_path, core, read, off, nodata)
            for core, read, off in tile_windows(width, height, tile, halo)]
    log.info("Processing %d tiles of %d px with a %d px halo", len(jobs), tile, halo)

    # --- Pass 1: every tile independently, in parallel ---
    with ProcessPoolExecutor(max_workers=workers) as pool:
        results = list(pool.map(process_tile_local, jobs))
    log.info("Local pass complete: %d tiles", len(results))

    by_key = {r["key"]: r for r in results}

    # --- Build the dependency graph from the edge exports ---
    # A tile depends on any neighbour that exports accumulation into it.
    deps = defaultdict(set)
    for (row, col), r in by_key.items():
        for side, amount in r["exports"].items():
            if amount <= 0:
                continue
            nbr = {
                "N": (row - tile, col), "S": (row + tile, col),
                "W": (row, col - tile), "E": (row, col + tile),
            }[side]
            if nbr in by_key:
                deps[nbr].add((row, col))

    # --- Topological order: a tile is ready when all its feeders are done ---
    order, done = [], set()
    remaining = set(by_key)
    while remaining:
        ready = [k for k in remaining if deps[k] <= done]
        if not ready:
            # Only possible if the direction grids disagree across a seam.
            log.error("Tile dependency cycle across %d tiles — seams disagree; "
                      "re-condition with a wider halo", len(remaining))
            ready = sorted(remaining)
        for k in sorted(ready):
            order.append(k)
            done.add(k)
            remaining.discard(k)
    log.info("Resolved a processing order over %d tiles, max depth %d",
             len(order), len({len(deps[k]) for k in order}))

    # --- Pass 2: finalise in order, adding inbound totals as seeds ---
    profile.update(dtype="float64", count=1, compress="lzw", tiled=True, nodata=0)
    with rasterio.open(out_path, "w", **profile) as dst:
        for key in order:
            r = by_key[key]
            inbound = sum(
                by_key[f]["exports"].get(_side_between(f, key, tile), 0.0)
                for f in deps[key]
            )
            final = r["acc"] + inbound
            row, col = key
            dst.write(final.astype("float64"), 1,
                      window=Window(col, row, final.shape[1], final.shape[0]))
            log.info("tile (%d,%d): local max %.0f + inbound %.0f",
                     row, col, r["local_total"], inbound)

    log.info("Wrote tiled accumulation: %s", out_path)
    return {"tiles": len(order), "output": out_path}


def _side_between(feeder, receiver, tile):
    """Which edge of the feeder tile touches the receiver tile."""
    fr, fc = feeder
    rr, rc = receiver
    if rr == fr + tile:
        return "S"
    if rr == fr - tile:
        return "N"
    if rc == fc + tile:
        return "E"
    return "W"


# --- Example usage ---
# run_tiled_accumulation(
#     dem_path="state_1m_cog.tif",
#     out_path="state_acc.tif",
#     tile=4096, halo=128, workers=6,
# )

Validation Protocol

Seam artefacts are the failure mode this whole design exists to prevent, and they are easy to test for.

• Accumulation continuity across seams. Sample the accumulation raster in a one-cell strip either side of every tile boundary. The distribution of the ratio should be centred on 1.0 with a narrow spread. A cluster of ratios near zero on one side is an unseeded tile.
• Whole-domain total. The maximum accumulation at the domain outlet should equal the total number of contributing cells. Tiled runs that dropped an exchange come out low, and the shortfall is exactly the missing tiles’ area.
• Catchment boundaries against tile boundaries. Intersect the delineated catchment edges with the tile grid. A real divide crosses a tile boundary at an arbitrary angle; an artefact runs along it. Any divide segment more than a few hundred metres long that is collinear with a seam is a defect.
• Single-tile reference. Process one moderate sub-basin both ways — tiled and as a single window — and difference the results. They should agree cell for cell. This is the only check that tests the whole chain at once, and it is worth the runtime.

Speedup is bounded by the dependency graph long before it is bounded by hardware. Measuring where the time actually goes on a real run explains why adding workers stops helping.

Common Failure Modes and Optimization

• map_overlap for accumulation. Produces a raster where every tile’s accumulation restarts near zero at its upstream edge. The numbers look plausible in isolation, which is what makes this dangerous.
• Per-tile conditioning without a halo. The filled surfaces disagree at the seam by a few centimetres, the directions disagree, and catchment boundaries follow the seam. Symptom: suspiciously straight divides on a north–south grid.
• Tiles cut across major valleys. Maximises the dependency depth and destroys concurrency. Where the basin structure is known, align tiles to sub-basin boundaries — the delineation then parallelises almost perfectly.
• Halo too narrow for the largest depression. A depression straddling a boundary is filled to different levels on each side. Measure the largest depression’s extent on a sample tile and set the halo above it.
• Peak memory measured on the input. The input is float32; the peak is roughly five times that. Size from the table above, not from the file size.
• Object storage without range reads. Reading a window from a non-COG GeoTIFF on object storage transfers the whole file. Convert to COG first; the conversion pays for itself on the second tile.

When to Use This vs. Alternatives

Process in one window whenever the DEM fits, with room for the working copies. A 10 m DEM of a 5 000 km² basin is 50 million cells — about a gigabyte at peak — and tiling it adds complexity and seam risk for no benefit.

Split by basin instead of by grid when the basins are known and independent. Delineating each HUC separately is embarrassingly parallel with no exchange at all, because no water crosses a basin divide. This is almost always better than grid tiling when it is available, and it is the pattern the pipeline orchestration page builds its DAGs around.

Coarsen first when the analysis does not need the resolution. Resampling a 1 m grid to 10 m before delineation reduces the work a hundredfold and, for regional questions, changes the answer by a fraction of a percent — see choosing between 10 m and 1 m DEM resolution. Tile only when the fine resolution is genuinely required over the full extent.

Frequently Asked Questions

How wide should the halo be for tiled flow accumulation?

For depression filling and flow direction, a halo of a few dozen cells is usually enough, because both are local operations whose dependency extends only as far as the depression being filled. For flow accumulation there is no safe fixed halo: a single cell can drain a strip that runs the full width of the tile, so the dependency is topological rather than spatial. Accumulation needs an edge-exchange pass, not a wider margin.

Can Dask parallelise watershed delineation directly?

Dask parallelises the embarrassingly parallel parts well — reading, reprojection, per-tile conditioning, per-tile flow direction, and any per-cell arithmetic. It cannot express flow accumulation as a simple map_overlap, because accumulation has an unbounded, data-dependent dependency across tiles. The usual pattern is Dask for everything local and a separate ordered pass for accumulation and delineation.

Why does my tiled delineation show straight lines along tile boundaries?

Because each tile was conditioned in isolation, so the filled surface differs slightly on either side of the seam and the flow directions disagree there. The catchment boundary then follows the seam, because that is where the two surfaces meet. Conditioning with an overlap margin and discarding the margin afterwards removes the artefact; conditioning the mosaic removes it entirely.

Related Topics

• Watershed Delineation & Catchment Synchronization — the parent section: outlets, partitioning, topology and the synchronisation rules tiling must not break
• Pipeline Orchestration: DEM to Watershed — the scheduler layer that runs these tiles as a managed DAG
• Basin Partitioning Strategies — splitting by basin instead of by grid, and when that is the better parallelisation
• Boundary Topology Validation — the checks that catch a seam artefact once catchments are polygons
• Spatial Resolution Tradeoffs — the coarsening option that removes the need to tile at all