Tiled Inference with Overlapping Windows

To run a convolutional model over a scene without tile seams, read a window that is larger than the window you write:

from rasterio.windows import Window

HALO, SIZE = 64, 256
write = Window(col, row, SIZE, SIZE)
read = Window(col - HALO, row - HALO, SIZE + 2 * HALO, SIZE + 2 * HALO)

pred = model(src.read(window=read, boundless=True, fill_value=0))
dst.write(pred[HALO:-HALO, HALO:-HALO], 1, window=write)

Every output pixel was predicted with a full neighbourhood around it, so nothing on a tile boundary is different from anything else. This page expands on running model inference over large rasters within Raster Machine Learning & Model Inference.

Why the Halo Exists

A convolutional network computes each output pixel from a neighbourhood — its receptive field. At the edge of an input tile that neighbourhood is missing, so the framework pads with zeros and the model predicts from data that is half fabricated. The result is a band of unreliable predictions around every tile, which becomes a visible grid across the mosaic.

Prediction reliability across a tile Across a tile the prediction quality is flat through the centre and falls away within roughly half a receptive field of each edge, where the convolutions see zero padding rather than real imagery. The halo places those degraded bands outside the region that is written, so only the flat central part of the profile reaches the output. Quality falls within half a receptive field of any edge halo halo written region — full context everywhere 256 px The red curve is prediction quality; the dashed lines are where the halo is cropped away. flat through here

The halo is pure overhead — those pixels are read, predicted and thrown away — but the overhead is bounded and predictable. With a 256-pixel tile and a 64-pixel halo the read is 384², which is 2.25 times the write area. Doubling the tile to 512 brings that ratio down to 1.56 for the same halo, which is why larger tiles are the cheaper lever.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Windows, boundless reads, block-aligned writes
numpy >=1.23 Tile arrays and the crop arithmetic
torch >=2.1 The segmentation model (any framework works identically) 
pip install "rasterio>=1.3.0" "numpy>=1.23" "torch>=2.1"

Complete Working Example

from dataclasses import dataclass

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


@dataclass(frozen=True)
class HaloTile:
    read: Window
    write: Window
    crop: tuple[int, int, int, int]      # top, bottom, left, right


def halo_tiles(width: int, height: int, size: int, halo: int):
    """Write windows partition the scene; read windows carry a clipped halo."""
    for row in range(0, height, size):
        for col in range(0, width, size):
            w, h = min(size, width - col), min(size, height - row)
            r0, c0 = max(0, row - halo), max(0, col - halo)
            r1, c1 = min(height, row + h + halo), min(width, col + w + halo)
            yield HaloTile(
                read=Window(c0, r0, c1 - c0, r1 - r0),
                write=Window(col, row, w, h),
                crop=(row - r0, r1 - (row + h), col - c0, c1 - (col + w)),
            )


def infer_with_halo(predict, src_path: str, dst_path: str, *,
                    size: int = 256, halo: int = 64, scale: float = 10_000.0) -> None:
    """predict: (bands, H, W) float32 -> (H, W) uint8 class array."""
    with rasterio.open(src_path) as src:
        profile = src.profile | {
            "count": 1, "dtype": "uint8", "nodata": 0, "compress": "deflate",
            "predictor": 2, "tiled": True, "blockxsize": 512, "blockysize": 512,
        }
        profile.pop("photometric", None)
        fill = src.nodata if src.nodata is not None else 0

        with rasterio.open(dst_path, "w", **profile) as dst:
            for tile in halo_tiles(src.width, src.height, size, halo):
                arr = src.read(window=tile.read, boundless=True, fill_value=fill)
                x = arr.astype("float32") / scale
                pred = predict(x)

                top, bottom, left, right = tile.crop
                h, w = pred.shape
                core = pred[top:h - bottom, left:w - right]

                # No prediction where there was no observation
                nodata = (arr[0] == fill)[top:h - bottom, left:w - right]
                core = np.where(nodata, 0, core).astype("uint8")

                assert core.shape == (tile.write.height, tile.write.width)
                dst.write(core, 1, window=tile.write)

The assertion before the write is worth keeping in production. It costs nothing and it is the difference between a shifted output discovered during QA and one discovered by a user six months later.

Variant Patterns

1. Sub-tiling a large read for a fixed-input model

Some exported models accept only one input size. The read window can still be large — split it into fixed-size sub-tiles after reading, which keeps the number of file reads low while satisfying the model.

One read, four model calls A single large windowed read covers a region that is then divided into four fixed-size inputs for a model that cannot accept a variable shape. The file is touched once instead of four times, which matters most when the source is remote, and the four predictions are reassembled before the halo is cropped. Read once, predict in fixed-size pieces one read: 896 x 896 1 HTTP range set, 1 decode pass 384 x 384 384 x 384 384 x 384 384 x 384 4 model calls, batched into one device transfer

2. Halo zero for a pixel-based model

A model with a one-pixel receptive field needs no halo at all, and setting one wastes 125% of the read. Make the halo a parameter rather than a constant so the same loop serves both model classes.

halo = 0 if contract["chip_size"] == 1 else contract["chip_size"] // 4

3. Writing probabilities alongside classes

Keeping the winning probability costs one more band and makes the output re-thresholdable later, as described in exporting and serving model outputs.

import numpy as np

probs = softmax_output[:, top:h - bottom, left:w - right]
cls = probs.argmax(axis=0).astype("uint8") + 1
conf = np.clip(np.round(probs.max(axis=0) * 254) + 1, 1, 255).astype("uint8")
dst.write(cls, 1, window=tile.write)
dst.write(conf, 2, window=tile.write)

Budgeting the Overhead

The halo turns every tile into more work than the output it produces, and it is worth knowing how much before a continental run starts.

Wasted read fraction by tile size With a 64 pixel halo a 128 pixel tile reads four times the pixels it writes, a 256 pixel tile reads 2.25 times, a 512 pixel tile 1.56 times and a 1024 pixel tile 1.27 times. Halving the halo to 32 pixels brings a 256 pixel tile down to 1.56. The curve flattens beyond about 512, which is where memory rather than overhead becomes the binding constraint. Pixels read per pixel written 4.0x 1.0x halo 64 halo 32 halo 16 128 px 256 px 512 px 1024 px Tile size, at a fixed halo. The gain from going past 512 is small; the memory cost is not.

The arithmetic is ((size + 2 * halo) / size) ** 2. At 128 pixels with a 64-pixel halo you read four pixels for every one written — three quarters of the run is waste. At 512 it is 1.56, and at 1024 it is 1.27. Since inference over a whole archive is usually read-bound, that ratio is close to a direct multiplier on wall clock and on egress charges.

Memory pushes the other way. Peak residency is batch × bands × (size + 2 * halo)² × 4 bytes, plus whatever activations the model holds, which for a U-Net is several times the input. A batch of eight 384² six-band tiles is 28 MB of input and perhaps ten times that in activations — comfortable on any modern GPU. The same batch at 1152² is 250 MB of input and will not fit alongside the activations on a 16 GB card.

The practical resolution is to fix the tile at 512 with a halo of 64, giving a 640² read and a 1.56 ratio, then tune the batch size to the device. That combination is close to optimal for the common case of a four-level segmentation network over six bands, and it aligns naturally with the 512-pixel internal blocks that a well-written Cloud-Optimized GeoTIFF already uses — so every read maps onto whole blocks rather than straddling them.

Common Errors

The output is offset by exactly the halo

The crop used a constant halo instead of the per-tile values. Edge tiles have an asymmetric halo, so the crop must come from the tile, not from the configuration.

ValueError: could not broadcast input array

The prediction is a different shape from the write window, usually because the model changed the spatial dimensions — an odd input size with stride-2 downsampling will do it. Pad the read to a multiple of the model’s downsampling factor.

Predictions appear over nodata regions

The model was run on padded pixels and its output was written. Mask with the source nodata after cropping, as the example does, so the output footprint matches the input.

Frequently Asked Questions

Q: How do I choose the halo width? Take half the model’s receptive field and round up to a multiple of the downsampling factor. For a U-Net with four pooling stages the receptive field is typically around 128 pixels, so a halo of 64 is the usual answer.

Q: Is it cheaper to use bigger tiles or a smaller halo? Bigger tiles. The wasted work is the halo ring, whose area grows linearly with tile edge while the useful area grows quadratically, so doubling the tile roughly halves the overhead fraction. Increase the tile until memory becomes the constraint, then stop.

Q: What happens at the corners of the scene? There is no data to pad with, so the halo is clipped and the crop shrinks accordingly. Those output pixels genuinely have less context than the rest, which is unavoidable; boundless reads keep the array shape fixed so the model still runs.