Python LiDAR

Computing RMSEz Against Survey Checkpoints

TL;DR: With ΔZ = Z_lidar − Z_checkpoint for each checkpoint, RMSEz = sqrt(mean(ΔZ²)). Always report it alongside the mean error (bias) and the standard deviation, because RMSEz² ≈ bias² + sd²: a large RMSEz with a small standard deviation is a systematic offset, usually a datum problem, not noisy data. A bootstrap over checkpoints gives a confidence interval that tells readers how much to trust the number.

# Context and Motivation

This guide is part of Vertical Accuracy Assessment for LiDAR. RMSEz is the single number every accuracy report leads with, and it is easy to compute. It is also easy to misread. The same 0.15 m RMSEz can mean “random noise of 15 cm, nothing systematic” or “a 14 cm vertical offset across the whole project with only 5 cm of noise” — and those call for completely different responses. The first is a property of the sensor and processing; the second is almost always fixable, typically by correcting a geoid model or datum realization.

Computing the decomposition takes three extra lines, and a bootstrap confidence interval takes a few more. Together they turn a single number into an interpretable result.

Same RMSEz, different stories Two right triangles whose hypotenuse is RMSEz of 0.15 metres. In dataset A the horizontal leg, bias, is 0.01 metres and the vertical leg, standard deviation, is 0.15 metres: random error. In dataset B the bias leg is 0.14 metres and the standard deviation leg is 0.05 metres: a systematic offset. The relation is RMSE squared equals bias squared plus standard deviation squared. sd 0.15 bias 0.01 A: random error RMSEz 0.15 m bias 0.14 sd 0.05 B: systematic offset RMSEz 0.15 m RMSE² = bias² + sd²

# Prerequisites and Assumptions

  • A table of checkpoints with LiDAR elevations already interpolated at each one, as produced by the vertical accuracy workflow or interpolating LiDAR elevations at checkpoints.
  • Both elevations in the same vertical datum, geoid model and units.
  • A land-cover label per checkpoint so open and vegetated terrain can be separated.
  • Python with NumPy and pandas.

# Step-by-Step Implementation

# Step 1 — Fix the sign convention

Define ΔZ = Z_lidar − Z_checkpoint and state it in the report. Positive values mean the LiDAR surface is above the surveyed ground.

# Step 2 — Compute the three statistics

Mean error (bias), sample standard deviation, and RMSEz, for open-terrain checkpoints.

# Step 3 — Check the decomposition

Verify that RMSEz² is close to bias² + sd² × (n − 1)/n. If bias dominates, stop and investigate the datum chain before reporting.

# Step 4 — Bootstrap a confidence interval

Resample checkpoints with replacement a few thousand times and take the 2.5th and 97.5th percentiles of the resampled RMSEz.

# Step 5 — Report per group

Repeat for each land-cover class, flight block or region if the project spans several; a single project-wide number can hide a bad block.

# Complete Working Example

python
"""RMSEz with bias/sd decomposition and a bootstrap confidence interval."""
from __future__ import annotations

import numpy as np
import pandas as pd


def rmse_stats(dz: np.ndarray, n_boot: int = 5000, seed: int = 0) -> dict:
    dz = np.asarray(dz, dtype=float)
    n = dz.size
    bias = dz.mean()
    sd = dz.std(ddof=1)
    rmse = np.sqrt(np.mean(dz ** 2))
    rng = np.random.default_rng(seed)
    boot = np.sqrt(np.mean(rng.choice(dz, size=(n_boot, n), replace=True) ** 2, axis=1))
    lo, hi = np.percentile(boot, [2.5, 97.5])
    return {
        "n": n,
        "bias_m": round(bias, 3),
        "sd_m": round(sd, 3),
        "rmsez_m": round(rmse, 3),
        "rmsez_ci95_m": (round(lo, 3), round(hi, 3)),
        "bias_share": round(bias ** 2 / rmse ** 2, 2) if rmse > 0 else 0.0,
        "t_bias": round(bias / (sd / np.sqrt(n)), 2) if sd > 0 else float("inf"),
    }


if __name__ == "__main__":
    res = pd.read_csv("qa/checkpoint_results.csv")          # id, cover, z_check, z_lidar, status
    res = res[res.status == "ok"].copy()
    res["dz"] = res.z_lidar - res.z_check

    for cover, grp in res.groupby("cover"):
        s = rmse_stats(grp.dz.to_numpy())
        print(cover, s)
        if s["bias_share"] > 0.5 and abs(s["t_bias"]) > 3:
            print(f"  {cover}: bias explains {s['bias_share']:.0%} of RMSE² — check the vertical datum")

Example output for a project with a geoid mismatch:

text
open {'n': 42, 'bias_m': 0.121, 'sd_m': 0.047, 'rmsez_m': 0.13, 'rmsez_ci95_m': (0.117, 0.142), 'bias_share': 0.87, 't_bias': 16.68}
  open: bias explains 87% of RMSE² — check the vertical datum
vegetated {'n': 28, 'bias_m': 0.162, 'sd_m': 0.118, 'rmsez_m': 0.2, 'rmsez_ci95_m': (0.162, 0.241), 'bias_share': 0.66, 't_bias': 7.26}

After converting the checkpoints to the same geoid model as the LiDAR, the open-terrain result typically collapses to a bias near zero and an RMSEz close to the standard deviation.

How certain is the RMSEz? A histogram of 5,000 bootstrap RMSEz values for 42 open checkpoints, centred near 0.130 metres. The central 95 percent, from 0.117 to 0.142 metres, is shaded. A note says fewer checkpoints widen the interval, so reports should state the count and interval alongside the value. 0.117 0.142 95 % interval n = 42 bootstrap RMSEz, metres

# Key Parameter Table

Statistic Formula Tells you
Mean error (bias) mean(ΔZ) Systematic offset; datum, geoid, calibration
Standard deviation sd(ΔZ), ddof = 1 Random scatter around the bias
RMSEz sqrt(mean(ΔZ²)) Combined error; the headline number
bias share bias² / RMSEz² How much of RMSEz is systematic
t statistic of bias bias / (sd / √n) Whether the bias is distinguishable from zero
bootstrap 95 % interval percentiles of resampled RMSEz Uncertainty in the headline number

# Verification

  • Hand-check three rows. Recompute ΔZ for three checkpoints by hand from the source survey file; sign mistakes are common.
  • Identity check. RMSEz² should equal bias² + sd² × (n − 1)/n to rounding.
  • Stability to one point. Remove each checkpoint in turn and recompute RMSEz; if one point moves it by more than 10–20 percent, report that point explicitly.
python
dz = res[res.cover == "open"].dz.to_numpy()
n = dz.size
assert np.isclose(np.mean(dz**2), dz.mean()**2 + dz.var(ddof=1) * (n - 1) / n)
loo = np.array([np.sqrt(np.mean(np.delete(dz, i) ** 2)) for i in range(n)])
print("most influential checkpoint changes RMSEz by", round(np.abs(loo - np.sqrt(np.mean(dz**2))).max(), 3), "m")

# Gotchas and Edge Cases

Mixed vertical datums. Checkpoints delivered as ellipsoid heights and LiDAR in orthometric heights differ by the geoid undulation — tens of metres in some places, and never zero. Even two geoid models (for example GEOID12B and GEOID18 in the US) can differ by several centimetres.

Units. Survey files in US survey feet against LiDAR in metres produce a ΔZ that scales with elevation. A bias that grows with height is the tell-tale sign.

Vegetated checkpoints in RMSEz. Pooling vegetated and open checkpoints inflates RMSEz and violates the assumptions behind the 95 % NVA factor. Keep them separate.

A bias that grows with height A scatter of checkpoint error against checkpoint elevation. The points lie along a rising straight line through the origin rather than scattering around zero, meaning the error is proportional to elevation. That pattern indicates a unit mismatch such as feet against metres, not a constant datum offset. checkpoint elevation ΔZ error ∝ elevation: a units problem

Too few checkpoints. With ten checkpoints, the bootstrap interval is wide and one outlier dominates. Report the interval; it is the honest way to say the sample is small.

# Frequently Asked Questions

How do I calculate RMSEz for LiDAR?

Subtract each surveyed checkpoint elevation from the LiDAR elevation interpolated at the same location, square the differences, average them and take the square root. Use open-terrain checkpoints for the non-vegetated figure.

What is the difference between RMSEz and standard deviation?

Standard deviation measures scatter around the mean error; RMSEz measures scatter around zero, so it also includes any systematic bias. When bias is near zero they are almost equal; when bias is large, RMSEz is much larger.

My RMSEz is high but the errors look consistent. What does that mean?

A consistent offset — a large mean error with small scatter — usually comes from a datum, geoid model or unit mismatch between checkpoints and LiDAR. Fix the reference frame before drawing conclusions about data quality.

Should I report a confidence interval for RMSEz?

It is good practice, especially with fewer than about 50 checkpoints. A bootstrap interval shows how much the headline number could move with a different sample of checkpoints.

Related in Vertical Accuracy Assessment for LiDAR

↑ Vertical Accuracy Assessment for LiDAR overview --- title: "Vertical Accuracy Assessment for LiDAR" description: "Measure and report LiDAR vertical accuracy in Python: surveyed... Interpolating LiDAR Elevations at Checkpoints --- title: "Interpolating LiDAR Elevations at Checkpoints" description: "Get the LiDAR elevation at a surveyed checkpoint the way accuracy... Reporting NVA and VVA Accuracy --- title: "Reporting NVA and VVA Accuracy" description: "Compute and report non-vegetated and vegetated vertical accuracy for LiDAR... Measuring Swath-to-Swath Relative Accuracy --- title: "Measuring Swath-to-Swath Relative Accuracy" description: "Check internal consistency between overlapping LiDAR flightlines...