Merging Tiles with rasterio.merge

To assemble adjacent rasters into one mosaic, pass their datasets to merge with an explicit nodata value and a deliberate ordering:

import rasterio
from rasterio.merge import merge

sources = [rasterio.open(p) for p in sorted(tile_paths)]     # order decides overlap priority
mosaic, transform = merge(sources, nodata=0, method="first")

profile = sources[0].profile | {
    "height": mosaic.shape[1], "width": mosaic.shape[2],
    "transform": transform, "driver": "COG", "compress": "DEFLATE",
}
with rasterio.open("mosaic.tif", "w", **profile) as dst:
    dst.write(mosaic)
for src in sources:
    src.close()

Merging is the assembly step that precedes the blending questions in Seamless Mosaicking and Edge Blending.

Why This Arises in Remote Sensing Workflows

Satellite data arrives tiled, and study areas rarely respect tile boundaries. A catchment straddles two Sentinel-2 tiles; a district spans three Landsat paths; a national product is assembled from hundreds of scenes. Somewhere in every pipeline, adjacent rasters have to become one.

merge does the geometric part: it computes the union extent, allocates the output array, and copies each input into its place. What it does not do is make the result look continuous. Where two scenes overlap it applies a fixed rule — first wins, last wins, maximum wins — and the boundary between them is a hard edge. If the scenes differ radiometrically, that edge is visible, and no merge method hides it.

Separating those two concerns keeps the pipeline honest. Merge assembles; blending and histogram matching make the assembly look continuous, and they are separate steps with their own decisions.

What each merge method does in the overlap With first, the earlier dataset in the list keeps the overlap. With last, the later one does. With max, the brighter pixel wins, which favours cloud. With min, the darker wins, which favours shadow. None of them blends, so the join remains a hard edge in every case. Two overlapping scenes, four rules scene A scene B overlap method overlap gets first (default) scene A — list order decides last scene B max the brighter pixel — favours cloud min the darker pixel — favours shadow A custom callable can do better: mean over the overlap, or a weighted blend by distance to the edge. Whatever the rule, merge produces a hard boundary — feathering is a separate operation.

Environment & Setup

Package Version Why
rasterio ≥1.3.0 merge, dataset handling, COG output
numpy ≥1.23 Custom merge callables
shapely ≥2.0 Optional: footprint checks before merging 
pip install "rasterio>=1.3.0" "numpy>=1.23"

Complete Working Example

This function validates that the inputs can legitimately be merged, orders them by a quality key, merges with an explicit nodata, and writes a tiled output with overviews.

Compatibility requirements merge does not check merge assumes a shared CRS, resolution and dtype and produces a plausible-looking but wrong result when they differ. Checking them before the call converts three silent failures into one clear error. What to assert before calling merge property if it differs symptom in the output CRS not checked by merge one scene lands in the wrong place resolution not checked a stretched or squashed region dtype not checked silent truncation of values nodata declaration defaults to none fill pixels win the overlap band count raises the only one that fails loudly Four of these five fail silently; the wrapper in this page turns them into exceptions. 
import rasterio
from rasterio.enums import Resampling
from rasterio.merge import merge


def merge_tiles(
    paths: list[str],
    dst_path: str,
    *,
    nodata: float | int = 0,
    method: str = "first",
    quality_key=None,
) -> dict:
    """Merge compatible rasters into one mosaic, with overlap priority under your control."""
    if quality_key is not None:
        paths = sorted(paths, key=quality_key)      # best first: 'first' then keeps the best

    sources = [rasterio.open(p) for p in paths]
    try:
        crs = {src.crs for src in sources}
        res = {tuple(round(v, 6) for v in src.res) for src in sources}
        dtypes = {src.dtypes[0] for src in sources}
        if len(crs) != 1:
            raise ValueError(f"inputs span {len(crs)} CRSs — reproject onto one grid first")
        if len(res) != 1:
            raise ValueError(f"inputs have {len(res)} resolutions: {res}")
        if len(dtypes) != 1:
            raise ValueError(f"inputs have mixed dtypes: {dtypes}")

        mosaic, transform = merge(sources, nodata=nodata, method=method)

        profile = sources[0].profile | {
            "driver": "COG",
            "height": mosaic.shape[1],
            "width": mosaic.shape[2],
            "count": mosaic.shape[0],
            "transform": transform,
            "nodata": nodata,
            "compress": "DEFLATE",
            "blocksize": 512,
        }
    finally:
        for src in sources:
            src.close()

    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(mosaic)
        dst.build_overviews([2, 4, 8, 16], Resampling.average)
        dst.update_tags(ns="rio_overview", resampling="average")
        dst.update_tags(merged_from=str(len(paths)), merge_method=method)

    return {"shape": mosaic.shape, "transform": transform, "inputs": len(paths)}


if __name__ == "__main__":
    # Least cloudy first, so 'first' keeps the cleanest pixels in every overlap
    info = merge_tiles(tile_paths, "mosaic.tif", nodata=0, method="first",
                       quality_key=lambda p: cloud_fraction_lookup[p])
    print(info)

The three compatibility assertions are the point of the wrapper. merge does not check them, and each failure produces a plausible-looking output: mixed CRSs give a mosaic where one scene is in the wrong place, mixed resolutions give one that is stretched, and mixed dtypes give silent truncation.

Variant Patterns

1. A custom merge callable

When neither first nor last is right — for example, averaging the overlap rather than choosing — pass a callable.

import numpy as np


def mean_merge(merged_data, new_data, merged_mask, new_mask, **kwargs):
    """Average where both contribute; take whichever exists elsewhere."""
    both = ~merged_mask & ~new_mask
    only_new = merged_mask & ~new_mask

    merged_data[both] = (merged_data[both] + new_data[both]) / 2
    merged_data[only_new] = new_data[only_new]
    merged_mask[only_new] = False


mosaic, transform = merge(sources, nodata=0, method=mean_merge)

Averaging removes the hard edge but replaces it with a band where the values are neither scene’s, which is visible whenever the two differ radiometrically — the case that motivates Histogram Matching Across Scenes.

2. Merging into a fixed output grid

By default the mosaic covers the union of the inputs. For a tiled product you usually want a predetermined extent so outputs abut exactly.

mosaic, transform = merge(
    sources,
    bounds=(300_000, 9_880_000, 400_000, 9_980_000),   # the product tile, not the data extent
    res=(10.0, 10.0),
    nodata=0,
)

Fixing bounds and resolution makes the output reproducible and lets two runs produce files that align pixel for pixel, which is what a downstream mosaic of mosaics requires.

3. Merging more data than fits in memory

Two ways to assemble a large mosaic rasterio.merge allocates the whole output array in memory, so a national mosaic can exceed available RAM before any pixel is written. Iterating the output's blocks and merging only the inputs that intersect each block keeps peak memory at one block, at the cost of more open calls. Peak memory: whole mosaic versus one block merge(sources) entire output array in RAM 40,000 × 40,000 uint16 ≈ 3.2 GB block-wise merge one 512² block ≈ 0.5 MB at a time Same output, same values — only the allocation strategy differs. 
import rasterio
from rasterio.merge import merge
from rasterio.windows import Window

with rasterio.open(dst_path, "w", **profile) as dst:
    for _, window in dst.block_windows(1):
        bounds = rasterio.windows.bounds(window, dst.transform)
        hits = [p for p in paths if footprints[p].intersects(box(*bounds))]
        if not hits:
            continue
        srcs = [rasterio.open(p) for p in hits]
        block, _ = merge(srcs, bounds=bounds, res=dst.res, nodata=profile["nodata"])
        for s in srcs:
            s.close()
        dst.write(block, window=window)

The footprint index is what keeps this from reopening every input for every block; building it once from the STAC items costs nothing, as in Querying STAC Catalogs Programmatically.

Checking the Mosaic

Three checks catch the mistakes that survive a successful merge.

Extent and alignment: the output transform’s origin should be a multiple of the resolution, and the extent should match what you expect from the input footprints. An origin at an odd offset means the inputs were not on a shared grid.

Nodata behaviour: sample a few pixels along a join. A pixel that is valid in one input and fill in the other must take the valid value; if it is fill, the nodata declaration did not reach merge.

Continuity across joins: take a transect crossing a seam and plot the values. A step indicates a radiometric difference that merging cannot fix, and points to the histogram-matching and feathering steps in Removing Seams in Multi-Scene Mosaics with Feathering.

Common Errors

MergeError: Datasets have different CRS

Reproject onto a single grid first. Merging across CRSs is not something to work around — the result would be geometrically wrong.

The mosaic is enormous and mostly empty

One input has an outlying footprint — often a scene from the wrong tile or a mislabelled CRS — so the union extent balloons. Check src.bounds for every input before merging.

Memory blows up on a large mosaic

merge allocates the whole output. Use the block-wise pattern above, or merge in groups and then merge the groups.

Frequently Asked Questions

Q: What decides which scene wins in an overlap? The method argument. The default first keeps whichever dataset appears earliest in the list, so ordering the inputs is a real decision — usually least cloudy first, or most recent first.

Q: Why does my mosaic have black stripes along the joins? The inputs declare no nodata, so fill pixels are treated as data and win in the overlap. Pass nodata explicitly to merge, and make sure the source files declare it too.

Q: Can I merge rasters in different CRSs? No. merge assumes a common CRS and resolution, and mixing them produces a misaligned result rather than an error. Reproject onto one grid first, then merge.