Mosaicking CHM Tiles Without Edge Seams in Python

A canopy height mosaic is where every earlier compromise becomes visible at once. Tiles that were classified independently show a height step at their shared edge; tiles whose grids are a fraction of a pixel apart force a resample that smears the whole surface; tiles with inconsistent NoData punch black rectangles through the forest. This guide covers assembling tiles into one regional raster whose seams cannot be found, as the final assembly step of Canopy Height Model Creation within Canopy Height Modeling & Terrain Extraction.

When to use this approach

Build a VRT first, always. It is free, it fails loudly when the grids disagree, and it can be turned into a GeoTIFF later with one command. Reach for feathering only as a last resort — on a height surface a feathered seam does not fix a bad edge, it just makes it harder to find.

The three conditions for a seamless mosaic

A mosaic is seamless when three things hold, and each has its own failure signature.

1. Identical grid. Same CRS, same cell size, and origins that differ by a whole number of cells. Otherwise every read resamples.
2. Consistent NoData. One value, declared in the profile of every tile, and never a real elevation like 0.
3. No conflicting overlap. Either the tiles abut exactly, or where they overlap the values agree — which they will only if the halo discipline from Tiling and Buffering Large LiDAR Surveys with PDAL was followed.

Minimal reproducible example

Validate first, then build. The validation is what turns a silent resample into an error message.

import glob

import rasterio


def assert_mosaicable(tifs: list[str], tol: float = 1e-6) -> dict:
    """Fail loudly if tiles do not share a CRS, a cell size, and a common grid."""
    ref_res = ref_crs = ref_nodata = None
    for path in sorted(tifs):
        with rasterio.open(path) as src:
            res = (abs(src.transform.a), abs(src.transform.e))
            if ref_res is None:
                ref_res, ref_crs, ref_nodata = res, src.crs, src.nodata
            assert src.crs == ref_crs, f"{path}: CRS {src.crs} != {ref_crs}"
            assert abs(res[0] - ref_res[0]) < tol and abs(res[1] - ref_res[1]) < tol, \
                f"{path}: cell size {res} != {ref_res}"
            assert src.nodata == ref_nodata, f"{path}: nodata {src.nodata} != {ref_nodata}"
            # Origin must be a whole number of cells from the reference origin.
            offx = (src.bounds.left / ref_res[0]) % 1.0
            offy = (src.bounds.top / ref_res[1]) % 1.0
            assert min(offx, 1 - offx) < 1e-3 and min(offy, 1 - offy) < 1e-3, \
                f"{path}: origin is off-grid by {offx:.3f}, {offy:.3f} cells"
    return {"n": len(tifs), "res_m": ref_res[0], "crs": str(ref_crs),
            "nodata": ref_nodata}


print(assert_mosaicable(sorted(glob.glob("chm_tiles/*.tif"))))

With that passing, the mosaic itself is one command. -vrtnodata and -srcnodata must both be set or GDAL will treat NoData cells as real heights:

gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 \
    -resolution highest -r nearest \
    chm_mosaic.vrt chm_tiles/*.tif

gdal_translate -of COG -co COMPRESS=DEFLATE -co PREDICTOR=2 \
    -co BLOCKSIZE=512 -co OVERVIEWS=IGNORE_EXISTING \
    chm_mosaic.vrt chm_mosaic.tif

-r nearest is deliberate. Any other resampler will interpolate across NoData at tile edges, which manufactures a soft height ramp exactly where you are trying to prove there is no seam.

When you need custom overlap logic — taking the maximum where two halos meet, for instance, so a crown clipped by one tile is recovered from the other — do it in Python:

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


def max_of_overlaps(old_data, new_data, old_nodata, new_nodata, **kwargs):
    """rasterio.merge method: keep the taller value where tiles overlap."""
    valid = ~new_nodata
    np.copyto(old_data, np.maximum(old_data, new_data),
              where=valid & ~old_nodata)
    np.copyto(old_data, new_data, where=valid & old_nodata)


def merge_chm(tifs: list[str], out_path: str, nodata: float = -9999.0) -> None:
    srcs = [rasterio.open(p) for p in sorted(tifs)]
    try:
        mosaic, transform = merge(srcs, nodata=nodata, method=max_of_overlaps)
        profile = srcs[0].profile
        profile.update(height=mosaic.shape[1], width=mosaic.shape[2],
                       transform=transform, nodata=nodata,
                       compress="deflate", predictor=2, tiled=True)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(mosaic)
    finally:
        for s in srcs:
            s.close()

Taking the maximum is the right rule for a canopy surface specifically: a cell is only lower in one tile than another when that tile saw fewer returns there, so the taller value is the better-observed one. The same rule would be wrong for a terrain model, where the minimum is the defensible choice.

Parameter reference

Expected output and verification

Extent and grid checks prove the tiles fit together. Proving there is no value seam takes one more step: compare the distribution of horizontal height differences along the boundary columns against the same statistic in the tile interiors. A seam is a spike that exists only at the boundary.

import numpy as np
import rasterio


def seam_score(mosaic_path: str, boundary_xs: list[float], sample_rows: int = 2000) -> dict:
    """Compare |dz/dx| at tile boundaries against the interior baseline."""
    with rasterio.open(mosaic_path) as src:
        arr = src.read(1, masked=True).astype("float32")
        res = abs(src.transform.a)
        cols = [int((x - src.bounds.left) / res) for x in boundary_xs]

    step = max(1, arr.shape[0] // sample_rows)
    grad = np.abs(np.diff(arr[::step, :], axis=1))
    interior = float(np.ma.median(grad))

    edges = [float(np.ma.median(grad[:, c - 1:c + 1])) for c in cols
             if 1 <= c < grad.shape[1] - 1]
    worst = max(edges) if edges else 0.0
    return {"interior_median_dz": interior,
            "worst_boundary_median_dz": worst,
            "ratio": worst / interior if interior else float("inf")}


report = seam_score("chm_mosaic.tif", boundary_xs=[501000.0, 502000.0, 503000.0])
assert report["ratio"] < 1.5, f"seam detected at a tile boundary: {report}"

A ratio near 1.0 means the boundary looks exactly like the forest around it, which is the definition of seamless. A ratio above about 1.5 means the tiles were classified with different parameters or without a halo, and the fix is upstream — no amount of mosaicking technique will repair it.

Common pitfalls

• Setting -vrtnodata and forgetting -srcnodata. The VRT then declares a NoData value that none of its sources actually use, so the holes stay filled with −9999 and every statistic is poisoned.
• Using bilinear resampling in gdalbuildvrt. It interpolates across NoData at tile edges and creates a smooth ramp that looks like a genuine canopy gradient.
• Mosaicking before validating alignment. GDAL will happily produce output from misaligned tiles, resampling silently; the result looks fine and is blurred everywhere.
• Materialising the GeoTIFF too early. A VRT costs nothing and can be rebuilt in seconds when one tile is reprocessed. A 200 GB GeoTIFF cannot.
• Treating a seam as a cosmetic problem. A height step at a boundary propagates into crown segmentation, biomass, and fuel layers. Fix the tile processing, then re-mosaic.

Frequently Asked Questions

Related

• Canopy Height Model Creation — produces the tiles this page assembles
• Tiling and Buffering Large LiDAR Surveys with PDAL — the upstream discipline that makes seamlessness possible
• Calculating Canopy Cover from CHM in Python — a statistic that a bad seam will quietly distort
• Reading Cloud Optimized GeoTIFFs with rasterio Windows — how the finished mosaic gets read efficiently

Up: Canopy Height Model Creation · Canopy Height Modeling & Terrain Extraction