How to Merge and Mosaic Rasters in Python with Rasterio

Problem statement

Large raster data arrives as tiles. The Copernicus 30 m elevation model, for instance, is one GeoTIFF per degree of latitude and longitude, named like Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif, and every analysis you want to run needs one continuous surface. Here are the four tiles that cover the Pennines between 53° and 55°N (the setup is under The data used on this page).

The one-liner exists:

from pathlib import Path

import rasterio
from rasterio.merge import merge

paths = sorted(Path("tiles").glob("*.tif"))
mosaic, transform = merge([rasterio.open(p) for p in paths])
print(f"{len(paths)} tiles → {mosaic.shape[2]:,} × {mosaic.shape[1]:,} cells, {mosaic.nbytes / 1e6:,.0f} MB in memory")

4 tiles → 4,800 × 7,200 cells, 138 MB in memory

and it fails in four predictable ways. At national scale it is a MemoryError, because merge builds the whole output in memory — four tiles already cost the figure above. 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, 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()
print(f"mosaic.tif: {profile['width']:,} × {profile['height']:,}, {Path('mosaic.tif').stat().st_size / 1e6:,.0f} MB on disk")

mosaic.tif: 4,800 × 7,200, 115 MB on disk

Audit tiles, normalise NoData, choose the overlap rule

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 = [v.item() if hasattr(v, "item") else v for v in 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:4326']
  ✓ dtype    ['float32']
  ✓ bands    [1]
  ✓ res      [0.00041666666666666664]
  ✓ nodata   [None]

Every check passes, and the one that deserves a second look is the NoData line: none of these tiles declares a NoData value at all. On a tile set where some tiles used 0 and others -9999, that line is what produces seams — merge treats an undeclared fill value as real data and paints it over its neighbours. The res column hides something too: it shows only the x resolution, and Copernicus pixels are not square — step 1 of example 1 finds out what that costs.

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. These tiles show the trap exactly:

import numpy as np

with rasterio.open("mosaic.tif") as src:
    z = src.read(1)

print(f"exactly 0: {(z == 0).mean():.1%} of cells — the sea")
print(f"below 0:   {(z < 0).sum():,} cells, down to {z.min():.1f} m")

exactly 0: 4.0% of cells — the sea
below 0:   8,843 cells, down to -12.2 m

Declaring nodata=0 here would turn the sea into holes, and a flood model or a coastline would silently go with it — and the cells below zero show that 0 is not even the floor of the real values. When you do know a value is a fill, 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. (The Copernicus tiles here have no overlap to resolve: each tile's bounds start half a cell past the whole degree, so neighbours meet edge to edge.) 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:

sources = [rasterio.open(p) for p in paths]
part, part_transform = merge(sources, bounds=(-2.5, 53.5, -1.5, 54.5))      # only this window, in the tiles' CRS
print(f"bounded merge: {part.shape[2]:,} × {part.shape[1]:,} cells, {part.nbytes / 1e6:,.0f} MB")

bounded merge: 2,400 × 3,600 cells, 35 MB

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=3,        # 3 for floating point; 2 for integers
    tiled=True, blockxsize=512, blockysize=512,
    BIGTIFF="IF_SAFER",                     # >4 GB output
)
with rasterio.open("mosaic_ovr.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")
for name in ("mosaic.tif", "mosaic_ovr.tif"):
    print(f"{name:15} {Path(name).stat().st_size / 1e6:6,.0f} MB")

mosaic.tif         115 MB
mosaic_ovr.tif     121 MB

The floating-point predictor and the overviews pull in opposite directions on size; the overviews are what make a zoomed-out view fast. 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.

A merge that will not run out of memory

The data used on this page

Four tiles of the Copernicus 30 m digital elevation model (GLO-30), downloaded from its public bucket on AWS into tiles/: 53–55°N, 3–1°W, the Pennines from the Peak District to Teesdale, with sea in two corners: the Irish Sea in the west and, in the north-east, the North Sea off County Durham, which holds most of the zero-valued cells. For example 3, delivery/SD_bng.tif is one of those tiles reprojected to British National Grid — the kind of tile a supplier working in the national grid would send.

from pathlib import Path

import numpy as np
import rasterio
import requests
from rasterio.warp import Resampling, calculate_default_transform, reproject

TILE = "https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{0}_DEM/Copernicus_DSM_COG_10_{0}_DEM.tif"
Path("tiles").mkdir(exist_ok=True)
for t in ("N53_00_W003_00", "N53_00_W002_00", "N54_00_W003_00", "N54_00_W002_00"):
    dst = Path("tiles") / f"Copernicus_DSM_COG_10_{t}_DEM.tif"
    if not dst.exists():
        dst.write_bytes(requests.get(TILE.format(t), timeout=600).content)

Path("delivery").mkdir(exist_ok=True)
with rasterio.open("tiles/Copernicus_DSM_COG_10_N54_00_W003_00_DEM.tif") as src:
    transform, width, height = calculate_default_transform(src.crs, "EPSG:27700", src.width, src.height,
                                                           *src.bounds, resolution=30)
    profile = src.profile | {"crs": "EPSG:27700", "transform": transform, "width": width, "height": height}
    with rasterio.open("delivery/SD_bng.tif", "w", **profile) as dst:
        reproject(rasterio.band(src, 1), rasterio.band(dst, 1), resampling=Resampling.bilinear)
print(sorted(p.name for p in Path("tiles").glob("*.tif")))

['Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif', 'Copernicus_DSM_COG_10_N53_00_W003_00_DEM.tif', 'Copernicus_DSM_COG_10_N54_00_W002_00_DEM.tif', 'Copernicus_DSM_COG_10_N54_00_W003_00_DEM.tif']

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(paths, "mosaic_windowed.tif")
with rasterio.open("mosaic.tif") as a, rasterio.open("mosaic_windowed.tif") as b:
    print(f"in-memory merge: {a.width:,} × {a.height:,}, bounds {tuple(round(v, 3) for v in a.bounds)}")
    print(f"windowed merge:  {b.width:,} × {b.height:,}, bounds {tuple(round(v, 3) for v in b.bounds)}")

  rows       0 / 4,800
  rows   2,048 / 4,800
  rows   4,096 / 4,800
in-memory merge: 4,800 × 7,200, bounds (-3.0, 53.0, -1.0, 55.0)
windowed merge:  4,800 × 4,800, bounds (-3.0, 53.0, -1.0, 55.0)

The two should be the same raster, and they are not. res = srcs[0].res[0] takes the x resolution and uses it for both axes — fine for square pixels, wrong for these: Copernicus tiles north of 50° are 1.5 arc-seconds wide and 1 arc-second tall. The windowed version therefore has a third fewer rows, each squeezed from a taller strip of ground. Use both resolutions:

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

        profile = srcs[0].profile.copy()
        profile.update(width=width, height=height, transform=transform,
                       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))
                    part, _ = merge(srcs, bounds=rasterio.windows.bounds(win, transform), res=(xres, yres))
                    dst.write(part, window=win)
    finally:
        for s in srcs:
            s.close()
    return dst_path


merge_windowed(paths, "mosaic_windowed.tif")
with rasterio.open("mosaic.tif") as a, rasterio.open("mosaic_windowed.tif") as b:
    print(f"windowed merge now {b.width:,} × {b.height:,}; identical to the in-memory merge: "
          f"{a.transform.almost_equals(b.transform) and np.array_equal(a.read(1), b.read(1))}")

windowed merge now 4,800 × 7,200; identical to the in-memory merge: True

Peak memory is one block × block window rather than the whole mosaic, so an output far larger than memory 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: hundreds of open datasets are hundreds of 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:

import subprocess

subprocess.run(["gdalbuildvrt", "-q", "mosaic.vrt", *map(str, paths)], check=True)   # or osgeo.gdal.BuildVRT
with rasterio.open("mosaic.vrt") as src:
    print(f"mosaic.vrt: {Path('mosaic.vrt').stat().st_size / 1e3:.1f} kB, "
          f"{src.width:,} × {src.height:,}, {src.count} band")        # behaves as one raster
    window = src.read(1, window=rasterio.windows.from_bounds(-2.1, 53.9, -1.9, 54.1, transform=src.transform))
    print(f"a window across all four tiles: {window.shape}, {window.min():.0f}–{window.max():.0f} m")

mosaic.vrt: 2.5 kB, 4,800 × 7,200, 1 band
a window across all four tiles: (720, 480), 84–504 m

The VRT is a few kilobytes, whatever the size of the tiles behind it. 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


tiles_mixed = paths + [Path("delivery/SD_bng.tif")]
unified, target = unify_crs(tiles_mixed, "work")
print([Path(p).name for p in unified])

target CRS: EPSG:4326
  reprojecting SD_bng.tif from EPSG:27700
['Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif', 'Copernicus_DSM_COG_10_N53_00_W003_00_DEM.tif', 'Copernicus_DSM_COG_10_N54_00_W002_00_DEM.tif', 'Copernicus_DSM_COG_10_N54_00_W003_00_DEM.tif', 'SD_bng.tif']

The odd tile is the setup's copy of one Copernicus tile in British National Grid, standing in for a supplier who delivers in the national grid. Two decisions here are deliberate. The target is the majority CRS, so the common case reprojects the odd tiles rather than the 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.

Merging is closer to a spatial join

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.

Memory, mixed NoData and seams

• 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 hundreds of 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.