How do I measure the right halo width instead of guessing?

Label the depressions on a sample of representative tiles, measure each one's bounding-box extent in cells, and take a high percentile — the 99th is a reasonable target. The halo must exceed that so no depression straddling a tile boundary is filled to two different levels. On a 1 metre LiDAR grid in low relief the 99th percentile is often 60 to 120 cells; in steep terrain it can be under 20.

Does a wider halo ever fix a flow accumulation seam?

No. Accumulation at a cell depends on every cell upstream of it, and that set is bounded by the terrain rather than by a distance. A cell one row inside a tile boundary can drain a valley four tiles away, so no halo width contains the dependency. The fix is an ordered edge exchange: each tile reports the accumulation crossing its boundaries and downstream tiles are reprocessed with those totals as seeds.

What is the overhead of a halo, and when does it stop being worth it?

The halo area is the difference between the padded tile and the core tile, so for a tile of n cells square with halo h the overhead is roughly 4h/n. At n equals 4096 and h equals 128 that is 12 percent of extra computation. At n equals 1024 with the same halo it is 50 percent, which is usually the point to increase the tile size rather than accept the waste.

Choosing Tile Size and Halo Width for Flow Accumulation

Both numbers are usually chosen by copying whatever the last project used, and both are measurable in a few minutes. Getting them wrong produces either a pipeline that wastes half its compute on redundant halo cells or one that leaves visible seams in the delineation. This guide covers the measurement, as part of the tiled and parallel watershed processing topic within watershed delineation and catchment synchronization.

Prerequisites

A representative sample of the DEM — three or four tiles spanning the terrain types present is enough.
scipy.ndimage for connected-component labelling of depressions.
The per-worker memory budget, stated as a number.

Core Technique: Measure the Halo from the Depressions

A halo is wide enough when no depression that straddles a tile boundary extends beyond it. Wider than that buys nothing; narrower produces two different fill levels on either side of the seam, and the flow directions disagree there.

The measurement is direct: label the closed depressions in a sample area, take each one’s extent in cells, and read off a high percentile.

Annotated Code Example

python

import logging

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

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

# Peak bytes per cell across the whole conditioning + routing chain: input
# float32, filled float32, working copy float32, direction uint8,
# accumulation float64. Measured, not assumed.
PEAK_BYTES_PER_CELL = 21


def measure_halo(
    dem_path: str,
    sample_windows: int = 4,
    sample_size: int = 2048,
    percentile: float = 99.0,
) -> dict:
    """
    Measure the halo a DEM actually needs, from its depression extents.

    Samples several windows, fills each, labels the cells that were raised,
    and reports the extent distribution of those depressions.
    """
    with rasterio.open(dem_path) as src:
        width, height = src.width, src.height
        nodata = src.nodata if src.nodata is not None else -9999.0

        # Spread the samples across the raster rather than clustering them,
        # so a single terrain type does not dominate the estimate.
        step_r = max(1, (height - sample_size) // max(1, sample_windows))
        step_c = max(1, (width - sample_size) // max(1, sample_windows))
        extents = []

        for i in range(sample_windows):
            r = min(height - sample_size, i * step_r)
            c = min(width - sample_size, i * step_c)
            if r < 0 or c < 0:
                continue
            patch = src.read(1, window=Window(c, r, sample_size, sample_size))
            patch = patch.astype(np.float32)

            grid = rd.rdarray(patch.copy(), no_data=nodata)
            rd.FillDepressions(grid, epsilon=False, in_place=True)
            raised = np.asarray(grid) - patch

            # Any cell the fill raised was inside a depression. Label the
            # connected groups and measure each one's bounding box.
            labels, n = ndimage.label(raised > 1e-6)
            if n == 0:
                continue
            for sl in ndimage.find_objects(labels):
                extents.append(max(sl[0].stop - sl[0].start,
                                   sl[1].stop - sl[1].start))
            log.info("window (%d,%d): %d depressions", r, c, n)

    if not extents:
        log.warning("No depressions found in the sample — a halo of 32 cells "
                    "is a safe floor for flow direction alone")
        return {"halo_cells": 32, "n_depressions": 0}

    arr = np.asarray(extents)
    p = float(np.percentile(arr, percentile))
    # Round up to a power of two so the halo aligns with typical block sizes.
    halo = int(2 ** np.ceil(np.log2(max(p, 8))))
    log.info("%d depressions sampled: median %d, p%.0f %d, max %d cells → halo %d",
             len(arr), int(np.median(arr)), percentile, int(p), int(arr.max()), halo)
    return {
        "halo_cells": halo,
        "n_depressions": len(arr),
        "median_extent": int(np.median(arr)),
        "percentile_extent": int(p),
        "max_extent": int(arr.max()),
    }


def choose_tile_size(memory_budget_gb: float, halo: int,
                     max_overhead: float = 0.15) -> int:
    """
    Largest tile that fits the budget, and smallest that keeps halo waste low.

    Returns a tile size in cells, rounded down to a multiple of 512 so it
    aligns with typical internal blocking.
    """
    budget_cells = (memory_budget_gb * 1e9) / PEAK_BYTES_PER_CELL
    # (tile + 2*halo)^2 must fit the budget.
    max_tile = int(np.sqrt(budget_cells)) - 2 * halo
    # Halo overhead is about 4*halo/tile; invert for the minimum tile.
    min_tile = int(4 * halo / max_overhead)

    if max_tile < min_tile:
        log.warning("Budget of %.1f GB cannot hold a tile large enough to keep "
                    "halo overhead under %.0f %% with a %d-cell halo — either "
                    "raise the budget or accept the overhead",
                    memory_budget_gb, 100 * max_overhead, halo)
        tile = max(512, (max_tile // 512) * 512)
    else:
        tile = max(512, (max_tile // 512) * 512)

    overhead = 4.0 * halo / tile
    padded = tile + 2 * halo
    peak_gb = (padded ** 2) * PEAK_BYTES_PER_CELL / 1e9
    log.info("Tile %d cells, halo %d → padded %d, peak %.2f GB per worker, "
             "halo overhead %.0f %%", tile, halo, padded, peak_gb, 100 * overhead)
    return tile


# --- Example usage ---
# m = measure_halo("state_1m_cog.tif", sample_windows=4)
# tile = choose_tile_size(memory_budget_gb=6.0, halo=m["halo_cells"])

Parameter Reference

Quantity	How to get it	Typical value
Halo	99th percentile depression extent, rounded to a power of two	32–256 cells
Tile size	From the memory budget, floor to a multiple of 512	2048–8192 cells
Peak bytes per cell	Sum of every simultaneous array’s dtype	≈ 21
Halo overhead	4 × halo / tile	Keep under 15 %
Accumulation halo	Not applicable	Use an ordered exchange

Worked Example: Two Terrains, Two Answers

The same measurement on two datasets gives very different halos, which is the argument for measuring at all.

	Steep mountain, 1 m LiDAR	Low-relief coastal plain, 1 m LiDAR
Depressions per 2048² sample	480	14 200
Median extent	6 cells	11 cells
99th percentile extent	18 cells	96 cells
Halo chosen	32	128
Tile at 6 GB budget	8 192	7 680
Halo overhead	1.6 %	6.7 %

The mountain dataset needs almost no halo, because depressions there are small pits from vegetation noise. The coastal plain, with its wide flat depressions behind roads and levees, needs four times as much — and running the mountain’s halo on it would leave a filled surface that disagrees across every seam.

The two knobs pull against each other under a fixed memory budget: a wider halo forces a smaller core tile, which raises the overhead the halo was meant to justify. Plotting the feasible pairs makes the compromise explicit.

Gotchas and Edge Cases

A halo copied from another project. The two rows above differ by a factor of four on the same sensor and resolution. Terrain, not sensor, sets the number.
Measuring on one sample window. A single window over uniform terrain understates the tail. Spread the samples.
Using the maximum depression extent. The largest depression is often a genuine lake, which the halo should not attempt to contain. A high percentile plus the mosaic-level fill pass is the right combination.
Applying the halo to accumulation. No width works. The accumulation stage needs the ordered exchange described in the parent topic.
Tile size not a multiple of the internal block size. Every read then straddles blocks — see writing Cloud-Optimized GeoTIFFs.
Budget stated per machine. With eight workers a 32 GB machine gives each worker 4 GB, not 32. The tile calculation takes the per-worker figure.
Ignoring the epsilon flag when measuring. Filling with epsilon=True raises a thin gradient across every flat, so almost every cell registers as raised and the depression labels merge into one enormous blob. Measure with epsilon=False.

Tiled and Parallel Watershed Processing — the parent topic and the ordered exchange that accumulation requires
Processing Large DEMs with Dask and rioxarray — where these numbers become the chunks and depth arguments
Writing Cloud-Optimized GeoTIFFs for Hydrology Pipelines — the internal blocking the tile size should align with
DEM Pit Filling Algorithms — the operation whose reach the halo is sized against

Related topics