How to Merge and Mosaic Rasters in Python with Rasterio #

Problem statement #

National raster data arrives as tiles. You download 240 GeoTIFFs named SJ80NE.tif, SJ80NW.tif and so on, and every analysis you want to run needs one continuous surface.

The one-liner exists:

from rasterio.merge import merge
mosaic, transform = merge([rasterio.open(p) for p in paths])

and it fails in four predictable ways:

MemoryError

because merge builds the whole output in memory. Or it succeeds and the result has visible seams, because the tiles have different NoData values. Or the overlap regions are wrong, because the default picks the first tile rather than the best one. Or it silently produces garbage because one tile is in a different CRS and merge does not check.

Mosaicking is not a hard operation. It is an operation with four preconditions that nothing verifies for you.

Quick answer #

Check first, then merge:

import rasterio
from rasterio.merge import merge
from pathlib import Path

paths = sorted(Path("tiles").glob("*.tif"))
sources = [rasterio.open(p) for p in paths]

# the four preconditions
assert len({s.crs for s in sources}) == 1,      "tiles have different CRS"
assert len({s.dtypes[0] for s in sources}) == 1, "tiles have different dtypes"
assert len({s.count for s in sources}) == 1,     "tiles have different band counts"
assert len({s.res for s in sources}) == 1,       "tiles have different resolutions"

mosaic, transform = merge(sources, nodata=sources[0].nodata, method="first")

profile = sources[0].profile.copy()
profile.update(height=mosaic.shape[1], width=mosaic.shape[2],
               transform=transform, compress="deflate", tiled=True,
               blockxsize=512, blockysize=512)

with rasterio.open("mosaic.tif", "w", **profile) as dst:
    dst.write(mosaic)

for s in sources:
    s.close()

Step-by-step solution #

1. Audit the tiles before merging anything #

import rasterio, pandas as pd
from pathlib import Path

def audit_tiles(folder, pattern="*.tif"):
    rows = []
    for p in sorted(Path(folder).glob(pattern)):
        with rasterio.open(p) as src:
            rows.append({
                "file": p.name, "crs": str(src.crs), "dtype": src.dtypes[0],
                "bands": src.count, "res": src.res[0], "nodata": src.nodata,
                "w": src.width, "h": src.height,
            })
    df = pd.DataFrame(rows)
    for col in ["crs", "dtype", "bands", "res", "nodata"]:
        vals = df[col].unique()
        mark = "✓" if len(vals) == 1 else "✗"
        print(f"  {mark} {col:<8} {list(vals)[:4]}{' …' if len(vals) > 4 else ''}")
    return df

df = audit_tiles("tiles/")

  ✓ crs      ['EPSG:27700']
  ✓ dtype    ['int16']
  ✓ bands    [1]
  ✓ res      [25.0]
  ✗ nodata   [-9999.0, None, 0.0]

That last line is the one that produces seams. Three tiles have no NoData tag and two use 0, so merge will treat their fill values as real data and paint them over their neighbours.

2. Normalise NoData before merging #

merge uses nodata to decide which cells are "empty" and therefore eligible to be filled from another tile. Get it wrong and the result is a mosaic with rectangular blocks of the wrong value.

NODATA = -9999

def normalise_nodata(paths, out_dir, nodata=NODATA):
    """Give every tile the same NoData value, converting old fills."""
    out_dir = Path(out_dir); out_dir.mkdir(parents=True, exist_ok=True)
    fixed = []
    for p in paths:
        with rasterio.open(p) as src:
            arr = src.read()
            profile = src.profile.copy()
            if src.nodata is not None and src.nodata != nodata:
                arr[arr == src.nodata] = nodata
            elif src.nodata is None:
                pass          # nothing to convert — but see the note below
            profile.update(nodata=nodata)
        dst = out_dir / p.name
        with rasterio.open(dst, "w", **profile) as out:
            out.write(arr)
        fixed.append(dst)
    return fixed

The src.nodata is None branch is deliberately a no-op, because guessing is worse than declaring. A tile with no NoData tag may genuinely be full — a land tile with data everywhere — or it may be using 0 as an undeclared fill. Only the data's provenance tells you which, and assuming 0 means "empty" on an elevation raster erases sea level. When you do know, set it explicitly:

profile.update(nodata=nodata)
arr[arr == 0] = nodata            # only when you know 0 was the fill

3. Choose how overlaps are resolved #

mosaic, transform = merge(sources, method="first")

"First" is a defensible default only if the tiles are genuinely interchangeable. When one tile is newer, sort so it comes last and use "last":

sources = sorted(sources, key=lambda s: s.tags().get("ACQUISITION_DATE", ""))
mosaic, transform = merge(sources, method="last")     # newest wins in overlaps

A custom callable gets the destination and source arrays and can implement anything:

import numpy as np

def prefer_valid_then_mean(old_data, new_data, old_nodata, new_nodata, **kwargs):
    """Average where both tiles have data, otherwise take whichever is valid."""
    both = ~old_nodata & ~new_nodata
    only_new = old_nodata & ~new_nodata
    old_data[both] = (old_data[both] + new_data[both]) / 2
    old_data[only_new] = new_data[only_new]

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

Averaging the overlap removes the visible seam that "first" leaves wherever two tiles disagree slightly — common with tiles processed at different times.

4. Bound the output extent when you do not need everything #

merge defaults to the union of all inputs, which for a national tile set is a very large raster. Restrict it:

mosaic, transform = merge(
    sources,
    bounds=(320000, 380000, 400000, 440000),      # only this window
    res=25,
    nodata=-9999,
)

bounds is the single most effective defence against MemoryError here, because it limits the output array rather than the input reads.

5. Write with tiling and overviews #

A mosaic is a file people will pan around in. Write it so that is cheap:

profile.update(
    height=mosaic.shape[1], width=mosaic.shape[2], transform=transform,
    compress="deflate", predictor=2,        # predictor=2 helps integer data
    tiled=True, blockxsize=512, blockysize=512,
    BIGTIFF="IF_SAFER",                     # >4 GB output
)
with rasterio.open("mosaic.tif", "w", **profile) as dst:
    dst.write(mosaic)
    dst.build_overviews([2, 4, 8, 16, 32], rasterio.enums.Resampling.average)
    dst.update_tags(ns="rio_overview", resampling="average")

BIGTIFF="IF_SAFER" costs nothing and prevents the failure where a mosaic crosses 4 GB and the classic TIFF offsets overflow. Overviews use resampling — average for continuous data, mode or nearest for categorical.

Code examples #

Example 1: a merge that will not run out of memory #

rasterio.merge builds the output array in RAM. For a mosaic larger than memory, write it window by window instead:

import numpy as np, rasterio
from rasterio.merge import merge
from rasterio.windows import Window
from pathlib import Path

def merge_windowed(paths, dst_path, block=2048, nodata=-9999):
    srcs = [rasterio.open(p) for p in paths]
    try:
        assert len({s.crs for s in srcs}) == 1, "mixed CRS"
        res = srcs[0].res[0]
        xs = [b for s in srcs for b in (s.bounds.left, s.bounds.right)]
        ys = [b for s in srcs for b in (s.bounds.bottom, s.bounds.top)]
        bounds = (min(xs), min(ys), max(xs), max(ys))
        width  = int(round((bounds[2] - bounds[0]) / res))
        height = int(round((bounds[3] - bounds[1]) / res))
        transform = rasterio.transform.from_origin(bounds[0], bounds[3], res, res)

        profile = srcs[0].profile.copy()
        profile.update(width=width, height=height, transform=transform,
                       nodata=nodata, compress="deflate", tiled=True,
                       blockxsize=512, blockysize=512, BIGTIFF="IF_SAFER")

        with rasterio.open(dst_path, "w", **profile) as dst:
            for row in range(0, height, block):
                for col in range(0, width, block):
                    win = Window(col, row,
                                 min(block, width - col), min(block, height - row))
                    win_bounds = rasterio.windows.bounds(win, transform)
                    part, _ = merge(srcs, bounds=win_bounds, res=res, nodata=nodata)
                    dst.write(part, window=win)
                print(f"  rows {row:>7,} / {height:,}")
    finally:
        for s in srcs:
            s.close()
    return dst_path

merge_windowed(sorted(Path("tiles").glob("*.tif")), "mosaic.tif")

Peak memory is one block × block window rather than the whole mosaic, so a 200,000 × 200,000 output is no harder than a small one. merge(..., bounds=win_bounds) reads only the tiles that intersect the window, so the cost per block stays flat as the tile count grows.

The try/finally is not decoration: 240 open datasets are 240 file handles, and leaking them across a loop is how a batch job hits too many open files.

Example 2: a virtual mosaic — no output file at all #

Often you do not need a mosaic file; you need to read the tile set as if it were one raster. A VRT does that with a small XML file:

from osgeo import gdal
from pathlib import Path

paths = [str(p) for p in sorted(Path("tiles").glob("*.tif"))]
vrt = gdal.BuildVRT("mosaic.vrt", paths, options=gdal.BuildVRTOptions(
    resolution="highest", srcNodata=-9999, VRTNodata=-9999,
))
vrt = None          # flush to disk

import rasterio
with rasterio.open("mosaic.vrt") as src:
    print(src.width, src.height, src.count)    # behaves as one raster
    window = src.read(1, window=rasterio.windows.from_bounds(
        320000, 380000, 340000, 400000, transform=src.transform))

The VRT is a few kilobytes and takes seconds to build against hundreds of gigabytes of tiles. Everything downstream — clipping, zonal statistics, sampling — works against it unchanged.

Use a VRT when the tiles are stable and local. Materialise a real mosaic when you need a single portable file, when the tiles live on slow or remote storage, or when the same region will be read many times and the per-tile open overhead starts to dominate.

Example 3: merging tiles that do not share a CRS #

The precondition check fails; the fix is to reproject the odd ones out first:

from collections import Counter
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pathlib import Path

def unify_crs(paths, work_dir, resampling=Resampling.nearest):
    """Reproject any tile that disagrees with the majority CRS."""
    work_dir = Path(work_dir); work_dir.mkdir(parents=True, exist_ok=True)
    crs_of = {}
    for p in paths:
        with rasterio.open(p) as src:
            crs_of[p] = src.crs
    if None in crs_of.values():
        missing = [p.name for p, c in crs_of.items() if c is None]
        raise ValueError(f"tiles with no CRS: {missing}")

    target = Counter(crs_of.values()).most_common(1)[0][0]
    print(f"target CRS: {target}")

    out = []
    for p in paths:
        if crs_of[p] == target:
            out.append(p)
            continue
        dst = work_dir / p.name
        print(f"  reprojecting {p.name} from {crs_of[p]}")
        with rasterio.open(p) as src:
            transform, width, height = calculate_default_transform(
                src.crs, target, src.width, src.height, *src.bounds)
            profile = src.profile.copy()
            profile.update(crs=target, transform=transform, width=width, height=height)
            with rasterio.open(dst, "w", **profile) as d:
                for i in range(1, src.count + 1):
                    reproject(source=rasterio.band(src, i),
                              destination=rasterio.band(d, i),
                              src_transform=src.transform, src_crs=src.crs,
                              src_nodata=src.nodata,
                              dst_transform=transform, dst_crs=target,
                              dst_nodata=src.nodata, resampling=resampling)
        out.append(dst)
    return out, target

Two decisions here are deliberate. The target is the majority CRS, so the common case reprojects the two odd tiles rather than the 238 good ones. And a tile with no CRS raises rather than being assumed to match — an unlabelled tile is not evidence of anything, and quietly treating it as the majority CRS is how a mosaic ends up with one block of data from the wrong place. The full reprojection details are in how to reproject a raster.

Explanation #

Merging rasters looks like concatenation and is closer to a spatial join. The reason is that the output grid does not exist yet, and every input cell has to be placed on it.

When tiles share a CRS and resolution and their origins differ by a whole number of cells, the placement is exact: each input cell maps to exactly one output cell and no values change. This is the case for a properly produced national tile set, and it is why merging such a set is fast and lossless.

When origins are offset by a fraction of a cell — which happens as soon as tiles come from different processing runs, or after independent reprojections — the placement is not exact, and the merge resamples. Values change, and the change is invisible unless you check. This is the same alignment problem described in how to reproject a raster: a shared CRS does not imply a shared grid.

Overlaps are the second structural issue. Most tile sets overlap deliberately, because a buffer lets each tile be processed independently without edge artefacts — the same read-wide, write-narrow pattern used when splitting a large layer into tiles. In the overlap, two tiles both have valid data, and they usually disagree slightly. method="first" resolves this by tile order, which is arbitrary; the visible result is a seam wherever the disagreement exceeds the display's colour resolution. Averaging removes the seam, "max" or "min" implement a compositing rule, and a callable implements whatever the data actually requires.

NoData is what makes merging possible at all. A tile is a rectangle, its data usually is not, and the cells that are not data must be recognisable as such or they will overwrite a neighbour's real values. This is why the NoData audit matters more here than anywhere else in raster work: a mosaic is precisely the operation where one tile's fill value lands on top of another tile's measurements.

Finally, there is the question of whether to mosaic at all. A mosaic is a large file that duplicates data you already have, and it must be rebuilt whenever a tile is updated. A VRT is a pointer that costs kilobytes and updates instantly. Modern formats push this further: a Cloud-Optimized GeoTIFF mosaic with internal overviews serves windowed reads efficiently over HTTP, and STAC catalogues describe tile sets without materialising them at all. Build the mosaic when you need one file to hand to someone or one file to read many times; otherwise, read the tiles through a VRT and skip the copy.

Edge cases or notes #

• merge holds the whole output in memory. Use bounds= to limit it, or write window by window as in Example 1.
• Mixed NoData produces rectangular seams. Normalise before merging; the audit in step 1 catches it.
• A tile with nodata=None is treated as fully valid, so its fill blocks paint over neighbours.
• method="first" depends on argument order, which depends on your sort. Sort deliberately.
• Output over 4 GB needs BIGTIFF="IF_SAFER". Without it, the write fails or produces an unreadable file.
• Hundreds of open datasets exhaust file handles. Close them in a finally, or use a context-manager stack.
• res= forces a resolution and resamples tiles that disagree — check the method before relying on it.
• Overviews are not inherited from the tiles. Build them on the mosaic.
• gdalbuildvrt + gdal_translate -co TILED=YES is the standard command-line route and streams throughout.
• rio merge and gdal_merge.py exist but are less flexible than the Python API for anything beyond "first wins".

Internal links #

• The raster data model explained — NoData, dtype and the transform that must agree
• How to reproject a raster in Python with Rasterio — unifying CRS before a merge
• Raster resampling explained — what happens when grids do not align
• How to split a large layer into tiles — the operation this reverses
• How to batch process rasters with Rasterio — processing tiles before merging them
• How to merge many shapefiles into one file — the vector equivalent
• OSError: too many open files in a batch GIS job — what 240 open datasets cause
• How to reduce GIS file size in Python — compression choices for a large mosaic

FAQ #

Why does my mosaic have visible seams? #

The tiles have different NoData values, so one tile's fill is being treated as data and painted over its neighbour. Audit nodata across all tiles and normalise before merging.

Can I merge tiles in different coordinate systems? #

Not directly — merge does not check or reproject. Reproject the minority tiles to the majority CRS first, as in Example 3.

How do I merge rasters too large for memory? #

Write the output window by window, calling merge(sources, bounds=window_bounds) per block, or build a VRT and skip materialising the mosaic entirely.

What does method="first" mean? #

In an overlap, the first source in the list that has valid data for a cell wins. Since that depends on your argument order, sort deliberately or use an explicit rule.

Should I build a mosaic or a VRT? #

A VRT when the tiles are stable and local — it costs kilobytes and builds in seconds. A real mosaic when you need one portable file, or when per-tile open overhead dominates repeated reads.

Why is my merged output enormous? #

merge defaults to the union of all inputs. Pass bounds= to restrict it, and always write with compress="deflate" and tiled=True.

Do I need to rebuild overviews after merging? #

Yes. Overviews belong to a file, not to the data, so a new mosaic has none until you call build_overviews.