Stitching Prediction Tiles without Seams

If cropping a halo still leaves a visible join, blend instead of cutting: accumulate weighted probabilities from every overlapping tile and normalise at the end.

import numpy as np

taper = np.hanning(size)[:, None] * np.hanning(size)[None, :]   # 0 at the edge, 1 in the middle
prob_acc[:, r0:r1, c0:c1] += probs * taper
weight_acc[r0:r1, c0:c1] += taper
# ... after every tile has contributed ...
classes = (prob_acc / np.maximum(weight_acc, 1e-6)).argmax(axis=0).astype("uint8") + 1

Because the taper falls to zero at each tile’s edge, the pixels a tile is least sure about contribute least, and the transition between tiles is continuous rather than abrupt. This page extends running model inference over large rasters, part of Raster Machine Learning & Model Inference.

What a Seam Actually Is

A seam is a discontinuity in the decision function, not in the imagery. Two adjacent output pixels came from different model calls; one had context to its right, the other had zero padding. Where the model was near a decision boundary, that difference flips the class, and the flip aligns perfectly with the tile grid.

Three ways to join tiles With abutting tiles the class flips along the shared edge and a hard line appears. Cropping a halo removes the line because both sides were predicted with full context. Weighted blending goes further by fading each tile's contribution to zero at its edge, so even an unstable model produces a continuous field. abutting cropped halo weighted blend hard line on the tile grid no line; both sides had context gradual, even for an unstable model Lower row: the predicted probability along a line crossing the join.

This is why the seam test in the parent page works on the output alone: real class boundaries do not know where your tile grid is, so any concentration of transitions at multiples of the tile size is an artefact by construction.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Taper construction and accumulation buffers
rasterio >=1.3.0 Streaming the output band of rows to disk
scipy >=1.11 Optional smoother window functions 
pip install "numpy>=1.23" "rasterio>=1.3.0" "scipy>=1.11"

Complete Working Example

The buffer cannot hold a whole scene, so the accumulator covers a band of tile rows and flushes the rows that no future tile will touch.

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


def cosine_taper(size: int, flat: float = 0.5) -> np.ndarray:
    """Separable window: flat through the middle, cosine fall to zero at the edge."""
    ramp = int(size * (1 - flat) / 2)
    w = np.ones(size, dtype="float32")
    edge = 0.5 * (1 - np.cos(np.linspace(0, np.pi, ramp, dtype="float32")))
    w[:ramp], w[-ramp:] = edge, edge[::-1]
    return (w[:, None] * w[None, :]).astype("float32")


def blended_inference(predict, src_path: str, dst_path: str, *,
                      size: int = 512, stride: int = 384, n_classes: int = 4) -> None:
    """Overlapping tiles blended by a cosine taper, flushed row band by row band."""
    taper = cosine_taper(size)

    with rasterio.open(src_path) as src:
        profile = src.profile | {"count": 1, "dtype": "uint8", "nodata": 0,
                                 "compress": "deflate", "tiled": True}
        profile.pop("photometric", None)

        with rasterio.open(dst_path, "w", **profile) as dst:
            for row in range(0, src.height, stride):
                band_h = min(size, src.height - row)
                acc = np.zeros((n_classes, band_h, src.width), dtype="float32")
                wsum = np.zeros((band_h, src.width), dtype="float32")

                for col in range(0, src.width, stride):
                    w = min(size, src.width - col)
                    win = Window(col, row, w, band_h)
                    arr = src.read(window=win, boundless=True,
                                   fill_value=src.nodata or 0).astype("float32")
                    probs = predict(arr / 10_000.0)          # (classes, h, w)

                    t = taper[:band_h, :w]
                    acc[:, :, col:col + w] += probs * t
                    wsum[:, col:col + w] += t

                out = (acc / np.maximum(wsum, 1e-6)).argmax(axis=0).astype("uint8") + 1
                flush_h = min(stride, band_h)                 # rows no later tile touches
                dst.write(out[:flush_h], 1,
                          window=Window(0, row, src.width, flush_h))
Tapers that sum to a constant Three overlapping tile weight profiles are drawn, each flat through its centre and falling as a cosine to zero at its edges. Where they overlap their sum stays flat, so dividing by the accumulated weight recovers an unbiased average everywhere. A taper narrower than the overlap would leave the sum dipping between tiles. Per-tile weight and the accumulated total sum tile 1 tile 2 tile 3 tile 4 Check it numerically: the interior of wsum should vary by well under five per cent.

Two subtleties are worth naming. The taper must never be exactly zero across a whole pixel column, or wsum is zero there and the normalisation divides by the epsilon instead — hence the flat centre and the clamp. And the flush height is stride, not size: rows beyond the stride will receive contributions from the next band of tiles and must not be written yet.

Choosing Stride, Taper and Buffer Size

Three parameters interact, and the sensible region of the space is small.

Stride and taper combinations A stride equal to the tile size gives no overlap, so blending is impossible and only cropping can help. A stride of three quarters of the tile costs about seventy per cent more compute and removes seams for most models. A stride of half the tile costs four times the compute and is only justified for a model that is genuinely unstable at its border. Overlap buys smoothness at a quadratic price stride = size no overlap at all 1.0x compute blending impossible use a cropped halo instead stride = 0.75 x size 25% overlap each way 1.8x compute seams gone for most models the default worth starting from stride = 0.5 x size every pixel predicted 4x 4.0x compute smoothest possible result rarely worth the bill Compute multiplier is (size / stride) squared — the cost of overlap compounds in both axes. Run the seam test at 0.75 first; only move to 0.5 if it still fails.

The taper’s flat fraction controls how much of each tile is treated as fully trusted. A flat fraction of 0.5 means the central half of the tile contributes at full weight and the outer quarter on each side ramps down — a good match for a stride of 0.75. Making the taper narrower than the overlap wastes the overlap; making it wider than the overlap leaves regions where the total weight dips, producing a soft version of the very banding you are trying to remove.

Buffer memory is the third constraint, and it is the one that bites on large scenes. A four-class float32 accumulator over a 10,980-pixel-wide band of 512 rows is 90 MB, plus the weight buffer. That is fine. The same accumulator over a whole scene would be 1.9 GB, which is why the band-and-flush structure exists rather than a simpler two-pass approach.

An alternative that avoids buffers entirely is worth knowing about: predict into separate per-tile files and mosaic them afterwards with feathering, reusing the machinery in removing seams in multi-scene mosaics with feathering. It costs disk and a second pass, but it parallelises across machines trivially, which the in-memory accumulator does not.

Common Errors

A faint grid appears at the stride spacing rather than the tile spacing

The taper and the stride disagree, so total weight dips periodically. Check that wsum is roughly constant across the interior: wsum.std() / wsum.mean() should be well under 0.05.

MemoryError allocating the accumulator

The buffer was sized to the scene rather than to a band of rows. Accumulate over size rows, flush stride rows, and move on.

The blended output is blurrier than the cropped one

Blending averages probabilities, which softens boundaries by design. If crisp edges matter more than continuity, use a cropped halo and leave the probabilities alone.

Frequently Asked Questions

Q: When is cropping enough and blending unnecessary? Whenever the halo is at least half the receptive field, cropping alone removes the seam completely — every written pixel had full context. Blending matters when the halo has to be small for memory reasons, or when the model is genuinely unstable near its input border.

Q: Can I blend class labels instead of probabilities? No. Averaging class codes is meaningless — the mean of forest and water is not a class. Blend the probability or logit tensor and take the argmax once at the end, or use a per-pixel vote if probabilities are unavailable.

Q: How much memory does an accumulation buffer need? Classes times rows times columns times four bytes for a float32 buffer, which is far too much for a whole scene. Accumulate over a band of tile rows, flush completed rows to disk, and keep only the rows still being written to.