⚙️ How-To LiDAR & point clouds Terrain & elevation LiDAR & Point Clouds · step 14 of 14 Updated

How to Measure Building Heights from LiDAR in Python

Problem statement

Building height from lidar sounds like a lookup: take the building footprint, take the maximum return inside it, subtract the ground. Four things make it harder.

The building class is often absent. A real 3DEP survey used here classifies 22.9% ground and 77.1% unclassified — no building class at all. Any code filtering classification == 6 returns nothing.

The ground under a building is never observed. The laser cannot see through a roof, so the DTM there is interpolated from the surroundings.

"Height" is ambiguous. Eaves height, ridge height, mean roof height and maximum height are different numbers, and a pitched roof separates them by metres.

Footprints and returns disagree at the edge. A footprint is a wall outline; a roof overhangs it, and mixed pixels at the boundary contain both roof and ground.

Quick answer

Work from normalised heights and a robust statistic, buffered inward:

import numpy as np
import geopandas as gpd
from shapely.geometry import Point


def building_heights(x, y, height_above_ground, footprints,
                     inset_m=1.0, min_points=10, percentile=95):
    """Height per footprint, from returns inside an inward-buffered outline."""
    points = gpd.GeoDataFrame(
        {"h": height_above_ground},
        geometry=gpd.points_from_xy(x, y), crs=footprints.crs)

    inner = footprints.copy()
    inner["geometry"] = inner.geometry.buffer(-inset_m)
    inner = inner[~inner.geometry.is_empty]

    joined = gpd.sjoin(points, inner[["geometry"]], predicate="within")
    stats = joined.groupby("index_right")["h"].agg(
        n="size", p50="median",
        pxx=lambda s: np.percentile(s, percentile), pmax="max")

    stats = stats[stats["n"] >= min_points]
    print(f"  {len(stats)} of {len(footprints)} footprints have "
          f"{min_points}+ returns after a {inset_m} m inset")
    return stats

The inward buffer removes the mixed edge pixels. The percentile rather than the maximum removes the single spurious high return — a bird, an aerial, a chimney.

A pitched-roof building showing eaves height, mean roof height, 95th-percentile height and maximum height as four different numbers.
Four defensible heights for one building. Say which you computed.

Step-by-step solution

1. Normalise the point cloud first 

height_above_ground = z - dtm[row, col]

Every building height is a height above ground, and the ground under a building is interpolated. That interpolation is the largest source of error in the whole calculation — see the explanation below.

2. Do not rely on a building class

Measured on the survey used here:

class 1 unclassified   77.13%
class 2 ground         22.87%
class 7 low noise       0.00%

If a building class exists, use it. If not, use footprints from a vector source — OpenStreetMap, a national mapping agency, a cadastre — and select returns geometrically.

3. Buffer the footprint inward

A footprint is the wall line at ground level. The roof usually overhangs it, and any return within a metre of the edge may be roof, wall, gutter or ground.

An inward buffer of 0.5–1.5 m removes the ambiguous ring. On small buildings it can remove the entire polygon, so drop empty geometries and record how many.

4. Choose the statistic, and name it 

maximum          the highest return: a chimney, an aerial, a bird
95th percentile  robust to single outliers, close to the ridge
median           roughly mid-roof; on a flat roof, the roof
minimum          eaves, or the ground if the inset was too small

For a pitched roof the median and the maximum can differ by several metres. For a flat roof they nearly coincide. Report which you used, and prefer a percentile to the raw maximum.

5. Require a minimum number of returns

A footprint with three returns gives a "height" that is one measurement with no redundancy. At 16.9 points per square metre, a 100 m² building should have around 1,700 returns — an order of magnitude fewer means something is wrong with the footprint, the alignment or the flight coverage.

A building footprint with a ring of mixed roof-and-ground returns at its edge, removed by an inward buffer.
The edge ring contains roof, wall and ground. An inward buffer is cheaper than trying to classify it.

Code examples

Example 1 — heights with the uncertainty carried through 

import numpy as np
import geopandas as gpd


def height_table(x, y, z, dtm, transform, cell, footprints,
                 inset_m=1.0, min_points=10):
    """Per-building heights plus the evidence: counts, spread, ground support."""
    left, top = transform.c, transform.f
    h, w = dtm.shape
    col = np.clip(((x - left) / cell).astype(int), 0, w - 1)
    row = np.clip(((top - y) / cell).astype(int), 0, h - 1)
    ground_z = dtm[row, col]
    above = z - ground_z

    ok = np.isfinite(above)
    points = gpd.GeoDataFrame(
        {"above": above[ok], "ground": ground_z[ok]},
        geometry=gpd.points_from_xy(x[ok], y[ok]), crs=footprints.crs)

    inner = footprints.reset_index(drop=True).copy()
    inner["geometry"] = inner.geometry.buffer(-inset_m)
    lost = int(inner.geometry.is_empty.sum())
    inner = inner[~inner.geometry.is_empty]

    joined = gpd.sjoin(points, inner[["geometry"]], predicate="within")
    grouped = joined.groupby("index_right")

    table = grouped["above"].agg(
        returns="size",
        h_median="median",
        h_p95=lambda s: float(np.percentile(s, 95)),
        h_max="max",
        h_iqr=lambda s: float(np.percentile(s, 75) - np.percentile(s, 25)),
    )
    table["ground_range"] = grouped["ground"].agg(lambda s: float(s.max() - s.min()))
    table["area_m2"] = inner.geometry.area
    table["returns_per_m2"] = table["returns"] / table["area_m2"]

    enough = table["returns"] >= min_points
    print(f"  {len(footprints)} footprints, {lost} removed by the "
          f"{inset_m} m inset, {int(enough.sum())} with {min_points}+ returns")
    print(f"  median height p95 {table.loc[enough, 'h_p95'].median():.1f} m, "
          f"median IQR {table.loc[enough, 'h_iqr'].median():.1f} m")

    steep = table["ground_range"] > 2.0
    if steep.any():
        print(f"  ! {int(steep.sum())} buildings sit on ground varying by "
              "more than 2 m — height depends on which ground you mean")
    return table[enough]

The ground_range column is the one that stops silly answers. A building on a slope has no single ground level, so its "height" depends on whether you measure from the uphill or downhill side — and the difference can exceed the building's own height.

Example 2 — deriving footprints when you have none 

import numpy as np
from scipy.ndimage import binary_opening, binary_closing, label


def candidate_buildings(chm, cell=1.0, min_height=2.5, min_area_m2=30.0,
                        max_roughness=1.5):
    """Flat, tall, compact blobs in a CHM. A first pass, not a product."""
    tall = np.isfinite(chm) & (chm > min_height)
    tall = binary_closing(binary_opening(tall, np.ones((3, 3))), np.ones((3, 3)))

    labels, n = label(tall)
    keep = np.zeros_like(labels, bool)
    accepted = 0

    for i in range(1, n + 1):
        blob = labels == i
        area = blob.sum() * cell ** 2
        if area < min_area_m2:
            continue
        values = chm[blob]
        roughness = float(np.percentile(values, 90) - np.percentile(values, 10))
        if roughness > max_roughness:        # vegetation: rough on top
            continue
        keep |= blob
        accepted += 1

    print(f"  {n} blobs above {min_height} m, {accepted} kept "
          f"(area >= {min_area_m2} m², roughness <= {max_roughness} m)")
    return keep, labels

The roughness test is what separates roofs from trees: a roof is a plane or two planes, so its height distribution is narrow, while a tree crown is rough over its own diameter. It fails on pitched roofs with a steep pitch and on dense flat-topped hedges — which is why this is a first pass rather than a building layer.

Example 3 — a sanity check against known heights 

import numpy as np


def validate_heights(table, known, height_col="h_p95", known_col="height_m"):
    """Compare against surveyed or attributed heights, if you have any."""
    merged = table.join(known[[known_col]], how="inner")
    if merged.empty:
        print("  no overlap with the reference heights")
        return None

    residual = merged[height_col] - merged[known_col]
    print(f"  {len(merged)} buildings with a reference height")
    print(f"  bias {residual.mean():+.2f} m, "
          f"RMSE {np.sqrt((residual ** 2).mean()):.2f} m")
    print(f"  p5 {residual.quantile(0.05):+.2f} m, "
          f"p95 {residual.quantile(0.95):+.2f} m")

    if abs(residual.mean()) > 1.0:
        print("  ! systematic bias — check the height definition and the "
              "vertical datum on both sides")
    return residual

A systematic bias of a metre or two almost always means a definition mismatch — your 95th percentile against their eaves height — or a vertical datum difference. Neither is a lidar problem, and both are invisible without a reference.

Explanation

Why the ground under a building is the biggest error

A roof is opaque, so no pulse reaches the ground beneath it. Every DTM cell inside a footprint is interpolated from the ground around the building.

On flat terrain that interpolation is excellent — the ground under a house on a level street really is at street level. On a slope it is a guess, and the guess is a plane fitted across the building.

The measured proxy is the ground_range column: the spread of interpolated ground heights under one footprint. Where that exceeds a couple of metres, "the height of the building" is not a well-defined quantity, and any single number needs a stated convention (uphill, downhill or mean).

Why the maximum is the wrong statistic

The maximum return inside a footprint is the highest thing the laser saw there. That is often a chimney, a satellite dish, a lightning conductor or a bird.

A single spurious return sets the height, and there is no redundancy at all. The 95th percentile keeps the ridge while discarding the top 5% of returns, which on a building with a thousand returns is fifty points — enough to absorb every aerial.

For a flat roof the two are nearly identical, which is why the maximum survives in code that was only ever tested on flat roofs.

Why footprints and lidar disagree at the edge

A cadastral footprint is the wall line at ground level. The roof overhangs it by anything from zero to a metre, and gutters project further.

So the ring just inside a footprint contains roof, and the ring just outside also contains roof. Returns in that band are a mixture, and no threshold separates them cleanly.

An inward buffer discards the ambiguous band. It costs area — and on a small building it can cost the whole polygon, which is why the count of removed footprints belongs in the output.

Why a height class would not solve it

Even with a building class present, all four problems remain: the ground under the roof is still interpolated, the height definition is still ambiguous, edge returns are still mixed, and the class itself has errors at exactly the boundaries where it matters.

Classification helps with selection. It does not make the measurement well-defined, and the definition is the part that most often makes two organisations' building heights disagree.

Interpolated ground under a building on flat terrain against a slope, where the height depends on which side you measure from.
On a slope, "the height of the building" is not a well-defined quantity without a stated convention.

Edge cases or notes

  • Do not assume a building class exists. This survey has none.
  • Buffer footprints inward by 0.5–1.5 m and count what that removes.
  • Use a percentile, not the maximum. The maximum is a chimney or a bird.
  • Name the height definition — eaves, mean roof, p95, maximum.
  • Report the interpolated ground spread under each footprint.
  • Require a minimum return count; at 16.9 points per square metre a 100 m² roof should have around 1,700.
  • Buildings on slopes have no single height. State the convention.
  • Check the vertical datum on both sides before comparing with reference heights.

FAQ

How do I measure building height from lidar?

Normalise the point cloud against a DTM, select returns inside an inward-buffered footprint, and take a high percentile of the height above ground.

What if my point cloud has no building class?

Use vector footprints and select returns geometrically. Many surveys classify only ground — the one measured here is 77% unclassified.

Should I use the maximum height?

No. The maximum is whatever the highest object was — a chimney, an aerial, a bird. The 95th percentile keeps the ridge and discards single outliers.

Why buffer the footprint inward?

Because a footprint is the wall line and the roof overhangs it. The ring at the edge contains a mixture of roof, wall and ground returns.

How accurate is the ground under a building?

It is interpolated, because no pulse reaches it. Good on flat terrain, a guess on a slope. Report the spread of interpolated ground under each footprint.

Why do my heights disagree with the council's?

Usually a definition mismatch — eaves against ridge, or mean roof against maximum — or a vertical datum difference. Check both before blaming the lidar.

How many returns do I need per building?

At 16.9 points per square metre, a 100 m² roof should have around 1,700. An order of magnitude fewer means a footprint, alignment or coverage problem.