3D structure from laser ranging—measuring topography and vegetation
2026-02-27
The discovery in 2010 that LiDAR surveys from aircraft could detect the outlines of Maya cities hidden beneath the jungle canopy of the Yucatán Peninsula changed how archaeologists thought about the ancient lowland Maya. Conventional survey methods, even with experienced field teams cutting transects through dense vegetation, had missed entire city centres. LiDAR sees through the canopy gaps to the ground, recording returns from both the top of the vegetation and the terrain surface beneath. The difference — the canopy height model — revealed pyramid mounds, causeways, agricultural terraces, and reservoirs that were invisible to aerial photography and impractical to find on foot. A single airborne LiDAR survey, flown in hours, mapped more ruins than decades of ground survey.
LiDAR measures distance by timing the travel of laser pulses: emit a pulse, wait for the return, multiply the elapsed time by the speed of light and halve it for the one-way distance. An airborne system fires up to a million pulses per second, combining range measurements with GPS-derived position and IMU-derived orientation to place each return as a 3D point in an absolute coordinate system. The result — the point cloud — may contain tens or hundreds of points per square metre, capturing the vertical structure of vegetation and the shape of terrain beneath it with centimetre-scale precision. Processing a point cloud requires separating ground returns from vegetation returns (a classification problem), interpolating the ground points to a smooth Digital Terrain Model, and computing the Canopy Height Model as the difference between the top-of-canopy surface and the terrain. This model derives the ranging physics, implements the ground classification algorithm, and builds the derived products.
What’s the forest canopy height and how does ground elevation vary beneath it?
LiDAR (Light Detection and Ranging) measures distance by timing laser pulse returns.
How it works: 1. Emit laser pulse (typically 905 or 1064 nm) 2. Pulse reflects off surface 3. Measure round-trip time 4. Calculate distance: d = ct/2 5. Combine with position (GPS) and orientation (IMU) to get 3D coordinates
Key advantages: - Active sensor (day/night operation) - Penetrates vegetation canopy (multiple returns) - Centimeter-level vertical accuracy - Dense sampling (1-100 points/m²)
Platforms: - Airborne (aircraft, helicopter): Large area mapping - UAV/drone: Flexible, lower cost - Mobile (vehicle-mounted): Corridor mapping - Terrestrial (tripod): Ultra-high density, limited extent
Applications: - High-resolution DEMs (bare earth topography) - Forest structure (canopy height, biomass) - Urban 3D models (buildings, infrastructure) - Archaeology (detect subtle features under vegetation) - Coastal mapping (beaches, dunes) - Powerline corridor management
Basic principle:
d = \frac{c t}{2}
Where: - d = distance to target (m) - c = speed of light (3 × 10⁸ m/s) - t = round-trip time (s) - Factor of 2: Light travels to target and back
Example:
Target 100 m away:
t = \frac{2d}{c} = \frac{200}{3 \times 10^8} = 6.67 \times 10^{-7} \text{ s} = 667 \text{ ns}
Requires nanosecond timing precision.
Laser pulse encounters multiple surfaces:
Example - forest: 1. Top of canopy (first return) 2. Mid-canopy branches (intermediate returns) 3. Ground (last return)
Discrete return systems:
Record up to 4-7 returns per pulse.
Full-waveform systems:
Record entire return signal, extract arbitrary number of returns.
Use: - First return → canopy height map (DSM - Digital Surface Model) - Last return → ground elevation (DTM - Digital Terrain Model) - All returns → vegetation structure
Each point has: - X, Y, Z coordinates (3D position) - Intensity (return signal strength) - Return number (1st, 2nd, …, last) - Classification (ground, vegetation, building, …) - RGB color (if fused with imagery)
Typical dataset:
1 km² at 10 points/m² = 10 million points
File formats: - LAS/LAZ (binary, compressed) - ASCII (text, large files)
Definition:
\text{CHM} = \text{DSM} - \text{DTM}
Where: - DSM = Digital Surface Model (top surface, first returns) - DTM = Digital Terrain Model (bare ground, last returns)
Result: Height above ground for each location.
Applications: - Forest inventory (tree height, volume, biomass) - Habitat mapping - Fire fuel load estimation
Sensor-centric coordinates → ground coordinates
Range and angles: - r = measured range - \theta = scan angle - \phi = roll angle
Sensor frame to ground frame:
\begin{bmatrix} X \\ Y \\ Z \end{bmatrix} = \begin{bmatrix} X_0 \\ Y_0 \\ Z_0 \end{bmatrix} + \mathbf{R} \begin{bmatrix} r\sin\theta\cos\phi \\ r\sin\theta\sin\phi \\ r\cos\theta \end{bmatrix}
Where: - (X_0, Y_0, Z_0) = GPS position of sensor - \mathbf{R} = rotation matrix from IMU (orientation)
Errors propagate:
GPS error (~5 cm vertical), IMU error (~0.01°), range error (~1 cm) combine.
Total vertical accuracy: Typically 5-15 cm (airborne), 1-5 cm (terrestrial).
Challenge: Identify which points are ground vs. vegetation/buildings.
Progressive morphological filter:
Algorithm: 1. Create initial grid at coarse resolution 2. Find minimum elevation in each cell (likely ground) 3. Fit surface to minima 4. Calculate residuals for all points 5. Points within threshold of surface → ground 6. Repeat at finer resolutions
Threshold selection:
t = s \times w + z_{threshold}
Where: - s = slope-dependent factor - w = window size - z_{threshold} = base threshold (~0.15 m)
Steep terrain: Larger threshold needed.
Cloth Simulation Filter (CSF):
Alternative method: 1. Invert point cloud (flip upside down) 2. Simulate cloth falling onto inverted surface 3. Points contacting cloth = ground
Advantages: Robust in complex terrain.
Height percentiles:
For all points in area (e.g., 20m plot): - Sort by height above ground - Extract percentiles (25th, 50th, 75th, 95th)
Example: - 25th percentile = 5 m (understory) - 50th percentile = 12 m (mid-canopy) - 95th percentile = 25 m (canopy top)
Canopy cover:
\text{Cover} = \frac{N_{\text{veg}}}{N_{\text{total}}} \times 100\%
Where: - N_{\text{veg}} = returns above threshold (e.g., 2m) - N_{\text{total}} = all returns
Leaf Area Index (LAI) proxy:
Related to return density and penetration ratio.
\text{LAI} \approx -\ln\left(\frac{N_{\text{ground}}}{N_{\text{total}}}\right)
Tree segmentation:
Problem: Calculate canopy height from LiDAR returns.
Point cloud data (simplified, 10 returns):
| Point | X (m) | Y (m) | Z (m) | Return |
|---|---|---|---|---|
| 1 | 100 | 200 | 345.2 | 1st |
| 2 | 100.5 | 200 | 342.8 | 2nd |
| 3 | 100 | 200 | 338.5 | Last |
| 4 | 101 | 201 | 344.9 | 1st |
| 5 | 101 | 201 | 338.7 | Last |
| 6 | 102 | 202 | 346.5 | 1st |
| 7 | 102.5 | 202 | 343.1 | 2nd |
| 8 | 102 | 202 | 339.2 | Last |
| 9 | 103 | 203 | 338.9 | 1st |
| 10 | 103 | 203 | 338.8 | Last |
Note: Points 9-10 are ground returns (no canopy).
Calculate ground elevation, canopy height, and metrics.
Step 1: Identify ground points
Last returns typically hit ground (if penetration occurs).
Ground elevations (last returns): - Point 3: 338.5 m - Point 5: 338.7 m - Point 8: 339.2 m - Point 10: 338.8 m
Also point 9 (first return = last return → open area).
Step 2: Estimate ground surface
Average ground elevation:
Z_{\text{ground}} = \frac{338.5 + 338.7 + 339.2 + 338.8 + 338.9}{5} = 338.8 \text{ m}
(In practice, would interpolate spatially)
Step 3: Calculate canopy heights
Height above ground = Z - Z_ground:
Step 4: Summary statistics
Canopy heights: 6.4, 4.0, 6.1, 7.7, 4.3 m
Step 5: Canopy cover
Vegetation points (height > 2m): 5 points
Total points: 10
\text{Cover} = \frac{5}{10} \times 100\% = 50\%
Summary: - Ground elevation: ~338.8 m - Canopy height range: 4-7.7 m - Mean canopy height: 5.7 m - Canopy cover: 50%
Interpretation: Moderate canopy (6-8 m typical for small trees/shrubs), 50% open.
Below is an interactive LiDAR point cloud analyzer.
<label>
Terrain type:
<select id="terrain-type">
<option value="flat-forest">Flat with forest</option>
<option value="slope-forest">Sloped with forest</option>
<option value="urban">Urban (buildings)</option>
<option value="bare">Bare ground</option>
</select>
</label>
<label>
Point density (pts/m²):
<input type="range" id="point-density" min="1" max="20" step="1" value="10">
<span id="density-val">10</span>
</label>
<label>
View:
<select id="view-type">
<option value="profile">Cross-section profile</option>
<option value="chm">Canopy Height Model</option>
<option value="intensity">Intensity map</option>
</select>
</label>
<div class="lidar-info">
<p><strong>Total points:</strong> <span id="total-points">--</span></p>
<p><strong>Ground points:</strong> <span id="ground-points">--</span> (<span id="ground-pct">--</span>%)</p>
<p><strong>Mean ground elev:</strong> <span id="mean-ground">--</span> m</p>
<p><strong>Max canopy height:</strong> <span id="max-height">--</span> m</p>
<p><strong>Canopy cover:</strong> <span id="canopy-cover">--</span>%</p>
</div><canvas id="lidar-canvas" width="700" height="400" style="border: 1px solid #ddd;"></canvas>Observations: - Green points show canopy/building tops (first returns) - Light green shows intermediate vegetation returns - Brown points show ground surface (last returns) - Forest scenarios show clear vertical stratification - Urban buildings show flat tops and sharp transitions - Point density affects detail captured - Bare ground shows all returns at surface - Multiple returns enable seeing through vegetation to ground
Key insights: - LiDAR penetrates canopy to measure ground beneath forest - Multiple returns reveal 3D vegetation structure - Dense point clouds enable centimeter-scale terrain mapping - Classification separates ground from non-ground features - Combination of first and last returns produces DSM and DTM
Plot-level inventory:
Traditional field measurement: Labor-intensive, sample-based.
LiDAR approach: - 100% coverage - Tree-level metrics automatically - Consistent, objective
Biomass estimation:
Allometric relationships:
\text{Biomass} = a \times H^b \times \text{Crown Area}^c
Where parameters calibrated with field plots.
Accuracy: ±10-15% for plot-level biomass.
Commercial value:
Timber volume estimation without expensive field crews.
High-resolution DEMs reveal subtle features:
Landslides: - Scarps (elevation discontinuities) - Hummocky deposits - Lateral margins
Active faults: - Offset streams - Scarps - Deformation patterns
Glacial features: - Moraines - Drumlins - Meltwater channels
Example - Puget Lowland, Washington:
Airborne LiDAR through forest canopy reveals: - Unmapped faults - Thousands of landslides - Seismic hazard reassessment
Jungle-covered sites:
Traditional surveys miss features under dense vegetation.
LiDAR penetrates canopy:
Example - Angkor, Cambodia: - Massive urban complex revealed - Irrigation networks - Previously unknown temples
Example - Maya cities, Central America: - Pyramids, platforms, causeways - Agricultural terraces - Entire city layouts
Non-invasive: No excavation required for mapping.
High-resolution DTMs improve hydraulic modelling:
Bare earth elevation essential for: - Flow direction - Inundation extent - Flood depth
1m LiDAR DTM vs 30m SRTM:
Order of magnitude better accuracy in flood prediction.
Critical infrastructure:
Levee heights measured to ±5cm enables failure risk assessment.
Near-infrared laser (905, 1064 nm) absorbed by water.
Problem:
No returns from water surface (appears as data gaps).
Impact: - Rivers/lakes unmapped - Intertidal zones missing - Flooded areas invisible
Solution: - Bathymetric LiDAR (green laser, 532 nm penetrates water) - Combine with other data sources - Mark as “no data” rather than incorrect elevation
Some canopies too dense for penetration:
Result: Ground points sparse or absent.
Impact: DTM interpolation errors.
Solution: - Higher pulse density - Multiple flight lines (different angles) - Acknowledge limitation in data gaps
Off-nadir points have: - Lower point density - Reduced penetration - Geometric distortion
Typical maximum scan angle: ±15-20° from nadir
Solution: - Flight planning (adequate overlap) - Angle-dependent filtering
Automated ground classification fails at: - Steep slopes - Sharp terrain breaks - Low vegetation (difficult to distinguish from ground) - Cars, boulders (classified as ground incorrectly)
Impact: - DTM artifacts - Bridges included in terrain (should be excluded)
Solution: - Manual editing for critical areas - Multiple classification algorithms - Validation against field data
Discrete return systems: Record peak locations only.
Full-waveform: Records entire return signal.
Advantages:
More information: - Weak returns detectable (undergrowth) - Return width (target roughness) - Multiple returns per object
Canopy density:
Waveform integral proportional to leaf area encountered.
Ground in dense canopy:
Gaussian decomposition extracts ground return from composite signal.
Challenge: Data volume (100× larger than discrete).
Research applications: - Detailed forest structure - Biomass estimation - Snow depth (dual-wavelength)
Position vector:
\mathbf{p} = \begin{bmatrix} x \\ y \\ z \end{bmatrix}
Distance between points:
d = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2 + (z_2 - z_1)^2}
Alternative representation:
\begin{align} x &= r \sin\theta \cos\phi \\ y &= r \sin\theta \sin\phi \\ z &= r \cos\theta \end{align}
Where: - r = radial distance (range) - \theta = zenith angle (from vertical) - \phi = azimuth angle
LiDAR measures (r, \theta, \phi) then converts to (x,y,z).
Orientation from IMU:
\mathbf{R} = \mathbf{R}_z(\psi) \mathbf{R}_y(\theta) \mathbf{R}_x(\phi)
Where \psi, \theta, \phi = yaw, pitch, roll.
Transform point:
\mathbf{p}_{\text{ground}} = \mathbf{R} \mathbf{p}_{\text{sensor}} + \mathbf{t}
Where \mathbf{t} = translation (GPS position).