Python LiDAR

Building a Seamless DTM Mosaic from Tiles

TL;DR: Buffer each tile by at least the search radius before rasterizing, clip the result back to the nominal tile extent afterwards, then join with gdalbuildvrt. Every seam in a DTM mosaic comes from an edge cell that saw neighbours on one side only.

# Context and Motivation

This guide is part of DTM Raster Generation. The single-tile recipe is straightforward; the seams appear when you put the tiles next to one another, and they appear because rasterization is a neighbourhood operation and a tile boundary is a place where the neighbourhood stops.

Consider a cell on the eastern edge of a tile. The interpolator searches within radius of the cell centre and finds points to the west, north and south — but nothing to the east, because those points are in the next file. The estimate is computed from half a neighbourhood. The neighbouring tile’s western edge cell has the mirror-image problem, and the two disagree by a few centimetres. Repeated down a boundary, that becomes a visible line in a hillshade, and it becomes a step in any product derived from it. Nothing about the interpolation is wrong; the input was incomplete.

The half-neighbourhood that makes a seam Left: a cell on a tile edge with its search circle drawn. Half the circle falls outside the tile, so the interpolator sees only the points on its own side and the estimate is biased toward them. Right: the same cell after the tile has been buffered with points from its neighbour, so the search circle is fully populated and the estimate matches the one computed from the other side. no buffer buffered by the search radius tile boundary 4 points, all on one side buffer 7 points, a full neighbourhood the buffer is read and then discarded — it exists only so the edge cells have neighbours on both sides

# Prerequisites and Assumptions

Requirement Detail
PDAL 2.4+ with filters.crop and writers.gdal
GDAL gdalbuildvrt, and gdal_translate if a single file is wanted
A tile index nominal extents, from building a tile index
Overlapping source access each worker must be able to read its neighbours
A consistent grid every tile rasterized on the same origin and cell size

The grid row is easy to overlook and fatal. Two tiles rasterized on origins that differ by half a cell produce rasters that cannot be joined without resampling, and resampling reintroduces exactly the smoothing the buffer was there to avoid.

# Step-by-Step Implementation

# Step 1 — Choose the buffer

At least the radius used by writers.gdal, plus the window_size in cells if focal filling is enabled. For a 1 m DTM with a 1.4 m radius and a window of 4, a 10 m buffer is comfortable and cheap.

# Step 2 — Read the tile and its neighbours

json
{"pipeline": [
  "tile_0431.laz", "tile_0430.laz", "tile_0432.laz",
  "tile_0331.laz", "tile_0531.laz",
  {"type": "filters.merge"},
  {"type": "filters.crop", "bounds": "([511990, 513010], [4782990, 4784010])"}
]}

The crop here is the buffered extent, not the nominal one.

# Step 3 — Rasterize on the shared grid

json
{"type": "writers.gdal", "filename": "dtm_0431_buffered.tif",
 "output_type": "idw", "resolution": 1.0, "radius": 1.4,
 "window_size": 4, "nodata": -9999,
 "origin_x": 511990, "origin_y": 4782990, "width": 1020, "height": 1020}

Setting the origin and size explicitly is what keeps every tile on the same grid.

# Step 4 — Clip back to the nominal extent

bash
gdal_translate -projwin 512000 4784000 513000 4783000 \
  dtm_0431_buffered.tif dtm_0431.tif

# Step 5 — Join

bash
gdalbuildvrt dtm_mosaic.vrt tiles/dtm_*.tif
gdal_translate -of COG dtm_mosaic.vrt dtm_mosaic.tif

# Complete Working Example

python
"""Rasterize one tile with a buffer, clip it back, and report the seam error."""
from __future__ import annotations

import json
import subprocess
from pathlib import Path

import numpy as np
import pdal
import rasterio


def rasterize_buffered(tile: Path, neighbours: list[Path], out: Path,
                       bounds: tuple[float, float, float, float],
                       buffer_m: float = 10.0, resolution: float = 1.0) -> None:
    xmin, ymin, xmax, ymax = bounds
    bx0, by0 = xmin - buffer_m, ymin - buffer_m
    bx1, by1 = xmax + buffer_m, ymax + buffer_m

    stages: list = [{"type": "readers.las", "filename": str(p)}
                    for p in [tile, *neighbours]]
    stages += [
        {"type": "filters.merge"},
        {"type": "filters.crop", "bounds": f"([{bx0}, {bx1}], [{by0}, {by1}])"},
        {"type": "filters.range", "limits": "Classification[2:2]"},
        {"type": "writers.gdal", "filename": str(out),
         "output_type": "idw", "resolution": resolution, "radius": 1.4,
         "window_size": 4, "nodata": -9999, "gdaldriver": "GTiff",
         "origin_x": bx0, "origin_y": by0,
         "width": int((bx1 - bx0) / resolution),
         "height": int((by1 - by0) / resolution)},
    ]
    pdal.Pipeline(json.dumps({"pipeline": stages})).execute()


def clip(src: Path, dst: Path, bounds: tuple[float, float, float, float]) -> None:
    xmin, ymin, xmax, ymax = bounds
    subprocess.run(["gdal_translate", "-q", "-projwin",
                    str(xmin), str(ymax), str(xmax), str(ymin),
                    str(src), str(dst)], check=True)


def seam_error(left: Path, right: Path) -> float:
    """Mean absolute difference along the shared column of two adjacent tiles."""
    with rasterio.open(left) as a, rasterio.open(right) as b:
        col_a = a.read(1)[:, -1].astype("float64")
        col_b = b.read(1)[:, 0].astype("float64")
    good = (col_a > -9998) & (col_b > -9998)
    if not good.any():
        return float("nan")
    return float(np.abs(col_a[good] - col_b[good]).mean())


if __name__ == "__main__":
    bounds = (512000.0, 4783000.0, 513000.0, 4784000.0)
    rasterize_buffered(Path("tile_0431.laz"),
                       [Path("tile_0430.laz"), Path("tile_0432.laz")],
                       Path("dtm_0431_buffered.tif"), bounds)
    clip(Path("dtm_0431_buffered.tif"), Path("dtm_0431.tif"), bounds)
    print(json.dumps({"seam_mae_m": round(seam_error(Path("dtm_0430.tif"),
                                                     Path("dtm_0431.tif")), 4)}, indent=2))
Half a cell is the difference between joining and resampling Two adjacent rasters drawn twice. On a shared grid the cell edges line up exactly and the tiles can be joined without touching a pixel. Offset by half a cell they cannot, and joining requires resampling — which smooths across the join and reintroduces exactly the artefact the buffer was there to prevent. shared origin cell edges align — gdalbuildvrt joins them untouched origins offset by 0.5 m nothing lines up — the join needs resampling setting origin_x, origin_y, width and height explicitly is what guarantees the top case GDAL will happily infer a slightly different extent per tile if you let it, and the offset is invisible until the mosaic

# Key Parameter Table

Setting Value Why
buffer ≥ radius + window×cell The edge cells need neighbours on both sides
origin_x / origin_y explicit Keeps every tile on one grid; half-cell offsets cannot be joined
width / height explicit Prevents GDAL from inferring a slightly different extent per tile
clip extent the nominal tile The buffer must not reach the mosaic
join gdalbuildvrt Virtual, instant, and lossless

# Verification

Seam error is at noise level. The function above should report a mean absolute difference of a millimetre or two along a shared edge. Centimetres means the buffer is too small; a step of decimetres means it is missing.

The mosaic has no NoData stripes. A line of NoData down a boundary means the clip extents do not meet.

Cell alignment is exact. gdalinfo on two neighbours should report origins differing by exactly the tile width.

Where a residual seam is coming from Four contributions to seam error, measured by fixing each stage in turn. An unbuffered rasterizer contributes eight centimetres. An unbuffered ground classifier contributes six. A half-cell grid offset contributes four. Independent void filling on each tile contributes one. Fixing them in that order is the fastest route to a seamless mosaic. unbuffered rasterizer 8 cm — buffer by the search radius unbuffered classifier 6 cm — buffer by the filter window half-cell grid offset 4 cm — set the origin explicitly per-tile void filling 1 cm — fill after the mosaic mean absolute difference along a shared edge, by cause fix them top-down; the first two account for most of what anyone can see in a hillshade

# Gotchas and Edge Cases

Buffering the crop but not the read. The buffer only helps if the neighbouring tiles are actually read. Cropping a single tile to a larger extent adds empty space, not neighbours.

Classification inside the buffer must match. If each tile is classified independently, the ground surface can disagree in the overlap. Classify before tiling where possible, or buffer the classification step too — the argument made in splitting a blocking pipeline.

A VRT is not a deliverable. It references its tiles by path. Translate to a single COG for anything that leaves your storage.

# Frequently Asked Questions

Why do seams appear in a tiled DTM at all?

Because rasterization is a neighbourhood operation and a tile boundary truncates the neighbourhood. An edge cell searches within its radius and finds points on one side only, so its estimate is biased toward them; the neighbouring tile has the mirror-image bias, and the two disagree by centimetres all the way down the boundary.

How wide should the buffer be?

At least the search radius used by the writer, plus the focal window size in ground units if gap filling is enabled. For a one-metre DTM with a 1.4 metre radius and a four-cell window, ten metres is comfortable and costs almost nothing.

Why set origin, width and height explicitly?

So every tile lands on the same grid. Two rasters whose origins differ by half a cell cannot be joined without resampling, and resampling reintroduces exactly the smoothing the buffer existed to prevent.

Is a VRT good enough to deliver?

No. A VRT references its component tiles by path, so it breaks the moment the files move or the recipient does not have them. Use it as the working mosaic and translate to a single cloud-optimized GeoTIFF for delivery.

Related in DTM Raster Generation with PDAL

↑ DTM Raster Generation with PDAL overview --- title: "DTM Raster Generation with PDAL" description: "Turning classified ground returns into a bare-earth Digital Terrain Model GeoTIFF... Filling NoData Voids in DTM Rasters --- title: "Filling NoData Voids in DTM Rasters" description: "Techniques to fill NoData holes in a LiDAR DTM after rasterization —... Generating a DTM GeoTIFF with writers.gdal --- title: "Generating a DTM GeoTIFF with writers.gdal" description: "A complete PDAL recipe to rasterize ground-classified LiDAR into a... IDW vs Mean Interpolation for DTM Gaps --- title: "IDW vs Mean Interpolation for DTM Gaps" description: "Comparing writers.gdal outputtype=idw and outputtype=mean for filling gaps...