Merge Raster Tiles into a Mosaic in PyQGIS

Raster data rarely arrives as one file. National DEMs come as hundreds of 1 km or 10 km tiles, orthophotos as map sheets, satellite imagery as scenes cut along orbit paths. Analyses — slope, viewshed, zonal statistics across a catchment that spans six tiles — need a single continuous raster. Merging is simple in principle and full of small traps in practice: tiles in different CRSs, slightly different resolutions, inconsistent NoData values, overlaps where two tiles disagree.

This recipe belongs to Raster Analysis Workflows. It inventories a folder of tiles, checks they are compatible, merges them with gdal:merge into a compressed tiled GeoTIFF, handles NoData and overlaps deliberately, and compares the result with a virtual raster that merges nothing at all.

From tiles to one rasterA folder of tiles is first inventoried: CRS, pixel size, data type, band count and NoData for each. Incompatible tiles are reprojected or resampled. The compatible set is merged with gdal:merge into one GeoTIFF, with NoData declared so gaps between tiles stay transparent, and compression and tiling applied. Alternatively, a VRT references the tiles without copying them.Inventory, align, mergetilesdem_*.tifhundredsinventoryCRS · pixel sizedtype · NoDatagdal:merge → GeoTIFFone file, compressed, tiledor: VRTreferences tiles, no copy

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series, with the GDAL Processing provider (enabled by default).
  • A folder of raster tiles. The example uses 1 km DEM tiles named dem_*.tif.
  • Enough disk space for the output: a merged GeoTIFF is roughly the sum of the inputs, less with compression.

Inventory the tiles

Merging assumes every tile shares a CRS, pixel size, data type and band count. Checking that first takes seconds and prevents a merged raster that is subtly wrong.

from pathlib import Path
from collections import Counter
from qgis.core import QgsRasterLayer

tiles = sorted(Path("/data/dem/tiles").glob("dem_*.tif"))
info = []
for p in tiles:
    lyr = QgsRasterLayer(str(p), p.stem)
    if not lyr.isValid():
        print("cannot open", p.name)
        continue
    prov = lyr.dataProvider()
    info.append({
        "path": str(p), "crs": lyr.crs().authid(),
        "res": (round(lyr.rasterUnitsPerPixelX(), 6), round(lyr.rasterUnitsPerPixelY(), 6)),
        "dtype": prov.dataType(1), "bands": lyr.bandCount(),
        "nodata": prov.sourceNoDataValue(1) if prov.sourceHasNoDataValue(1) else None,
    })

for key in ("crs", "res", "dtype", "bands", "nodata"):
    print(key, Counter(str(i[key]) for i in info).most_common())

Breakdown: Printing a count of distinct values for each property shows at a glance whether the set is uniform. One line per property with a single value means the tiles are compatible; any property with two or more values names the problem. Mixed CRSs need reprojection first, as in batch reprojecting raster datasets. Mixed resolutions need resampling to a common pixel size, covered in resampling and aligning rasters. Mixed NoData values can be unified during the merge, as below.

Merge into one GeoTIFF

gdal:merge reads the tiles and writes one raster covering their combined extent. Declaring NoData and choosing creation options decide how usable the result is.

import processing

paths = [i["path"] for i in info]
out = processing.run("gdal:merge", {
    "INPUT": paths,
    "PCT": False,
    "SEPARATE": False,
    "NODATA_INPUT": -9999,
    "NODATA_OUTPUT": -9999,
    "DATA_TYPE": 5,                     # Float32
    "OPTIONS": "COMPRESS=DEFLATE|PREDICTOR=3|TILED=YES|BIGTIFF=IF_SAFER",
    "EXTRA": "",
    "OUTPUT": "/data/dem/dem_mosaic.tif",
})["OUTPUT"]

mosaic = QgsRasterLayer(out, "DEM mosaic")
print(mosaic.width(), "x", mosaic.height(), "pixels;", mosaic.extent().toString(0))

Breakdown: NODATA_INPUT tells the merge which input value means "no data", so those pixels do not overwrite valid data from another tile where tiles overlap. NODATA_OUTPUT declares the value in the result, so gaps between tiles are transparent rather than black. The data type enum follows GDAL's order — 5 is Float32 in the Processing dialog's list; check processing.algorithmHelp("gdal:merge") on your version. Deflate compression with predictor 3 shrinks floating-point elevation data substantially, tiling makes the file fast to display and clip, and BIGTIFF=IF_SAFER avoids failure when the output exceeds 4 GB.

Decide what happens where tiles overlap

Tiles often overlap by a few pixels, and the values there may differ slightly — different acquisition dates, different processing. gdal:merge lets the last tile win, so the order of the input list decides the result.

Last tile wins in overlapsTwo tiles overlap in a strip. With gdal:merge the tile listed later overwrites the earlier one in the overlap, unless the later tile's pixel is NoData. Ordering the list, for example by acquisition date so the newest is last, makes the outcome deliberate. A VRT follows the same rule with the last source on top.Order the list to choose the winnertile A (2019)tile B (2023)overlap: B winssort inputsoldest firstnewest last

import re

def acquisition_year(path):
    m = re.search(r"_(\d{4})\.tif$", path)
    return int(m.group(1)) if m else 0

ordered = sorted(paths, key=acquisition_year)      # newest last → newest wins
processing.run("gdal:merge", {
    "INPUT": ordered, "NODATA_INPUT": -9999, "NODATA_OUTPUT": -9999,
    "DATA_TYPE": 5, "OPTIONS": "COMPRESS=DEFLATE|PREDICTOR=3|TILED=YES",
    "OUTPUT": "/data/dem/dem_mosaic_newest.tif"})

Breakdown: Sorting by acquisition year, extracted here from the file name, puts the newest tiles last so their values win in overlaps. Any rule can be encoded the same way — highest accuracy class, lowest cloud cover. Pixels that are NoData in the later tile do not overwrite valid values from earlier ones, so a ragged-edged newer scene fills only where it has data. For blending rather than overwriting — averaging overlaps to hide seams — compute the overlap separately with the raster calculator, which is rarely worth the effort for DEMs but sometimes matters for imagery.

Check the result

A merged raster should have no unexpected holes, no visible seams and statistics consistent with the tiles.

from qgis.core import QgsRasterBandStats

stats = mosaic.dataProvider().bandStatistics(1, QgsRasterBandStats.All)
print(f"min {stats.minimumValue:.1f}  max {stats.maximumValue:.1f}  mean {stats.mean:.1f}")

block = mosaic.dataProvider().block(1, mosaic.extent(), 400, 400)
nodata_share = sum(block.isNoData(r, c) for r in range(400) for c in range(400)) / 160000
print(f"NoData in a 400×400 overview sample: {nodata_share:.1%}")

Breakdown: Minimum and maximum outside the tiles' range point to a NoData value that was not declared — a −9999 counted as an elevation. Sampling a low-resolution overview of the mosaic for NoData shows whether gaps are where expected (sea, outside the survey area) or appear as lines between tiles, which indicates misaligned tiles or inconsistent NoData. A hillshade of the mosaic, as in generating slope, aspect and hillshade, reveals seams no statistic shows.

Keep the mosaic reproducible

A mosaic is a derived product, and questions about it come later: which tiles went in, which won in overlaps, what NoData means. Recording that next to the file costs one small step.

import json
from datetime import datetime, timezone

record = {
    "created": datetime.now(timezone.utc).isoformat(timespec="seconds"),
    "inputs": ordered, "nodata": -9999, "data_type": "Float32",
    "overlap_rule": "later input wins; inputs sorted by acquisition year",
}
with open("/data/dem/dem_mosaic_newest.json", "w") as fh:
    json.dump(record, fh, indent=2)

Breakdown: A JSON sidecar listing inputs in merge order, the NoData value and the overlap rule lets anyone reproduce the mosaic or explain an odd value at a seam. It is also the natural place to note tile versions when a supplier reissues tiles. For mosaics rebuilt on a schedule, write the record in the same script that runs the merge, so the two can never disagree.

When a VRT is better than merging

Merging copies every pixel into a new file. A virtual raster (VRT) is a small XML file that references the tiles and presents them as one raster, with no copying.

Merge or VRTA merged GeoTIFF is a standalone copy: portable, fast to read, but duplicates storage and goes stale when tiles are updated. A VRT is a few kilobytes referencing the original tiles: no duplication, always current, but depends on the tiles staying in place and can be slower over many small files. Many workflows build a VRT for analysis and merge only for delivery.Copy everything, or reference itmerged GeoTIFFstandalone, portablefast readsdoubles storage, can go staleVRTkilobytes, no copyalways currenttiles must stay put

For analysis on a machine that holds the tiles, a VRT is usually the better choice: it is instant to build, needs no extra space and reflects tile updates automatically. For delivering a mosaic to someone else, or for very many small tiles where opening hundreds of files is slow, merge. Building a virtual raster covers the VRT route, and a common pattern combines both: build a VRT, check it, then translate it to a GeoTIFF in one step.

Merge only what an area needs

For a study area covering a few tiles out of hundreds, select the intersecting tiles first rather than merging the whole set and clipping afterwards.

from qgis.core import QgsProject, QgsRectangle

area = QgsProject.instance().mapLayersByName("catchment")[0]
target = area.extent()
target.grow(100)                                     # small margin in map units

needed = []
for i in info:
    ext = QgsRasterLayer(i["path"], "x").extent()
    if ext.intersects(target):
        needed.append(i["path"])
print(len(needed), "of", len(info), "tiles intersect the catchment")

Breakdown: Comparing each tile's extent with the study area's bounding box, plus a margin so edge pixels are not lost, selects the few tiles that matter. Merging those and then clipping by the catchment polygon is far faster than processing the national set. For repeated selections, build a tile index layer once with gdal:tileindex and select by location.

QGIS version compatibility

gdal:merge is available with the same parameters on QGIS 3.34 LTR, 3.40 LTR and QGIS 4; the data type enum and the creation OPTIONS string follow the bundled GDAL. Very long input lists are written to a temporary file list automatically on recent releases, avoiding command-line length limits on Windows.

Troubleshooting

  • Black stripes between tiles. NoData was not declared; set NODATA_INPUT and NODATA_OUTPUT.
  • The mosaic is shifted or blurry in parts. Tiles had different CRSs or resolutions; fix them before merging.
  • The merge fails on Windows with many files. Update QGIS, or merge in batches and merge the results.
  • The output is huge. Compression options were not set; use COMPRESS=DEFLATE and tiling.

Conclusion

Inventory tiles for CRS, resolution, data type and NoData before merging, order inputs so the right tile wins in overlaps, merge with declared NoData and compression into a tiled GeoTIFF, check statistics and seams, select only the tiles an area needs, and prefer a VRT for local analysis.

Frequently Asked Questions

Can I merge rasters with different numbers of bands? Not into one band stack with gdal:merge; make band counts consistent first, or use SEPARATE to put each input in its own band.

Does merging change pixel values? No, unless data types differ and values are converted, or overlaps are overwritten.

How do I merge imagery with colour tables? Set PCT to keep the colour table from the first input; all inputs should share it.

Can I merge rasters in different CRSs directly?gdalwarp can, by reprojecting during the merge; gdal:warpreproject with multiple inputs is the Processing route.