Computing Building Height from Lidar-Derived Surfaces

Height limits are among the easiest zoning standards to state and the hardest to measure, because “height” in an ordinance is a distance between two things the code defines and lidar does not: a grade reference at the bottom and a roof reference at the top. A normalised surface gives a number; whether that number is the height the ordinance means depends entirely on which references were used. This guide computes building height from a DSM and a DTM, implements the three grade definitions codes actually use, and reports the uncertainty — because a measurement compared against a limit needs to say when it is too close to call. It is the vertical half of height and FAR compliance logic.

Prerequisites

Step-by-step

Step 1: Confirm the two surfaces agree before subtracting them

A normalised height surface is DSM minus DTM, which is only meaningful if both are on the same grid, the same vertical datum and the same vertical unit. Mixed vertical datums are common — NAVD88 for one product and an ellipsoidal height for the other — and the difference in the contiguous United States runs to tens of metres, which produces heights that are absurd rather than subtly wrong. The subtler failure is a vertical unit mismatch, which produces heights wrong by 3.28 and entirely plausible.

import rasterio

def check_surfaces(dsm_path, dtm_path):
    with rasterio.open(dsm_path) as dsm, rasterio.open(dtm_path) as dtm:
        if dsm.crs != dtm.crs:
            raise ValueError(f"CRS mismatch: {dsm.crs} vs {dtm.crs}")
        if dsm.transform != dtm.transform or dsm.shape != dtm.shape:
            raise ValueError("grids are not aligned; resample to a common grid first")
        return {"crs": dsm.crs.to_string(), "res": dsm.res, "nodata": (dsm.nodata, dtm.nodata)}

The vertical datum is not carried in most raster headers, so it has to come from the product metadata and be recorded by hand in the run manifest. Where it is unavailable, treat the height as advisory rather than determinative — a point worth stating explicitly in the report rather than discovering during an appeal.

What must agree before two surfaces can be subtractedGrid alignment, horizontal CRS, vertical datum and vertical unit all have to match; only the first two are carried in a raster header.Same grid and transformchecked from the file; resample first if notSame horizontal CRSchecked from the fileSame vertical datumNOT in the header — from product metadata, by handSame vertical unitoften differs from the horizontal unit
A vertical unit mismatch produces heights wrong by 3.28 and entirely plausible — which is worse than the datum error that produces absurd ones.

Step 2: Implement the grade reference the ordinance names, not the convenient one

This is where automated height measurement most often disagrees with a surveyor, and the disagreement is definitional rather than technical. Codes use three grade references and they give materially different answers on a sloping site:

import numpy as np
from rasterstats import zonal_stats

def grade_reference(footprint, dtm_path, definition="average"):
    """Average grade, lowest adjacent grade, or the mean of the corner elevations."""
    if definition == "average":
        s = zonal_stats(footprint, dtm_path, stats=["mean"], all_touched=False)[0]
        return s["mean"]
    if definition == "lowest_adjacent":
        ring = footprint.buffer(1.5).difference(footprint)
        s = zonal_stats(ring, dtm_path, stats=["min"])[0]
        return s["min"]
    if definition == "corner_mean":
        pts = [footprint.exterior.interpolate(t, normalized=True) for t in (0, .25, .5, .75)]
        vals = [next(rasterio.open(dtm_path).sample([(p.x, p.y)]))[0] for p in pts]
        return float(np.mean(vals))
    raise ValueError(f"unknown grade definition {definition}")

On a level lot the three agree to within a few inches. On a site sloping ten feet across the footprint, lowest-adjacent grade produces a height five feet greater than average grade — enough to turn a compliant building into a violation, using the same lidar and the same roof. The definition is a property of the ordinance and belongs in the rule pack, next to the limit itself.

Step 3: Sample the roof robustly rather than taking the maximum

The naive roof reference is the maximum DSM value inside the footprint, and it is almost always wrong. Lidar returns from an antenna, a chimney, a lift overrun or a tree overhanging the roof all land inside the footprint and all exceed the roof. Most codes exclude such appurtenances explicitly.

def roof_reference(footprint, dsm_path, quantile=0.98, inset=1.0):
    """A high quantile of the inset footprint: above the roof plane, below the aerials."""
    core = footprint.buffer(-inset)
    if core.is_empty:
        core = footprint
    s = zonal_stats(core, dsm_path, stats=["percentile_98", "max", "count"],
                    add_stats=None)[0]
    return {"roof": s["percentile_98"], "max": s["max"], "cells": s["count"]}

The negative inset removes the edge cells, which straddle the wall and mix roof returns with ground returns; the quantile removes the appurtenances. Keeping the raw maximum alongside is worth doing, because a large gap between the 98th percentile and the maximum is a reliable signal that something — a tree, most often — is sitting over the footprint and the footprint may be stale.

The two ends of the measurement, and how codes define themGrade references and roof references compared by what they measure and where they disagree with a surveyor.What it takesWhere it goes wrongAverage gradeMean DTM inside the footprintUnderstates height on a slopeLowest adjacentgradeMinimum DTM in a ring outsideSensitive to a single low cellRoof as DSMmaximumThe highest return insideChimneys, aerials, overhanging treesRoof as an insetquantile98th percentile, footprint insetNeeds enough cells to be stable
On a site sloping ten feet across the footprint, lowest-adjacent grade gives a height five feet greater than average grade — same lidar, same roof.

Step 4: Compute the height and carry its uncertainty

The height is the difference, and the difference of two uncertain quantities is more uncertain than either. Lidar vertical accuracy is usually quoted as an RMSE in the region of 10 cm for the bare earth and rather worse on complex roofs; propagating it is a single line and it changes how the verdict should be phrased.

DSM_RMSE_FT = 0.5     # from the lidar report, not from optimism
DTM_RMSE_FT = 0.33

def building_height(footprint, dsm_path, dtm_path, definition="average"):
    grade = grade_reference(footprint, dtm_path, definition)
    roof = roof_reference(footprint, dsm_path)
    h = roof["roof"] - grade
    sigma = (DSM_RMSE_FT ** 2 + DTM_RMSE_FT ** 2) ** 0.5
    return {"height_ft": float(h), "grade_ft": float(grade),
            "roof_ft": float(roof["roof"]), "sigma_ft": float(sigma),
            "grade_definition": definition, "roof_max_ft": float(roof["max"])}

Step 5: Compare against the limit with a band, not a point

A measurement with an uncertainty of six-tenths of a foot cannot distinguish 34.8 feet from a 35-foot limit. Reporting that as a violation is indefensible; reporting it as compliant is equally arbitrary. The honest answer is three-valued.

def height_verdict(measure: dict, limit_ft: float, k: float = 2.0):
    h, s = measure["height_ft"], measure["sigma_ft"]
    if h - k * s > limit_ft:
        return "exceeds"
    if h + k * s < limit_ft:
        return "complies"
    return "indeterminate — within measurement uncertainty of the limit"
A measurement with uncertainty against a hard limitLidar height carries an uncertainty of several tenths of a foot, so a building close to the limit cannot be decided from the surface alone.Is the height further fromthe limit than twice itsuncertainty?noReturn indeterminate and route for surveya confident verdict here is not supportableyesIssue the verdictcomplies or exceeds, with the marginRecord the grade definition, roof rule, datum and lidar date either waya building finished after the flight has no roof in the data
The indeterminate band is not a hedge — it is the measurement naming the parcels that need a surveyed height.

The indeterminate band is not a hedge; it is the measurement telling you which parcels need a surveyed height. On a typical municipality it is a small percentage of buildings, and identifying them is far more valuable than issuing confident verdicts on all of them. The same three-valued pattern applies wherever a measured quantity meets a hard threshold, including the area measurements in implementing FAR checks with Shapely and GeoPandas.

Step 6: Record enough to defend the number

Height is the standard most likely to be contested, so the record needs to reconstruct the measurement completely.

record = {
    **measure,
    "dsm_source": dsm_manifest["content_sha256"],
    "dtm_source": dtm_manifest["content_sha256"],
    "lidar_collected": "2024-04-11",
    "vertical_datum": "NAVD88 (Geoid18)",
    "vertical_unit": "us-ft",
    "footprint_source": "county building footprints 2025-Q1",
    "roof_rule": "98th percentile of footprint inset 1.0 ft",
}

The lidar collection date is the field people forget and the first one an applicant asks about, because a building completed after the flight has no roof in the data at all — a case that shows up as a height of roughly zero and should be caught rather than reported.

Verification

Check that the surfaces are consistent, that heights are physically plausible, and that the footprint vintage matches the lidar.

h = heights["height_ft"]
assert (h > -1.0).all(), "negative height — DSM and DTM are swapped or datums differ"

no_building = h < 3.0
print(f"{no_building.sum()} footprints with no structure in the DSM "
      f"(demolished, or built after the {LIDAR_DATE} flight)")

implausible = h > 400.0
assert not implausible.any(), "heights above 400 ft — check the vertical unit"

# Ground-truth against surveyed heights where any exist.
resid = surveyed["height_ft"] - heights.loc[surveyed.index, "height_ft"]
print(f"vs surveyed: bias {resid.mean():+.2f} ft, RMSE {(resid ** 2).mean() ** 0.5:.2f} ft")

The residual against surveyed heights is the check worth building even from a handful of buildings. A bias near zero validates the grade definition; a consistent bias of a few feet usually means the code’s grade reference is not the one implemented, and that is much better found in a comparison of ten buildings than in an appeal.

Common Pitfalls

  • Taking the DSM maximum as the roof. Chimneys, aerials and overhanging trees all sit inside the footprint and all exceed the roof.
  • Mixed vertical datums. An ellipsoidal DSM against an orthometric DTM yields heights wrong by tens of metres — obvious once, invisible if only one product is ever wrong.
  • A vertical unit different from the horizontal one. Common in lidar products and it produces a height wrong by 3.28 with no other symptom.
  • Using average grade where the code says lowest adjacent grade. On a sloping site this understates height by several feet, systematically, in the direction that favours the applicant.
  • Ignoring the lidar date. New construction has no roof in an old DSM and returns a height near zero, which reads as compliant.
  • Reporting a point estimate against a hard limit. Buildings within the uncertainty band need a surveyed height, and saying so is the correct output.

Frequently Asked Questions

Can I use a normalised DSM product instead of subtracting myself?

Yes, and it saves a step, but you inherit whatever grade reference the vendor used — which is almost never the ordinance’s. An nDSM measures height above local ground at each cell, and no code defines height that way. Use it for screening and compute the ordinance’s definition for anything that produces a verdict.

What resolution is needed?

One metre is comfortable for building height; half a metre is better on small structures and on complex roofs. Below about two metres per cell, the inset-plus-quantile approach starts sampling too few cells on small footprints, and the count returned in step 3 is what tells you when that is happening.

How do I handle a building that spans two parcels?

Measure the building once, then attribute the result to whichever parcel the rule applies to; a shared structure across a lot line is usually a legal question that predates the pipeline. What the pipeline must not do is measure the footprint’s portion within each parcel separately, since the clipped portion’s roof quantile is not the building’s height.

Does this work for pitched roofs?

It measures the ridge, which is what most codes require for a flat or low-slope roof and what many codes explicitly do not require for a pitched one — a common provision measures to the midpoint between eave and ridge. That is implementable from the same surface by taking a lower quantile for the eave and averaging, but the definition has to come from the code rather than from the geometry.

How should the indeterminate cases be handled operationally?

Route them for a surveyed height, and report them as pending rather than as either outcome. The volume is manageable precisely because the band is narrow, and the alternative — a confident verdict on a measurement that cannot support one — is the failure mode that damages the credibility of the whole pipeline.

Is the uncertainty really that large?

Larger than most people assume. The quoted RMSE for a lidar product is for well-defined open terrain; roofs are neither, and the footprint’s own positional error contributes as well, since a footprint offset by a metre samples the wrong cells near the edges. The inset in step 3 mitigates that, and the residual check against surveyed heights is what turns the estimate into a measured quantity for your data rather than a manufacturer’s claim.

Where does the height feed next?

Into the envelope, alongside the setbacks. Codes that step a building back as it rises need both, and the two standards interact — a building can comply with the height limit at its centre and violate a sky-exposure plane at its street wall. Handling corner lots with two front setbacks covers the footprint side of that envelope.

Part of: Height and FAR compliance logic