Planetary-scale answers, unlocked.
A Hands-On Guide for Working with Large-Scale Spatial Data. Learn more.
Authors
LiDAR (light detection and ranging) data is a set of 3D points measured by a laser scanner, which times each light pulse's round trip to a surface. Each point carries x, y, and z coordinates, the echo's strength (intensity), its order (return number), and a class code such as ground, vegetation, building, or water. Airborne LiDAR data maps terrain, tree canopy, and buildings over whole regions. The U.S. Geological Survey (USGS) describes LiDAR data as the basis for elevation models with a vertical accuracy of about 10 cm. This article tests one claim with Portland data. The claim is that a 30 m radar map of the surface, minus the LiDAR ground height, gives the mean height of a forest and fails on towers.
The National Oceanic and Atmospheric Administration (NOAA) defines LiDAR as a remote sensing method that uses light in the form of a pulsed laser to measure ranges to the Earth. The range is half the round-trip time multiplied by the speed of light. Each measurement becomes one point.
The raw product is a point cloud: a loose set of points, with no grid and no lines between them. One survey samples every surface the pulses reached: ground, trees, roofs, power lines, and parked cars. Photogrammetry also builds point clouds, from overlapping photographs. LiDAR data adds several returns per pulse, which reach the ground through gaps in a canopy.
The laser scanner measures the range and angle of each pulse. A GNSS (Global Navigation Satellite System) receiver records the scanner's position by satellite, and an inertial measurement unit records its roll, pitch, and heading. Wehr and Lohr laid out these components and the positioning math in Airborne laser scanning: an introduction and overview (1999). In NOAA's description, topographic LiDAR uses near-infrared light to map land, and bathymetric LiDAR uses green light that penetrates water to map riverbeds and the seafloor.
Returns. A pulse spreads into a footprint, so part of the beam hits the canopy top, part hits branches, and part reaches the ground.
Discrete return vs full waveform. Most survey sensors record discrete returns, a few peaks per pulse. Full-waveform sensors record the whole returned signal, including the structure between peaks. Mallet and Bretar reviewed the method in Full-waveform topographic lidar: state-of-the-art (2009).
Point density. Density sets what a survey can resolve. The USGS quality levels require at least 2 pulses per square meter for QL2 and at least 8 for QL1. Both require 10 cm vertical accuracy on open ground, measured as RMSEz, the root mean square error of heights. A QL1 survey resolves roof planes and individual crowns that QL2 blurs.
NASA's ICESat-2 launched in September 2018. Its ATLAS (Advanced Topographic Laser Altimeter System) laser fires 10,000 times a second at 532 nm across six beams and counts single photons. Markus and colleagues set out its design in The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2): science requirements, concept, and implementation (2017). GEDI was installed on the International Space Station in December 2018. Its three full-waveform lasers produce eight ground tracks of 25 m footprints. Dubayah and colleagues described it in The Global Ecosystem Dynamics Investigation: high-resolution laser ranging of the Earth's forests and topography (2020).
LAS is the binary exchange format for point clouds. It is owned by the American Society for Photogrammetry and Remote Sensing (ASPRS) and maintained in the open. The Open Geospatial Consortium (OGC) has recognized LAS as a community standard since 2018. Each point record holds x, y, z, intensity, return number, number of returns, scan angle, GPS time, and a classification code. LAS 1.5, approved in August 2025, drops the legacy point formats 0 to 5 and keeps formats 6 to 10.
LAZ is LAS with lossless compression. Martin Isenburg's LASzip (2013) produces files 7 to 25 percent of the original size. Cloud Optimized Point Cloud (COPC) is a LAZ 1.4 file with its points sorted into nested cubes (an octree) and an index of where each cube sits. A client can download one area or level of detail over HTTP without the whole file.
The LAS specification defines the ASPRS standard point classes. These are the values for point formats 6 to 10:
Codes 8 and 12 and the range 23 to 63 are reserved, and 64 to 255 are free for user definitions.
Most analysts use LiDAR data after it is gridded into rasters, usually GeoTIFF:
The main public LiDAR archives are free:
WherobotsDB reads LiDAR-derived rasters directly and queries LiDAR points as tables.
Derived rasters. 3DEP publishes its 1 m DEMs as Cloud Optimized GeoTIFFs in a public bucket, and RS_FromPath opens them as out-db rasters. The file stays in the bucket, and a query reads only the tiles it needs:
SELECT RS_MetaData(r) AS md, RS_BandPixelType(r) AS pt, RS_BandNoDataValue(r) AS nd, RS_SRID(r) AS srid FROM (SELECT RS_FromPath('s3://prd-tnm/StagedProducts/Elevation/1m/Projects/OR_PortlandMetro_B24/TIFF/USGS_1M_10_x52y505_OR_PortlandMetro_B24.tif') AS r)
The tile from the 2024 Portland Metro survey holds 10,012 by 10,012 pixels of 1 m heights, stored as 32-bit floating-point numbers (float32). Its spatial reference identifier (SRID) is 26910, which is NAD83 (North American Datum of 1983) on Universal Transverse Mercator (UTM) zone 10N. It uses 512-pixel internal tiles, and −999999 marks missing values (NoData). A canopy height model or DSM stored as GeoTIFF loads the same way. RS_MapAlgebra subtracts one from the other in a single statement.
Point tables. The Wherobots documentation lists no LAS or LAZ reader, so point clouds are converted first. The Point Data Abstraction Library (PDAL) has an Arrow writer that writes LAS points to GeoParquet, one row per point. That file loads as a Havasu table, the WherobotsDB table format, with an ST_PointZ geometry for each point. ST_3DDWithin and ST_3DDistance then find points within a 3D distance of each other. A spatial join with ST_Intersects assigns every point to a parcel or building footprint.
The walkthrough uses Portland, Oregon, where these sources overlap:
OR_PortlandMetro_B24
x52y505
x52y504
WherobotsDB reads the rasters and the vector polygons in one SQL statement; most stacks pick one data model.
The Copernicus DEM is a surface model that includes trees and buildings. Surface minus ground measures how much height a 30 m radar model records over forest and towers.
The first query takes four Overture land use polygons: three forested parks in the West Hills and the Albina rail yard as open ground. It computes the mean of each raster inside each polygon, a zonal statistic. RS_ZonalStatsAll returns the sum and count of pixels per raster tile, so a zone spanning several tiles is weighted by its pixels in each.
-- Lidar ground, radar surface, and lidar-trained canopy height for four Portland zones WITH zones AS ( SELECT names.primary AS zone, ST_SetSRID(geometry, 4326) AS geom, ST_AreaSpheroid(geometry) / 10000 AS ha FROM wherobots_open_data.overture_maps_foundation.base_land_use WHERE bbox.xmax > -122.76 AND bbox.xmin < -122.66 AND bbox.ymax > 45.48 AND bbox.ymin < 45.56 AND names.primary IN ('Washington Park', 'Hoyt Arboretum', 'Marquam Nature Park', 'Albina Yard') ), dtm AS ( -- USGS 3DEP 1 m bare-earth DEM from the 2024 Portland Metro lidar survey (NAD83 / UTM 10N, NAVD88) SELECT RS_FromPath(concat('s3://prd-tnm/StagedProducts/Elevation/1m/Projects/OR_PortlandMetro_B24/TIFF/USGS_1M_10_', t, '_OR_PortlandMetro_B24.tif')) AS rast FROM VALUES ('x52y504'), ('x52y505') AS v(t) ), ground AS ( SELECT /*+ BROADCAST(d) */ z.zone, RS_ZonalStatsAll(d.rast, z.geom, 1, false, true) AS s FROM zones z CROSS JOIN dtm d ), surface AS ( -- Copernicus GLO-30, a radar digital surface model (EGM2008 heights) SELECT z.zone, RS_ZonalStatsAll(c.rast, z.geom, 1, false, true) AS s FROM zones z JOIN wherobots_open_data.copernicus_dem.glo_30m c ON ST_Intersects(c.footprint, z.geom) ), canopy AS ( -- Meta and WRI canopy height model, metres, trained on airborne lidar SELECT z.zone, RS_ZonalStatsAll(m.rast, z.geom, 1, false, true) AS s FROM zones z JOIN wherobots_open_data.meta_canopy_height.global_v2 m ON ST_Intersects(m.footprint, z.geom) ), agg AS ( SELECT 'ground' AS layer, zone, SUM(s.sum) AS total, SUM(s.count) AS n, MAX(s.max) AS mx FROM ground WHERE s.count > 0 GROUP BY zone UNION ALL SELECT 'surface', zone, SUM(s.sum), SUM(s.count), MAX(s.max) FROM surface WHERE s.count > 0 GROUP BY zone UNION ALL SELECT 'canopy', zone, SUM(s.sum), SUM(s.count), MAX(s.max) FROM canopy WHERE s.count > 0 GROUP BY zone ) SELECT z.zone, ROUND(z.ha, 1) AS ha, ROUND(MAX(CASE WHEN layer = 'ground' THEN total / n END), 2) AS lidar_ground_m, ROUND(MAX(CASE WHEN layer = 'surface' THEN total / n END), 2) AS radar_surface_m, ROUND(MAX(CASE WHEN layer = 'surface' THEN total / n END) - MAX(CASE WHEN layer = 'ground' THEN total / n END), 2) AS surface_minus_ground_m, ROUND(MAX(CASE WHEN layer = 'canopy' THEN total / n END), 2) AS mean_canopy_m, MAX(CASE WHEN layer = 'canopy' THEN mx END) AS max_canopy_m, MAX(CASE WHEN layer = 'ground' THEN n END) AS ground_px, MAX(CASE WHEN layer = 'surface' THEN n END) AS surface_px, MAX(CASE WHEN layer = 'canopy' THEN n END) AS canopy_px FROM zones z JOIN agg a ON a.zone = z.zone GROUP BY z.zone, z.ha ORDER BY mean_canopy_m DESC
It ran in 144 seconds. Washington Park holds 1,963,217 LiDAR ground pixels and 2,935 radar pixels.
Over the forested parks, radar surface minus LiDAR ground lands within 0.3 to 1.9 m of the canopy model's mean. Together, a 30 m surface model and a 1 m LiDAR terrain model give a forest's average height. They miss the tallest conifers, 49 to 50 m in the canopy model, which a 30 m cell averages with lower trees. The rail yard, with almost no canopy, reads −0.95 m. That baseline combines radar noise with the gap between two vertical datums, the reference surfaces that heights are measured from: NAVD88 for 3DEP and EGM2008 for GLO-30.
The second query repeats the subtraction for the downtown buildings listed in the sources table. It groups buildings by the source of their Overture height: USGS LiDAR, OpenStreetMap, or Microsoft's machine-learning heights. It compares each height with the highest radar pixel in the footprint minus the mean LiDAR ground under it. Each raster is clipped once to the box and copied to every worker, so each building reads two small rasters in memory. The query ran in 142 seconds.
-- Radar surface minus lidar ground at downtown Portland buildings, by Overture height source WITH aoi AS ( SELECT ST_SetSRID(ST_PolygonFromEnvelope(-122.695, 45.514, -122.670, 45.532), 4326) AS box ), b AS ( -- Footprints of 2,000 m2 or more with a height, fully inside the box SELECT id, height, CASE WHEN exists(sources, s -> s.dataset = 'USGS Lidar') THEN 'USGS Lidar' WHEN exists(sources, s -> s.property = '/properties/height' AND s.dataset = 'Microsoft ML Buildings') THEN 'Microsoft ML' ELSE 'OpenStreetMap' END AS height_source, ST_SetSRID(geometry, 4326) AS geom FROM wherobots_open_data.overture_maps_foundation.buildings_building WHERE bbox.xmin > -122.695 AND bbox.xmax < -122.670 AND bbox.ymin > 45.514 AND bbox.ymax < 45.532 AND height IS NOT NULL AND ST_AreaSpheroid(geometry) >= 2000 ), -- Clip each raster once to the box. RS_Clip(raster, band, geom, allTouched, noDataValue, crop) dtm AS ( -- USGS 3DEP 1 m bare-earth DEM from the 2024 Portland Metro lidar survey (NAVD88 heights) SELECT /*+ BROADCAST(a) */ RS_Clip(RS_FromPath(concat( 's3://prd-tnm/StagedProducts/Elevation/1m/Projects/OR_PortlandMetro_B24/TIFF/USGS_1M_10_', v.t, '_OR_PortlandMetro_B24.tif')), 1, a.box, true, -999999, true) AS rast FROM VALUES ('x52y504'), ('x52y505') AS v(t) CROSS JOIN aoi a ), dsm AS ( -- Copernicus DEM GLO-30, a radar digital surface model (EGM2008 heights) SELECT /*+ BROADCAST(a) */ RS_Clip(c.rast, 1, a.box, true, -32767, true) AS rast FROM wherobots_open_data.copernicus_dem.glo_30m c JOIN aoi a ON ST_Intersects(c.footprint, a.box) ), zs AS ( -- Zonal statistics per building on the small clipped rasters, broadcast to every partition SELECT /*+ BROADCAST(d) */ 'ground' AS layer, b.id, RS_ZonalStatsAll(d.rast, b.geom, 1, false, true) AS s FROM b CROSS JOIN (SELECT rast FROM dtm WHERE rast IS NOT NULL) d UNION ALL SELECT /*+ BROADCAST(d) */ 'surface', b.id, RS_ZonalStatsAll(d.rast, b.geom, 1, false, true) FROM b CROSS JOIN (SELECT rast FROM dsm WHERE rast IS NOT NULL) d ), per_building AS ( -- Highest radar surface pixel in the footprint minus the mean lidar ground under it SELECT b.id, b.height_source, b.height AS overture_height_m, MAX(CASE WHEN z.layer = 'surface' THEN z.s.max END) - SUM(CASE WHEN z.layer = 'ground' THEN z.s.sum END) / SUM(CASE WHEN z.layer = 'ground' THEN z.s.count END) AS radar_minus_ground_m FROM b JOIN zs z ON z.id = b.id WHERE z.s.count > 0 GROUP BY b.id, b.height_source, b.height ) -- One row per height source plus an all-buildings row (height_source is NULL). -- SELECT * FROM per_building returns each building. SELECT height_source, COUNT(*) AS buildings, ROUND(corr(overture_height_m, radar_minus_ground_m), 2) AS r, ROUND(percentile(overture_height_m, 0.5), 1) AS median_overture_m, ROUND(percentile(radar_minus_ground_m - overture_height_m, 0.5), 1) AS median_gap_m, COUNT_IF(overture_height_m >= 100) AS towers_100m, ROUND(percentile(CASE WHEN overture_height_m >= 100 THEN radar_minus_ground_m END, 0.5), 1) AS towers_radar_m FROM per_building WHERE radar_minus_ground_m IS NOT NULL GROUP BY GROUPING SETS ((height_source), ()) ORDER BY buildings DESC
It returned 117 buildings:
For low and mid-rise buildings whose Overture height came from USGS LiDAR, the radar surface tracks the height (r = 0.55). It sits a median 4.8 m low, consistent with 30 m cells that mix roofs with streets. The method fails on towers. A tower covers one or two 30 m cells. Synthetic aperture radar looks at the ground from the side, so tall facades are shifted in the image and hide what is behind them. The six towers of 100 m or more, all with OpenStreetMap heights, read a median 10.8 m above the LiDAR ground. The 34 Microsoft machine-learning (ML) heights, median 15.7 m, correlate little with the radar surface (r = 0.08). Overture's LiDAR-sourced heights came from 3DEP point clouds.
The same SQL runs per 3DEP tile across a county or a state. The Model Context Protocol (MCP) is a standard way to connect an AI agent to tools and data. An agent connected through the Wherobots MCP Server can write it from a plain-language request. The models people use every day were trained on text, documents, databases, and the internet. Spatial SQL over elevation rasters and footprints gives those models measured heights of real places, the basis of geospatial AI. Where no LiDAR survey exists, RasterFlow's Meta CHM v1 model estimates canopy height in meters from 60 cm aerial imagery. The imagery comes from the National Agriculture Imagery Program (NAIP), and the output is a raster that the same RS_ functions read.
After Wehr and Lohr's 1999 system overview, forestry led the applications. Lefsky, Cohen, Parker, and Harding showed in Lidar remote sensing for ecosystem studies (2002) that LiDAR measures canopy structure directly. Næsset's two-stage procedure (2002) introduced the area-based approach that forest management inventories use. A regression links tree measurements in field plots to LiDAR height statistics, then predicts every cell of the survey. Popescu, Wynne, and Nelson measured individual crown diameter (2003) to estimate tree volume and biomass. Khosravipour and colleagues removed a common artifact, false holes in the canopy surface, with pit-free canopy height models (2014).
Every DTM depends on separating ground returns from everything else. Axelsson's adaptive TIN (triangulated irregular network) densification (2000) builds a coarse triangle surface from the lowest local points. It adds points within angle and distance limits and rebuilds the surface each pass. Zhang, Chen, Whitman, and colleagues introduced the progressive morphological filter (2003). It removes cars, trees, and buildings with a minimum then maximum filter (an opening) whose window grows step by step. Sithole and Vosselman's experimental comparison of eight filters (2004) found that all did well on smooth rural terrain and all made errors in complex cities and on rough, vegetated slopes. Zhang, Qi, Wan, and colleagues then proposed the cloth simulation filter (2016). It turns the point cloud upside down, drops a simulated cloth on it, and keeps the points the cloth touches. It needs few parameters and ships as a plugin in the open-source CloudCompare.
Point clouds have no grid, so early networks voxelized or rendered them first. Qi, Su, Mo, and Guibas removed that step with PointNet (CVPR 2017), a network that reads raw points in any order. PointNet++ (NeurIPS 2017) added groups of nearby points at several scales. Thomas and colleagues' Kernel Point Convolution, KPConv (ICCV 2019), placed convolution kernels at points in space. Semantic3D (Hackel and colleagues, 2017) released over four billion hand-labeled terrestrial LiDAR points in eight classes. Guo and colleagues' survey of deep learning for 3D point clouds (2021) reviews the work that followed.
Airborne surveys are dense and local; spaceborne LiDAR is global and sparse. ICESat-2's land and vegetation product reports ground and canopy heights in 100 m segments along track. Neuenschwander and Pitts describe it in The ATL08 land and vegetation product for the ICESat-2 Mission (2019). GEDI samples canopy height in 25 m footprints. Potapov and colleagues extrapolated GEDI with Landsat to a 30 m global forest height map for 2019 (2021). Lang, Jetz, Schindler, and Wegner trained a deep network on GEDI and Sentinel-2 for a high-resolution canopy height model of the Earth at 10 m (2023). Tolan and colleagues reached 1 m with canopy height maps from Maxar RGB imagery (2024), imagery in red, green, and blue bands. Their vision transformer, trained on aerial LiDAR and calibrated with GEDI, had a reported mean absolute error of 2.8 m. Brandt and colleagues' CHMv2 (2026) rebuilt the model on the DINOv3 image network. They validated it against independent airborne LiDAR and against GEDI and ICESat-2 observations. It is the canopy layer in the Portland queries.
A 1,000 km² county at the QL1 minimum of 8 pulses per square meter is 8 billion pulses. The 3DEP archive on AWS held more than 12 trillion points when USGS announced it. Stoker and Miller tested whether data from many producers are accurate and consistent across 3DEP (2022). At that scale, analysis needs many machines working in parallel. Jia Yu, Jinxuan Wu, and Mohamed Sarwat introduced spatial processing on Apache Spark with GeoSpark (ACM SIGSPATIAL 2015). Yu, Zhang, and Sarwat described its partitioning, indexing, and join design in Spatial data management in Apache Spark: the GeoSpark perspective and beyond (GeoInformatica, 2019). GeoSpark became Apache Sedona. Its raster functions run the zonal statistics in the queries above in parallel. Raster processing at scale: the out-of-database architecture behind WherobotsDB explains why a zonal statistic reads only the Cloud Optimized GeoTIFF (COG) tiles under each zone.
Load LiDAR-derived elevation and canopy rasters into WherobotsDB with a Wherobots free trial at cloud.wherobots.com.
LiDAR data is a point cloud: millions of 3D points measured by a laser scanner that times the round trip of each light pulse. Each point has x, y, and z coordinates, an intensity, a return number, and a classification such as ground, vegetation, building, or water. It is the source of most high-resolution terrain, canopy height, and building height maps.
A scanner on an aircraft, drone, vehicle, tripod, or satellite emits laser pulses. It records the time each reflection takes to return, and the speed of light turns that time into a distance. A satellite positioning (GNSS) receiver records where the scanner was, and an inertial unit records which way it pointed. Together they place each return as a point with map coordinates.
A point cloud is any set of 3D points, and LiDAR data is a point cloud produced by a laser scanner. Photogrammetry also produces point clouds, from overlapping photographs. LiDAR data adds per-point intensity, return number, and multiple returns per pulse, which is how it maps the ground under vegetation.
LiDAR data is stored in LAS, the binary standard of the American Society for Photogrammetry and Remote Sensing (ASPRS). LAZ is LAS with lossless compression; the LASzip tool writes it at 7 to 25 percent of the LAS size. Cloud Optimized Point Cloud (COPC) is a LAZ 1.4 file with its points grouped into nested cubes, called an octree, so a client can download one area over HTTP. Terrain and canopy grids made from the points are written as GeoTIFF. Points can be converted to GeoParquet tables for SQL analysis.
USGS quality levels QL1 and QL2 both require 10 cm vertical accuracy on open ground, measured as the root mean square error of heights (RMSEz). QL1 needs at least 8 pulses per square meter and QL2 at least 2. Accuracy falls under dense canopy and on steep slopes, where fewer pulses reach the ground and ground filtering makes more errors.
Free LiDAR data for the United States is on the USGS 3D Elevation Program through LidarExplorer and The National Map, and on AWS as LAZ tiles and Entwine Point Tiles. NOAA’s Data Access Viewer serves coastal LiDAR, and OpenTopography hosts academic surveys worldwide. The Netherlands (AHN) and England (Environment Agency National LIDAR Programme) publish national LiDAR data as open data.
Yes. Free sources include USGS 3DEP for the United States, NOAA’s coastal surveys, OpenTopography’s research data, and several national programs.
Open a LAS, LAZ, or COPC file in a point cloud viewer such as CloudCompare, or in a browser viewer that streams COPC. USGS LidarExplorer and OpenTopography preview coverage before download. To view terrain or canopy height, load the derived GeoTIFF DEM or canopy height model in any GIS (geographic information system).
LiDAR measures distance with light. A scanner fires laser pulses many thousands of times a second and times how long each one takes to bounce back. Each measurement becomes a point in 3D, and together the points form a model of the ground, trees, and buildings below.
A digital surface model (DSM) is gridded from the first or highest returns and includes trees and roofs. A digital terrain model (DTM) is gridded from ground returns only. A canopy height model is the DSM minus the DTM: the height of vegetation above the ground. The same subtraction inside a building footprint gives building height.
Classification codes are ASPRS standard values stored with each LAS point. The common ones are 1 unclassified, 2 ground, 3 to 5 low, medium, and high vegetation, 6 building, 7 low noise, 9 water, 17 bridge deck, and 18 high noise. Values 64 to 255 are free for user definitions.
What is a Geofence? How Geofencing Works, With Examples
A geofence is a virtual boundary that triggers an action when a device enters or leaves it. How geofencing works, its uses and limits, and a test on ship data.
What is Map Matching? GPS Traces, Methods and Examples
Map matching places each GPS record on the road a vehicle drove. How it works, nearest-road snapping vs hidden Markov models, and a test on Ann Arbor drives.
Geospatial Big Data: Why Analytics Is Hard to Scale
Teams struggle to scale geospatial analytics for five reasons: data volume, mixed formats, many coordinate systems, spatial joins and pipeline upkeep.
share this article
Awesome that you’d like to share our articles. Where would you like to share it to:
