How to Mosaic Rasters: Join Tiles into One Raster in QGIS and Python

This guide has 2 routes. Pick one, or read both. Each route is complete on its own.

Satellite scenes, elevation models and land cover maps come in tiles. A district that crosses a tile edge needs two or four of them, and you want one raster to work with. Joining tiles into one is called mosaicking.

This guide builds a mosaic of the Copernicus DEM for the area around Ludhiana district in Punjab. The district sits where four 1 degree tiles meet, so it is a good test. The next guides use the result to clip a raster to a district and to reproject it. If you want to skip the download, the file is in the data zip linked in the guide on clipping.

What the tiles must have in common

A mosaic only works when the tiles agree on these things:

The CRS. Tiles in different coordinate systems must be reprojected first.
The cell size. Tiles with different cell sizes are resampled, and that can blur one of them.
The data type and the NoData value. A mix of integer and decimal tiles gives a decimal mosaic.
The grid. Tiles that sit on the same grid line up exactly. Tiles with a half-cell shift show a seam.
The number of bands.

Check this before you start. Tiles from one product, such as one DEM or one Sentinel-2 tile grid, always agree. Tiles from different products often do not.

Two ways to join: copy or point

A real mosaic writes a new file that holds all the values. It is self-contained, you can move it, and it needs disk space equal to the sum of the tiles.

A virtual raster, or VRT, writes only a small text file that lists the tiles and where each sits. It needs no disk space and it is instant, but it needs the original tiles to stay where they are. It is a good choice for looking at the data, and for feeding other tools. For a final product that you send to someone, make a real mosaic.

Route A: QGIS

Add the tiles to the project first, with Layer > Add Layer > Add Raster Layer.

For a real mosaic, open Raster > Miscellaneous > Merge. Choose the tiles under Input layers, and save the result under Merged. Open the advanced parameters if you want to set the NoData value to treat as empty in the input, and the value to assign to empty cells in the output. Leave Place each input file into a separate band unticked, otherwise you get a layer with one band per tile and no mosaic. Click Run.

For a virtual mosaic, open Raster > Miscellaneous > Build Virtual Raster. Choose the tiles, keep Place each input file into a separate band unticked, and save the result under Virtual with the extension .vrt.

Cut this DEM into a left and a right half and join the halves again with either tool, and the values come back identical to the original. Where tiles overlap, these GDAL tools take the value of the last tile in the list, so if the order matters, order the layers in the dialog. Not every tool agrees: the Python function below keeps the first value by default.

Route B: Python

The script reads only the part of each tile that your box touches, so it works even when the tiles are tens of megabytes each. You need Python with requests, numpy and rasterio.

import requests
import numpy as np
import rasterio
from rasterio.merge import merge

# Step 1: the area you want, as west, south, east, north in longitude and latitude.
# This box surrounds Ludhiana district, Punjab.
box = (75.35, 30.55, 76.34, 31.03)

# Step 2: ask the Earth Search catalogue which Copernicus DEM tiles touch the box
search = {"collections": ["cop-dem-glo-30"], "bbox": list(box), "limit": 20}
url = "https://earth-search.aws.element84.com/v1/search"
tiles = requests.post(url, json=search, timeout=60).json()["features"]
names = [t["id"].replace("Copernicus_DSM_COG_10_", "").replace("_DEM", "") for t in tiles]
print(len(tiles), "tiles:", sorted(names))

4 tiles: ['N30_00_E075_00', 'N30_00_E076_00', 'N31_00_E075_00', 'N31_00_E076_00']

Four tiles touch the box. The names give their south-west corners: N30_00_E075_00 is the tile from 30 N and 75 E, and so on.

# Step 1: the catalogue gives an s3:// address. A normal web address for the same file is what rasterio opens.
def web_address(s3_address):
    bucket, key = s3_address.replace("s3://", "").split("/", 1)
    return f"https://{bucket}.s3.amazonaws.com/{key}"

# Step 2: open the four tiles online. Nothing is downloaded yet.
files = [rasterio.open(web_address(t["assets"]["data"]["href"])) for t in tiles]

# Step 3: check that they agree. Every line should show one value only.
print("CRS:", {f.crs.to_string() for f in files})
print("cell size:", {round(f.res[0], 7) for f in files})
print("data type:", {f.dtypes[0] for f in files})
print("bands:", {f.count for f in files})

CRS: {'EPSG:4326'}
cell size: {0.0002778}
data type: {'float32'}
bands: {1}

One value on every line means the tiles match. If a line shows two values, fix that first.

# Step 1: join the tiles for your box as it is. Cells that no tile covers are filled with -9999.
plain, plain_transform = merge(files, bounds=box, nodata=-9999)
print("plain box: shape (bands, rows, columns):", plain.shape, "| cells with no data:", int((plain == -9999).sum()))

plain box: shape (bands, rows, columns): (1, 1728, 3564) | cells with no data: 5291

The four tiles cover the whole box, so there should be no empty cells at all. There are thousands: an empty last row and an empty last column. The cause is the grid of the tiles. The Copernicus DEM tiles start half a cell away from the whole degrees, at 74.99986 and not at 75, so the cell grid does not line up with the numbers in your box. When the box is not on the grid, rasterio.merge adds a row and a column that hold nothing, and it also resamples every cell by a fraction of a cell. The cure is to move the box onto the grid of the tiles, outwards by at most one cell.

import math

# Step 1: the grid of the tiles. Its origin is the top-left corner of the whole set of tiles, and the cell size is the same for all.
cell = files[0].res[0]
origin_x = min(f.bounds.left for f in files)
origin_y = max(f.bounds.top for f in files)

# Step 2: move each side of the box outwards to the nearest grid line
west = origin_x + math.floor((box[0] - origin_x) / cell) * cell
east = origin_x + math.ceil((box[2] - origin_x) / cell) * cell
north = origin_y - math.floor((origin_y - box[3]) / cell) * cell
south = origin_y - math.ceil((origin_y - box[1]) / cell) * cell
snapped = (west, south, east, north)
print("box on the grid of the tiles:", [round(v, 5) for v in snapped])

# Step 3: merge again with the snapped box
mosaic, transform = merge(files, bounds=snapped, nodata=-9999)
for f in files:
    f.close()
print("snapped box: shape (bands, rows, columns):", mosaic.shape, "| cells with no data:", int((mosaic == -9999).sum()))
print("height: lowest", round(float(mosaic.min()), 1), "m | highest", round(float(mosaic.max()), 1), "m")

box on the grid of the tiles: [75.34986, 30.54986, 76.34014, 31.03014]
snapped box: shape (bands, rows, columns): (1, 1729, 3565) | cells with no data: 0
height: lowest 214.3 m | highest 298.3 m

With the box on the grid there are no empty cells, and no cell was resampled: every value in the mosaic is a value of a tile. Always count the empty cells after a mosaic. An empty strip at the edge is easy to miss in a picture, and it spoils statistics and any later clip. The same rule holds for any tiles whose grid does not start on a round number, which is common for elevation data.

Save it

# Step 1: write the mosaic to a GeoTIFF. The CRS and the transform place it in the right position.
profile = {
    "driver": "GTiff", "height": mosaic.shape[1], "width": mosaic.shape[2], "count": 1,
    "dtype": "float32", "crs": "EPSG:4326", "transform": transform, "nodata": -9999,
    "compress": "deflate", "tiled": True,
}
with rasterio.open("dem_ludhiana.tif", "w", **profile) as dst:
    dst.write(mosaic)

# Step 2: read it back to check
with rasterio.open("dem_ludhiana.tif") as check:
    print("saved:", check.width, "x", check.height, "cells | CRS EPSG", check.crs.to_epsg(), "| left", round(check.bounds.left, 4), "| top", round(check.bounds.top, 4))

saved: 3565 x 1729 cells | CRS EPSG 4326 | left 75.3499 | top 31.0301

Open dem_ludhiana.tif in QGIS next to a basemap and look for seams along the lines at 76 E and 31 N, where the four tiles meet. There should be none.

Where tiles overlap

Sentinel-2 tiles overlap by about 10 km, and scenes from different dates overlap completely. When several rasters cover the same cell, the mosaic has to choose a value. rasterio.merge takes a method argument for this:

"first", the default, keeps the value of the first raster in the list that has data.
"last" keeps the value of the last one.
"min" and "max" keep the lowest or the highest.
"sum" adds the values, and "count" counts how many rasters have data.

For scenes of different dates, first and last are rarely what you want. A median over time, as in the NDVI time series, is the usual answer for cloud-free composites. For tiles of one product, which meet at an edge and do not overlap, the choice does not matter.

Common mistakes

Joining tiles in different CRSs, which gives a mosaic with strips in the wrong place.
Ticking Place each input file into a separate band in QGIS and getting a layer with one band per tile.
Moving the tiles after building a virtual raster, which leaves an empty layer.
Giving merge a box that is not on the grid of the tiles, which adds an empty row and column and resamples the cells, and not counting the empty cells.
Mixing products with different cell sizes and not noticing the resampling.
Forgetting the NoData value, so the filled -9999 values become real heights in the statistics.