Colorizing LiDAR from Orthophotos with PDAL

This page assigns red, green and blue values to LiDAR points from an aerial orthophoto with PDAL’s filters.colorization — matching the CRS of both datasets (EPSG:25832 here), scaling 8-bit imagery into LAS 16-bit colour, restricting colour to surfaces an orthophoto can actually see, and measuring the image-to-cloud misalignment that produces coloured roof edges on the ground.

Why you hit this

A coloured point cloud is far easier to read than one shaded by intensity or height, and many viewers and customers expect one. When the scan was flown without a camera, or its imagery is unusable, colouring from the national or municipal orthophoto is the obvious fix and a single PDAL filter does it. The result looks right from above and wrong everywhere else: walls painted with the colours of the pavement in front of them, tree crowns with a grey halo of road, roof edges repeated as a shadow on the street. Those are not bugs in PDAL; they are properties of orthophotos, and they are predictable enough to handle. The mesh-texturing equivalent of this problem is in aligning photogrammetry textures with point clouds.

Prerequisites

  • PDAL 2.6+ with Python bindings (pdal>=3.4), rasterio>=1.3, numpy>=1.24, scikit-image>=0.22.
  • LAZ in EPSG:25832 with heights in DHHN2016 (EPSG:7837), classified with at least ground (2), vegetation (3–5) and building (6).
  • An orthophoto GeoTIFF covering the same area, ideally a true orthophoto; its CRS must be known and its date close to the flight date.
  • Point data record format 7 or 8 output (LAS 1.4) so the file has RGB fields.

Step-by-Step

1. Confirm the two datasets share a CRS and overlap

python
import json
import pdal
import rasterio
from rasterio.warp import transform_bounds

with rasterio.open("dop20_32_691_5335.tif") as img:
    img_epsg = img.crs.to_epsg()
    img_bounds = img.bounds
    print("orthophoto:", img_epsg, img.res, img.dtypes, img.count, "bands")

info = pdal.Pipeline(json.dumps({"pipeline": ["tile_691_5335.laz", {"type": "filters.info"}]}))
info.execute()
meta = info.metadata["metadata"]["filters.info"]
bbox = meta["bbox"]
print("cloud:", meta["srs"]["horizontal"][:40], "…", bbox)

assert img_epsg == 25832, "reproject the orthophoto to EPSG:25832 first (gdalwarp)"
assert img_bounds.left <= bbox["minx"] and img_bounds.right >= bbox["maxx"], "orthophoto does not cover the tile in x"
assert img_bounds.bottom <= bbox["miny"] and img_bounds.top >= bbox["maxy"], "orthophoto does not cover the tile in y"

filters.colorization samples the raster at each point’s x and y in the point’s own coordinates; it does not reproject between the cloud and the image. An orthophoto delivered in a different CRS — EPSG:4258 geographic coordinates, or Gauss-Krüger EPSG:31468 in an older German archive — produces either an error or colours sampled from the wrong place entirely. Reproject the image once with gdalwarp -t_srs EPSG:25832 rather than the cloud.

2. Colourise with correct band mapping and scaling

python
pipeline = {
    "pipeline": [
        "tile_691_5335.laz",
        {
            "type": "filters.colorization",
            "raster": "dop20_32_691_5335.tif",
            "dimensions": "Red:1:256.0, Green:2:256.0, Blue:3:256.0",
        },
        {
            "type": "writers.las",
            "filename": "tile_691_5335_rgb.laz",
            "minor_version": 4, "dataformat_id": 7,
            "compression": "laszip", "forward": "all", "a_srs": "EPSG:25832+7837",
        },
    ]
}
n = pdal.Pipeline(json.dumps(pipeline)).execute()
print(f"{n:,} points colourised")

Each entry in dimensions is Name:band:scale. The LAS specification stores colour as 16-bit unsigned integers, and most viewers expect 8-bit imagery to be scaled into that range; multiplying by 256 maps 255 to 65,280. Leaving the scale at 1 produces a cloud that some viewers render almost black, because 255 out of 65,535 is 0.4% brightness. An orthophoto with a near-infrared fourth band is common — CIR products put infrared in band 1 — so check the band order in the image metadata before trusting 1, 2, 3.

Why roof colour lands on the street in a standard orthophoto A camera above and to the side of a building sees the roof displaced outward in a standard orthophoto, which was rectified to the terrain rather than to the building. The roof appears in the image over the street beside the building. Colourising ground points there gives them roof colours, and the wall points get pavement colours, while a true orthophoto rectified to a surface model places the roof correctly. camera building standard ortho: roof shown here true ortho: roof here displacement grows with height and view angle
A standard orthophoto is only geometrically correct at terrain height. Anything taller leans away from the image centre, and its colour is sampled by the wrong points.

3. Colour only what the orthophoto can see

An orthophoto is a view from above. It carries valid colour for surfaces facing up — ground, roofs, crowns — and no information at all about walls. Compute normals, and colour only points whose surface faces the sky.

python
pipeline = {
    "pipeline": [
        "tile_691_5335.laz",
        {"type": "filters.normal", "knn": 12},
        {"type": "filters.assign", "value": ["Red = 32896", "Green = 32896", "Blue = 32896"]},
        {
            "type": "filters.colorization",
            "raster": "dop20_32_691_5335.tif",
            "dimensions": "Red:1:256.0, Green:2:256.0, Blue:3:256.0",
            "where": "(Classification == 2 || Classification == 6 || Classification == 5) && NormalZ > 0.5",
        },
        {"type": "writers.las", "filename": "tile_691_5335_rgb.laz", "minor_version": 4,
         "dataformat_id": 7, "compression": "laszip", "forward": "all", "a_srs": "EPSG:25832+7837"},
    ]
}
pdal.Pipeline(json.dumps(pipeline)).execute()

The where option, available on most PDAL filters, applies the colourisation to a subset while leaving the other points untouched. Walls and low-confidence points first receive a neutral grey, so they are visibly “uncoloured” rather than silently carrying a colour from the pavement. NormalZ > 0.5 keeps surfaces within 60° of horizontal, which includes pitched roofs; filters.normal orients normals upward by default, so the sign is meaningful for airborne data.

4. Measure the image-to-cloud offset

Even a well-georeferenced orthophoto can be shifted relative to the LiDAR by several pixels. Rasterise the cloud’s intensity and compare it with the image’s luminance.

python
import numpy as np
from skimage.registration import phase_cross_correlation

res = 0.2
grid = {
    "pipeline": [
        "tile_691_5335.laz",
        {"type": "filters.range", "limits": "Classification[2:2]"},
        {"type": "writers.gdal", "filename": "intensity_020.tif", "resolution": res,
         "dimension": "Intensity", "output_type": "mean", "bounds":
         f"([{bbox['minx']},{bbox['maxx']}],[{bbox['miny']},{bbox['maxy']}])"},
    ]
}
pdal.Pipeline(json.dumps(grid)).execute()

with rasterio.open("intensity_020.tif") as a, rasterio.open("dop20_32_691_5335.tif") as ortho:
    lidar = a.read(1)
    window = rasterio.windows.from_bounds(*a.bounds, transform=ortho.transform)
    rgb = ortho.read((1, 2, 3), window=window, out_shape=(3, *lidar.shape)).astype("float32")
luma = 0.299 * rgb[0] + 0.587 * rgb[1] + 0.114 * rgb[2]

valid = np.isfinite(lidar) & (lidar > 0)
shift, error, _ = phase_cross_correlation(
    np.where(valid, lidar, 0), np.where(valid, luma, 0), upsample_factor=10
)
print(f"offset: {shift[1] * res:+.2f} m east, {-shift[0] * res:+.2f} m north (error {error:.3f})")

Ground-only intensity is compared because road markings, kerbs and field boundaries appear in both datasets without relief displacement. Phase correlation returns the shift in pixels, row first; a result within a pixel is good registration, and anything over a metre explains visible colour bleeding along every road edge. A consistent offset can be corrected by shifting the orthophoto’s geotransform before colourisation.

Run the correlation on several windows across a large sheet rather than once for the whole image. A single global shift is typical when the orthophoto and the LiDAR were referenced to different realisations of ETRS89 or processed with different geoid models; a shift that varies smoothly across the sheet points instead at the digital terrain model used to rectify the imagery, which displaces the image wherever that model disagrees with the LiDAR ground. The first case is fixed by one translation. The second cannot be fixed by moving the image and is a reason to request a true orthophoto rectified against a current surface model, or to accept coloured ground points near slopes as approximate and document it in the output metadata.

Colourisation pipeline with masking The classified LAZ gets normals, then every point is set to neutral grey. Colourisation is applied only where the class is ground, building or high vegetation and the normal faces upward. The output is written as LAS 1.4 point format 7. In parallel, a ground intensity raster is correlated with the orthophoto to measure any offset, which feeds back into the colourisation as a corrected raster. classifiedLAZ filters.normalknn = 12 assignneutral grey colorizationwhere: up-facing LAS 1.4format 7 intensity vs orthophase correlation corrected offset
The offset check runs once per orthophoto sheet and corrects the raster, so every tile colourised from that sheet benefits.

Expected Output & Verification

text
orthophoto: 25832 (0.2, 0.2) ('uint8', 'uint8', 'uint8', 'uint8') 4 bands
cloud: PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS8 … {'maxx': 692000.0, 'maxy': 5336000.0, 'minx': 691000.0, 'miny': 5335000.0, …}
18,441,902 points colourised
offset: +0.14 m east, -0.06 m north (error 0.412)

Beyond the offset, check the colour distribution by class. Ground and roof points should have a wide spread of colours; a cluster of pure black (0, 0, 0) means points fell outside the image or on its nodata border, and a large share of mid-grey on roofs means the normal threshold is excluding pitched roofs.

python
import laspy
import numpy as np

las = laspy.read("tile_691_5335_rgb.laz")
rgb = np.column_stack([las.red, las.green, las.blue])
for cls, name in ((2, "ground"), (6, "building"), (5, "high veg")):
    m = las.classification == cls
    grey = np.all(rgb[m] == 32896, axis=1).mean()
    black = np.all(rgb[m] == 0, axis=1).mean()
    print(f"{name:<9} {m.sum():>10,} pts | neutral {grey * 100:5.1f}% | black {black * 100:4.1f}%")
Share of points left neutral by class Bars show the share of points left neutral grey after masked colourisation. Ground is about two percent neutral, mostly steep embankments. Buildings are about thirty-five percent neutral, the walls. High vegetation is about eight percent neutral, the sides of crowns. Black nodata points are under a tenth of a percent in every class. ground: 2% — steep banks building: 35% — walls, as intended high vegetation: 8% — crown sides points left neutral grey after masked colourisation
A third of building points staying neutral is the mask working: those are wall points an orthophoto never observed.

Common Errors

filters.colorization: Unable to open raster. GDAL cannot read the image path, often a VRT whose relative source paths broke when the job ran from another directory. Use absolute paths inside the VRT, or pass a single GeoTIFF.

The whole cloud is nearly black in the viewer. The scale was left at 1, so 8-bit values were stored as 16-bit colour. Rerun with a scale of 256, or rescale in place with filters.assign and Red = Red * 256.

Colour stripes along one edge of every tile. The orthophoto sheet ends inside the tile and the points beyond it received black. Build a VRT mosaic of neighbouring sheets and colourise against that, or buffer the image extent assertion in step 1.

Frequently Asked Questions

Should I colourise before or after classification?

After. Classification drives the mask, and colour does not help ground filters. Keep colourisation as a late, re-runnable step so it can be repeated when a newer orthophoto arrives without touching classification.

Can I colour walls from street-level imagery instead?

Yes, but not with filters.colorization, which samples a single georeferenced raster in plan. Walls need projection from calibrated oblique or street-level cameras, which is a texture-mapping problem rather than a raster lookup.

Does the orthophoto date matter much?

For roofs and roads, little; for vegetation and anything mobile, a great deal. A summer flight coloured from a winter orthophoto gives leaf-on crowns the colour of bare branches and the ground beneath them. Record both dates in the output metadata.

Back to Texture Mapping Workflows for Digital Twins.