LiDAR and Point Clouds at Scale: PDAL, COPC, and Web Delivery

National elevation programmes have quietly produced the largest open geospatial datasets in existence. Australia’s Elvis holdings, the United States Geological Survey’s 3D Elevation Program, England’s Environment Agency survey, and the Netherlands’ national height model each run to tens or hundreds of terabytes of Light Detection and Ranging (LiDAR) point data, free to download.

Almost nobody uses them well. The typical workflow is to download a few tiles, open them in a desktop viewer, generate a hillshade, and never touch the points again. The points are where the value is: building footprints, canopy structure, powerline clearance, floor levels, change over time. All of it is derivable, and the tooling to do it at scale is mature and open.

# Formats: A Short Lineage

Four formats matter, and they form a clear progression.

LAS is the American Society for Photogrammetry and Remote Sensing binary standard. Every point has coordinates, intensity, a return number, and a classification code. It is uncompressed and consequently enormous.

LAZ is LAS compressed, typically to about one-seventh the size, with no loss. It has been the de facto distribution format for fifteen years. Its weakness is that it is a stream: to read the last point you decompress everything before it.

EPT, the Entwine Point Tile format, solved that by exploding a cloud into an octree of many small files plus a JSON index. It works, and it is a poor fit for object storage — a national dataset becomes millions of small files.

COPC, Cloud-Optimised Point Cloud, is the current answer and it is genuinely a step change. It is a valid LAZ file with an octree index embedded in its variable-length records. One file. Any LAZ reader can open it. A COPC-aware client can issue HTTP range requests to fetch only the octree nodes covering the area and level of detail it needs.

The parallel with Cloud-Optimised GeoTIFF is exact, and so is the consequence: the point cloud tile server era is over. Static file hosting plus a range-capable client is the whole architecture.

# PDAL: Pipelines, Not Scripts

PDAL is to point clouds what the Geospatial Data Abstraction Library (GDAL) is to rasters. Its distinguishing feature is that processing is expressed as declarative JSON — a reviewable, version-controllable artefact rather than a script whose behaviour you have to read to understand.

A pipeline is a list of stages. The first reads, the last writes, and everything between is a filter:

{
  "pipeline": [
    "input.laz",
    { "type": "filters.assign", "assignment": "Classification[:]=0" },
    { "type": "filters.elm" },
    { "type": "filters.outlier", "method": "statistical",
      "mean_k": 12, "multiplier": 2.2 },
    { "type": "filters.smrf", "ignore": "Classification[7:7]",
      "slope": 0.2, "window": 16, "threshold": 0.45, "scalar": 1.2 },
    { "type": "filters.hag_nn" },
    { "type": "writers.copc", "filename": "output.copc.laz" }
  ]
}

Executed with pdal pipeline classify.json. Reading it in order: reset any existing classification, flag extreme low noise, remove statistical outliers, classify ground using the Simple Morphological Filter, compute Height Above Ground for every point, and write a cloud-optimised archive.

Two of those stages deserve attention.

filters.smrf is the ground classifier, and its parameters are the ones you will actually tune. window is the largest non-ground object size in metres — set it to comfortably exceed your biggest building footprint. slope is the maximum terrain slope you expect as a ratio. threshold is the vertical tolerance in metres. On flat urban ground slope: 0.15, window: 18 works well; in steep forest slope: 0.5, window: 12 is closer.

filters.hag_nn computes Height Above Ground by nearest-neighbour interpolation of the ground returns. This single derived dimension unlocks most downstream analysis, because “how high is this point above the earth beneath it” is the question almost every application actually asks.

# Deriving Products

Once classified, rasters fall out directly. A bare-earth terrain model uses only ground returns:

pdal translate output.copc.laz dtm.tif \
  --readers.copc.filename=output.copc.laz \
  -f range --filters.range.limits="Classification[2:2]" \
  -w writers.gdal --writers.gdal.resolution=1.0 \
  --writers.gdal.output_type="idw" --writers.gdal.window_size=4

A canopy height model uses the Height Above Ground dimension on vegetation returns:

pdal translate output.copc.laz chm.tif \
  -f range --filters.range.limits="Classification[3:5]" \
  -w writers.gdal --writers.gdal.dimension="HeightAboveGround" \
  --writers.gdal.output_type="max" --writers.gdal.resolution=1.0

And building footprints come from clustering high non-ground returns:

{
  "pipeline": [
    "output.copc.laz",
    { "type": "filters.range", "limits": "Classification[6:6]" },
    { "type": "filters.cluster", "min_points": 40, "tolerance": 1.2 },
    { "type": "filters.groupby", "dimension": "ClusterID" },
    { "type": "filters.hexbin", "edge_size": 0.8, "threshold": 4 },
    { "type": "writers.ogr", "filename": "buildings.gpkg" }
  ]
}

Powerline clearance is the same shape of problem: filter to classification 14 (wire conductor), buffer the conductors, and query vegetation returns within the buffer whose Height Above Ground exceeds the clearance envelope. The utilities that pay consultants six figures for this analysis are paying for the specification, not the computation.

# Building a National Index

Where the work becomes interesting is at survey scale. A state-wide LiDAR holding is thousands of tiles, and the pattern that makes it tractable has three parts.

First, build a tile index so you can find the right files without opening them:

pdal tindex create --tindex holdings.gpkg \
  --filespec "/data/lidar/**/*.laz" \
  --fast_boundary -f GPKG

Second, convert each tile to COPC. This is embarrassingly parallel and a textbook case for the High Performance Computing (HPC) job array pattern described in Spatial Processing on HPC:

find /data/lidar -name '*.laz' | \
  parallel -j 16 'pdal translate {} /data/copc/{/.}.copc.laz \
    -w writers.copc --writers.copc.forward=all'

Third, query by geometry rather than by filename:

import geopandas as gpd, subprocess, json

idx = gpd.read_file("holdings.gpkg")
aoi = gpd.read_file("study_area.gpkg").to_crs(idx.crs)
needed = idx[idx.intersects(aoi.union_all())]
print(f"{len(needed)} of {len(idx)} tiles intersect the study area")

for path in needed.location:
    subprocess.run(["pdal", "pipeline", "classify.json",
                    "--readers.las.filename", path])

The storage arithmetic is worth stating plainly, because it determines whether a project is viable. A one-square-kilometre urban tile at 8 points per square metre is about 8 million points, roughly 100 MB as LAZ. A 40,000 km² region is therefore around 4 TB. National coverage at that density runs to 40–80 TB. On commodity object storage that is a few hundred dollars a month; on a managed platform with per-gigabyte processing fees it is not a project you will get approved.

# Delivering Points to a Browser

Three viable clients, for three different jobs.

Potree is the specialist. It is a mature WebGL point cloud viewer with measurement tools, clipping volumes, profile extraction, and elevation-based colour ramps. If your users’ job is to inspect and measure the cloud itself, use Potree and do not build anything.

deck.gl is the integrator. Its PointCloudLayer and Tile3DLayer render points alongside your other data layers in a single MapLibre canvas, which matters when the point cloud is context for something else:

import { Tile3DLayer } from '@deck.gl/geo-layers';
import { COPCLoader } from '@loaders.gl/las';

const layer = new Tile3DLayer({
  id: 'lidar',
  data: 'https://cdn.example.com/survey/site14.copc.laz',
  loader: COPCLoader,
  pointSize: 1.4,
  getPointColor: (d) => {
    const hag = d.HeightAboveGround ?? 0;
    return hag < 0.5  ? [110, 110, 110]      // ground
         : hag < 3.0  ? [ 60, 160,  70]      // low vegetation
         : hag < 12.0 ? [ 30, 110,  45]      // canopy
                      : [220, 180,  60];     // structures
  },
  onTilesetLoad: (t) => console.log(`${t.tiles.length} nodes`),
});

copc.js is the primitive. Where you need programmatic access — reading a specific octree node to answer a query rather than to draw it — it gives you range-request reads directly:

import { Copc, Getter } from 'copc';

const getter = Getter.http('https://cdn.example.com/survey/site14.copc.laz');
const copc = await Copc.create(getter);
const pages = await Copc.loadHierarchyPage(getter, copc, copc.info.rootHierarchyPage);
const node = pages.nodes['2-1-3-0'];              // depth-x-y-z
const view = await Copc.loadPointDataView(getter, copc, node);
const getZ = view.getter('Z');
console.log(`node holds ${view.pointCount} points`);

For all three, the delivery infrastructure is the same: the .copc.laz file on object storage, correct Cross-Origin Resource Sharing (CORS) headers, and Accept-Ranges: bytes enabled. Cloudflare R2 does this by default and charges nothing for egress to its edge. There is no server.

# Where This Fits

Point clouds are the acquisition layer for most of the geometry in a geospatial digital twin, the source for the terrain models that feed windowed raster processing, and the input to change detection when you have two epochs of the same area.

That last case is the one most worth pursuing and least often attempted. Two surveys of the same region, both classified with the same pipeline, differenced on Height Above Ground, will show you vegetation growth, subsidence, erosion, construction, and demolition — quantitatively, over a whole region, from data you did not pay for. The processing is a day’s work with the pipelines above. The reason it rarely happens is that the two epochs were delivered in different projections at different densities with different classification schemes, and nobody budgeted for the normalisation.

Budget for the normalisation.