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. The 3DEP survey of Manhattan used here has no building class at all, over one of the densest built-up areas in the world (step 2). 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


keep = classification != 7                                   # low noise
stats = building_heights(x[keep], y[keep], height_above_ground[keep], footprints)
print(stats.describe().loc[["50%"]].round(1).to_string())

  1683 of 1684 footprints have 10+ returns after a 1.0 m inset
         n   p50   pxx  pmax
50%  469.0  17.2  18.9  20.7

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.

Normalise, skip the building class, buffer inward

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:

codes, counts = np.unique(classification, return_counts=True)
for code, count in zip(codes, counts):
    print(f"class {code:2d}  {count:10,}  {count / len(classification):6.2%}")

class  1   4,462,972  36.00%
class  2   1,642,927  13.25%
class  7       2,599   0.02%
class  9      54,390   0.44%
class 10       1,200   0.01%
class 17   5,155,039  41.58%
class 18   1,077,348   8.69%

Classes 1, 2, 7 and 9 are unclassified, ground, low noise and water. 17 and 18 are this LAS 1.2 survey's own codes; none of the seven is a building.

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. Work out what to expect from the density:

density = len(x) / ((x.max() - x.min()) * (y.max() - y.min()))
print(f"{density:.1f} points per m²: a 100 m² roof should hold about {100 * density:,.0f} returns")

5.5 points per m²: a 100 m² roof should hold about 551 returns

An order of magnitude fewer means something is wrong with the footprint, the alignment or the flight coverage.

Heights with uncertainty carried through

The data used on this page

One 1,500 m tile of the USGS 3DEP survey of New York City after Hurricane Sandy (NY_New_York_CMGP_SANDY_LiDAR_15, tile 18TWL880150 — the Upper East Side between Central Park and the East River, flown 2013–14), and the city's own building footprints from NYC Open Data over the same area. The footprints carry a roof height, in feet, measured from building plans, oblique aerial imagery and street-level imagery — not from lidar — which makes them an independent check. The setup also builds the 1 m DTM, filled under buildings from the nearest observed ground, the CHM, and each point's height above ground.

import io
import os
import subprocess

import geopandas as gpd
import laspy
import numpy as np
import requests
from pyproj import Transformer

TILE = "18TWL880150"
URL = ("https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/NY_New_York_CMGP_SANDY_LiDAR_15/"
       f"NY_CMPG_2013/LAZ/USGS_LPC_NY_New_York_CMGP_SANDY_LiDAR_15_{TILE}.laz")
for _ in range(20):                                       # rockyweb is slow and drops transfers: resume
    if os.path.exists("tile.laz"):
        break
    if subprocess.run(["curl", "-sfL", "-C", "-", "--retry", "3", "-o", "tile.laz.part", URL]).returncode == 0:
        os.replace("tile.laz.part", "tile.laz")

las = laspy.read("tile.laz")                              # NAD83(2011) / UTM 18N, EPSG:6347, NAVD88 metres
x, y, z = np.asarray(las.x), np.asarray(las.y), np.asarray(las.z)
classification = np.asarray(las.classification)

if not os.path.exists("footprints.gpkg"):                 # NYC Open Data building footprints over the tile
    to_ll = Transformer.from_crs(6347, 4326, always_xy=True)
    (w, s), (e, n) = to_ll.transform(x.min(), y.min()), to_ll.transform(x.max(), y.max())
    r = requests.get("https://data.cityofnewyork.us/resource/5zhs-2jue.geojson", timeout=300, params={
        "$where": f"within_box(the_geom, {n}, {w}, {s}, {e})", "$limit": 50000})
    gpd.read_file(io.BytesIO(r.content)).to_file("footprints.gpkg")
footprints = gpd.read_file("footprints.gpkg").to_crs(6347)
footprints["height_m"] = footprints["height_roof"].astype(float) * 0.3048       # the field is in feet
print(f"{len(x):,} points; {len(footprints):,} footprints, "
      f"{footprints['height_roof'].astype(float).gt(0).sum():,} with a roof height")
footprints = footprints.reset_index(drop=True)

# a 1 m DTM from the lowest ground return per cell; no pulse reaches the ground under
# a roof, so those cells are filled from the nearest observed ground
from rasterio.transform import from_origin
from scipy.ndimage import distance_transform_edt

CELL = 1.0
left, top = x.min(), y.max()
width, height = int(np.ceil((x.max() - left) / CELL)), int(np.ceil((top - y.min()) / CELL))
transform = from_origin(left, top, CELL, CELL)
cols = np.clip(((x - left) / CELL).astype(int), 0, width - 1)
rows = np.clip(((top - y) / CELL).astype(int), 0, height - 1)
flat = rows * width + cols
ground = classification == 2
dtm = np.full(width * height, np.inf)
np.minimum.at(dtm, flat[ground], z[ground])
dtm = np.where(np.isinf(dtm), np.nan, dtm).reshape(height, width)
observed = np.isfinite(dtm)
dtm = dtm[tuple(distance_transform_edt(~observed, return_distances=False, return_indices=True))]
dsm = np.full(width * height, -np.inf)
np.maximum.at(dsm, flat, z)
chm = np.where(np.isinf(dsm), np.nan, dsm).reshape(height, width) - dtm
height_above_ground = z - dtm[rows, cols]
print(f"ground observed in {observed.mean():.1%} of 1 m cells; the rest interpolated")

12,396,475 points; 1,684 footprints, 1,681 with a roof height
ground observed in 47.0% of 1 m cells; the rest interpolated

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]


table = height_table(x[keep], y[keep], z[keep], dtm, transform, CELL, footprints)

  1684 footprints, 0 removed by the 1.0 m inset, 1683 with 10+ returns
  median height p95 18.9 m, median IQR 1.2 m
  ! 722 buildings sit on ground varying by more than 2 m — height depends on which ground you mean

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. On the Upper East Side 722 of 1,684 footprints sat on filled ground varying by more than 2 m — partly slope, partly a nearest-neighbour fill that gives a large building the heights of different streets on different sides.

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


roofs, labels = candidate_buildings(chm)
from rasterio.features import rasterize

on_footprint = rasterize(footprints.geometry, out_shape=chm.shape, transform=transform).astype(bool)
print(f"  of the candidate cells, {on_footprint[roofs].mean():.1%} fall on a city footprint; "
      f"of the footprint cells, {roofs[on_footprint].mean():.1%} were found")

  1257 blobs above 2.5 m, 187 kept (area >= 30.0 m², roughness <= 1.5 m)
  of the candidate cells, 61.9% fall on a city footprint; of the footprint cells, 92.9% were found

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. Against the city's footprints it found 92.9% of the building cells, but only 61.9% of what it kept was on a footprint: in Manhattan whole blocks of attached buildings merge into one blob, and the gaps between them come along.

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


known = footprints[(footprints["height_m"] > 0)
                   & (footprints["construction_year"].astype(float) < 2013)]   # built before the flight
for statistic in ("h_median", "h_p95", "h_max"):
    print(statistic)
    residual = validate_heights(table, known, height_col=statistic)

h_median
  1621 buildings with a reference height
  bias -1.33 m, RMSE 6.32 m
  p5 -6.36 m, p95 +1.59 m
  ! systematic bias — check the height definition and the vertical datum on both sides
h_p95
  1621 buildings with a reference height
  bias +1.71 m, RMSE 5.91 m
  p5 -0.99 m, p95 +5.73 m
  ! systematic bias — check the height definition and the vertical datum on both sides
h_max
  1621 buildings with a reference height
  bias +4.03 m, RMSE 8.24 m
  p5 +0.18 m, p95 +12.33 m
  ! systematic bias — check the height definition and the vertical datum on both sides

A systematic bias of a metre or two almost always means a definition mismatch or a vertical datum difference. Here the city's roof height sat between two lidar statistics: the median of the returns was 1.33 m below it and the 95th percentile 1.71 m above, with the maximum 4.03 m above. Both are relative heights, so no datum is involved; the gap is definition — the city measures to the roof, and the highest returns on a Manhattan roof are likely to be parapets, bulkheads and water tanks. Neither is a lidar problem, and both are invisible without a reference.

The ground under a building is the biggest error

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 — 43% of the buildings here — "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. Measured against the city's heights, the maximum ran 4.03 m high on average and more than 12 m high for one building in twenty.

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.

Missing classes, inward buffers and percentiles

• Do not assume a building class exists. This survey of Manhattan 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 5.5 points per square metre a 100 m² roof should have around 550.
• Buildings on slopes have no single height. State the convention.
• Check the vertical datum on both sides before comparing with reference heights.

Internal links

• DSM, DTM and CHM from LiDAR: how each surface is derived — the normalisation this depends on
• How to create a DTM from LiDAR ground points in Python — the interpolated ground under the roof
• How to create a canopy height model from LiDAR — the same normalisation for vegetation
• How to filter a point cloud by class and return number — when a building class does exist
• How to calculate zonal statistics in Python — the general per-polygon summary
• How to perform a spatial join in Python (GeoPandas) — assigning returns to footprints
• My LiDAR DTM has holes, spikes or terraces — the ground surface underneath
• Vertical datums explained: why your elevations are off — comparing against reference heights
• Where building heights come from, and how wrong they are — comparing the sources

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 Manhattan survey measured here has no building class at all.

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. Against New York City's roof heights, the lidar median ran 1.3 m low and the 95th percentile 1.7 m high. Check both before blaming the lidar.

How many returns do I need per building?

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