Point Cloud Data Analysis

Center for Geospatial Analytics at North Carolina State University

Corey White adapted from Justyna Jeziorska and Helena Mitasova

Objectives

• Compare lidar and UAS SfM point cloud characteristics
• Analyze point distribution by binning
• Extract bare ground and thin a point cloud
• Computing a DEM (DSM and DTM)
• Analyze vertical structure with voxel models and transects

• 4A compared surfaces; today we go one step back, to the points the surfaces are made from
• Two clouds of the same site, Mid Pines at Lake Wheeler: 2013 airborne lidar and a 2015 UAS SfM flight
• The lab (4B) bins, classifies, interpolates, and validates those two clouds in GRASS
• Every concept today has a GRASS tool next to it; the lab uses them in the same order

Two kinds of point clouds

• Section 1: what a point cloud is, and how a lidar cloud differs from an SfM cloud

What is a point cloud

Point cloud: a dense set of points with x, y, z coordinates in 3D space, measured directly by lidar or derived from overlapping images by structure from motion

UAS SfM point cloud of Mid Pines on the Hunt Library visualization wall.

• Same data structure from two very different sensors: lidar measures, SfM reconstructs (2A)
• Everything downstream, DEMs, DSMs, canopy models, starts from these points
• The picture is the 2015 Mid Pines SfM cloud the lab uses

Lidar acquisition

• The measured time of flight of each pulse is converted to a distance
• Georeferencing comes from the platform’s position (GNSS) and orientation (IMU), fused in the INS

• Active sensor: it emits the pulse and times the echo, so it works regardless of texture and light
• Position from GNSS plus attitude from the IMU place every point; the same two sensors a UAS uses for geotags (1A), with a far better IMU
• INS is the IMU plus the filter that integrates it with GNSS; spec sheets usually just say GNSS/IMU
• The scan lines in the 3D view are the scanner sweeping across the flight line

Multiple returns

A lidar pulse can penetrate the tree canopy and produce several returns: yellow first returns, dark brown second returns.

• One pulse, several echoes: the top of the crown, branches, and often the ground
• This is what lets lidar map the ground under trees; SfM cannot do it
• The lab filters on return number to build a canopy surface and a ground surface from one file

Returns in profile

• Most pulses give one return; a quarter give two, a few percent three or more
• Last return has the highest chance of being ground, but it is not ground: a dense crown stops the pulse
• Hence the classified ground class in modern data, and the caution in the lab about last returns

Lidar points

Lidar point: [x, y, z, r, i, c, …], a georeferenced position with the return number r, the intensity i of the echo, and a class c assigned by the provider

• Additional attributes depend on the sensor and the provider: RGB from a camera, scan angle and direction, GPS time, number of returns per pulse
• Classes follow the ASPRS LAS standard: 1 unclassified, 2 ground, 7 noise, 9 water, and so on

• Intensity is the strength of the echo: roads dark, roofs bright, useful for classification
• Class is what a previous processing step decided; it can be wrong
• The lab’s 2013 file uses classes 1, 2, 7, and 11; pdal info prints the counts

Lidar preview: swath overlap

• Parallel flight lines overlap at their edges; the overlap band has twice the density
• Density is not uniform, which matters when you pick a raster resolution
• The lab’s density maps show these bands as stripes

Lidar preview: canopy and walls

• Points distributed throughout the canopy
• No points on the walls of the building

• Nadir-looking scanner: roofs and ground, but a wall is parallel to the beam and gets nothing
• Points inside the crowns come from the multiple returns
• Keep both in mind for the SfM comparison next

SfM points

SfM point: [x, y, z, R, G, B], a position reconstructed from overlapping imagery with the color sampled from the photos

• No returns, no intensity: the camera sees one surface, the first one in view
• Classes only if a later step assigns them (Metashape and WebODM can classify ground points)
• Density depends on the GSD, the overlap, and the quality setting used in dense matching

• Same x, y, z, but the rest of the record is different: color instead of echo physics
• Nothing under the canopy; the leaves are the surface the camera sees
• The dense cloud quality setting from 2B decides how many points you get

UAS SfM preview

• Only the top of the tree canopy is captured
• The building is densely sampled, including the walls

• The mirror image of the lidar slide: walls yes, canopy interior no
• Oblique views from the overlapping photos see the walls
• Trees are solid shells; there is no ground under them in this cloud

UAS SfM preview in color

Much higher density of points, with RGB from the photos.

• Several times the lidar density, and every point carries a color
• Looks like an orthophoto in 3D; that is also its weakness, it only knows what the photos show
• The lab bins both clouds at 1 m to compare densities side by side

Point cloud formats

Format	Notes
ASCII x, y, z	Oldest and simplest; large, slow, no attributes beyond the columns
LAS	Binary industry standard from ASPRS: header, point records with x, y, z, intensity, returns, class
LAZ	Losslessly compressed LAS; a fraction of the size, read by most software
COPC	Cloud optimized LAZ with a spatial index, streamed straight from a web server

The ASPRS LAS 1.4 specification defines the record and the standard classes; the USGS Lidar Base Specification sets what public lidar must deliver.

• LAS is what you get from providers; LAZ is what you store; COPC is what you stream
• The lab opens the COPC version of the Mid Pines lidar in a browser before touching GRASS
• Waveform data (the full echo) has its own proprietary formats

Looking at the points

• Section 2: before making a surface, find out how the points are distributed

Processing steps

• Preview, and analyze the point distribution
• Filter outliers
• Extract bare earth points
• Classify: canopy, buildings, and so on
• Decimate (thin) the cloud
• Transform the points into surfaces or 3D objects

• The order matters: you cannot pick a resolution before you know the density
• Provider data arrives with some of this done (classes); check it rather than trust it
• The lab does the first, the fourth by a simple rule, and the last two

Binning

Binning: computing a statistic of the points that fall into each cell of a grid at a chosen resolution, with no interpolation; in GRASS, r.in.pdal with method=

• Number of points per cell: a map of point density
• Range or standard deviation of z within the cell: a map of vertical variability
• Shows data gaps and where artifacts will appear
• The basis for choosing a resolution the data can support

• The cheapest thing you can do with a point cloud, and the most informative
• Count per cell is the density map; z range per cell says whether the cell holds flat ground or a tree
• Lab Part 2 is exactly this at 1 m for both clouds

Point density: county lidar

All returns and bare earth from the 2013 county lidar, and a terrestrial scan: density rises along swath overlaps and close to a terrestrial scanner.

• Airborne: the diagonal band is the swath overlap, twice the density
• Bare earth: the white holes are buildings and dense canopy, no ground points there
• Terrestrial: density falls with distance from the tripod; the same instrument gives 100 points per cell nearby and none far away

Point density: Mid Pines lidar

All returns

Ground

• All returns: trees light up red because every pulse through a crown gives several echoes
• Ground only: about one point per square meter, zero under the building and in dense crowns
• The diagonal stripes in both are the scan pattern; they will show up in any surface you bin

Point density: Mid Pines SfM

High point densities around trees and building edges.

• Ten to twenty points per square meter on the fields, several times the lidar
• Highest where the photos have the most texture: crowns and edges; the matcher finds more there
• Density is not accuracy; the SfM cloud has no ground under those crowns at all

Vertical range within a cell

Lidar

UAS SfM

Lidar reaches through the canopy, so it gives a larger z range and better information about the trees than the denser UAS cloud.

• Range of z per cell: zero on bare ground, tree height under a crown
• Lidar’s range is the whole tree, crown to ground; SfM’s is only the crown’s shape
• Denser is not better here: the lidar knows the tree height, the SfM cloud does not

Outliers

• Lidar: birds, airborne particles, reflective materials
• SfM: errors in point matching
• Filter with local z-min, z-max, or range thresholds; outliers survive even in processed data

• The spike is one point, hundreds of meters up; a bird, or a multipath echo
• Providers filter most outliers (class 7), not all
• Binning with min or max makes them obvious; interpolation would turn the spike into a cone

Ground and thinning

• Section 3: separating ground from everything else, and reducing the cloud to what you need

Bare ground extraction

Ground classification: labeling the points that lie on the terrain so that a DEM can be built from them alone; the rest is vegetation, buildings, and noise

• Multiple returns help but are not required; the surface has to be sampled by enough points
• Multiscale curvature classification (Evans and Hudak 2007), v.lidar.mcc in GRASS
• Progressive morphological filter (Zhang et al. 2003)

• Both algorithms grow a ground surface and reject points that stick up from it by more than a curvature or slope tolerance
• MCC works scale by scale, so it keeps terraces and rejects crowns; it is the addon in the lab’s optional part
• Works on SfM clouds too, with no returns: the geometry is enough on open ground

Decimation

• Thinning the cloud by subsampling, to make the data manageable
• The thinning threshold follows from the features that must be preserved
• Count based: keep every n-th point; preserves variations in density
• Grid based: keep one point per cell; removes variations in density
• Distance and geometry based: keep points where the surface changes; more expensive

• A 5 GB cloud does not need 5 GB of points to make a 1 m DEM
• Count based is fast and keeps the swath overlap pattern; grid based evens it out
• v.in.pdal and r.in.pdal thin on import with the skip and preserve options

Decimation: count based

Full cloud

Every third point

Preserves the relative point densities.

• One point in three, everywhere; the density gradient of the terrestrial scan survives
• Right choice when you want the data to keep telling you where the sensor was strong

Decimation: grid based

Full cloud

One point per cell

Removes the variations in point density.

• One point per cell: the near field loses most of its points, the far field keeps all of them
• Right choice before interpolation, when uneven density would bias the surface

From points to surfaces

• Section 4: three ways to turn points into a DEM or DSM, and what each costs

Three ways to a surface

Method	How	Strengths	Weaknesses
Binning	Per cell statistic: mean, min, max, nearest	No import, fast, enough for many uses	Noisy, no-data cells where points are missing
TIN	Triangulate the points into a mesh	Variable resolution, standard in 3D design	Flat triangles on roads, dams across valleys, few analysis tools, hard to share
Interpolation	Fit a surface to the points and sample it on a raster	Resolution above the point spacing, smooth, every GIS tool applies	Result depends on the method; large rasters

• Binning is r.in.pdal with method mean, min, or max; the lab’s Part 3
• TIN is what Metashape’s mesh is, with the same problems in valleys
• Interpolation is v.surf.rst; the lab’s Part 4 at 0.5 m from 1 m point spacing
• With dense clouds most interpolation methods agree; the choice matters at the edges and under trees

Binning: Jockey’s Ridge 1999 lidar

Binned at 1 m: many null cells

Binned at 3 m

• 1999 data at about one point per 3 m: at 1 m most cells are empty
• 3 m is the resolution the density supports; binning works, but it is coarse
• This is why the density map comes first

Interpolation: Jockey’s Ridge 1999 lidar

Interpolated at 1 m from the same points.

• Same 1999 points, now a 1 m surface with no holes: the spline fills between points
• The detail on the dune is real; the roughness on the flanks is the vegetation in the first returns
• The price is that the surface is a model between the points, not a measurement

TIN from the Mid Pines SfM cloud

Low density

High density

• Metashape’s mesh at two face counts; the 2B face count setting in action
• Low density: large triangles smear the building edges; high density: edges are sharp but every bush is a spike field
• Either way the mesh is for visualization and volumes, not for terrain analysis

Smoothed TIN

• Smoothing removes the matching noise and also the sharp edges
• Good for a 3D model to look at; it has quietly changed the volume of every pile and building

Interpolated SfM surface

High density SfM point cloud imported to GRASS and interpolated with the regularized spline with tension.

• v.surf.rst on the dense cloud: a raster DSM you can run every GRASS tool on
• The trees are columns, not cones: the spline follows the crown shell and drops to the ground at the edge
• This is the surface the 4A DSM comparisons were made from

Lidar and SfM DSMs of Mid Pines

Lidar

UAS SfM

• Lidar at 1 point per square meter gives a smooth, slightly generalized surface
• SfM at 15 points per square meter resolves the plow rows; it also carries the matching noise
• The lab interpolates the lidar at 0.5 m; Topic 5 differences the two

Topographic analysis

• Section 5: slope, aspect, and curvature from point cloud surfaces, and why it is harder than it looks

Why derivatives are hard here

• Point cloud DEMs are noisy, so slope and curvature can reflect the noise or the scan pattern rather than the terrain
• At high resolution a landform is represented by tens or hundreds of cells, and a 3 x 3 window sees only a fraction of it
• The result is noisy derivatives, or derivatives biased toward the point distribution

• 4A ran r.slope.aspect on a 0.3 m DSM and the speckle was visible; that is this problem
• A terrace 5 m wide is 15 cells at 0.3 m; a 3 x 3 window is a tenth of the feature
• Two fixes: smooth the DEM first, or compute the derivatives from the points

Topographic analysis with splines

• The regularized spline with tension has explicit equations for its partial derivatives, so slope, aspect, and curvatures are computed in the same run as the interpolation (Mitasova et al. 2005)
• The parameters come from the original points, tens to hundreds of them per location, not from a 3 x 3 window
• Tension and smoothing tune the level of detail

In GRASS: v.surf.rst input=points elevation=dem slope=slope aspect=aspect pcurv=pcurv tension=20 smooth=1

• One command, one pass over the points: the DEM and its derivatives, consistent with each other
• Tension controls how stiff the surface is; smoothing how far it may depart from the points
• The lab runs this on the ground points and on the first returns, with slope and profile curvature

Tuning the level of detail

Same points, different tension: high tension follows every point and the scan pattern, low tension shows the landforms.

• Left: tension too high, the surface is tight to the points and the curvature maps the scan stripes
• Right: tension lower, the channels and ridges appear; the blobs are tree crowns in the first returns
• There is no single right value; it depends on the feature size you are after

Curvature and slope on a 1 m DEM

Profile curvature and slope draped over the 1 m DEM.

• Curvature finds the terrace edges and the channel; slope shows the valley sides
• Both computed by the spline from the ground points, so the terraces are continuous, not stair steps
• npmin in the lab controls how many points each segment uses; too few and segment seams appear

Vertical structure

• Section 6: when a surface is not enough, keep the third dimension

Voxel models

Voxel model: a 3D raster whose cells hold a statistic of the points inside each volume element, for example the point count; in GRASS, r3.in.xyz after setting a 3D region

Vertical profile through the point cloud, and the same profile as voxels classified by a 3D fragmentation index.

• A DSM keeps one number per column; the voxel model keeps the whole column
• Point count per voxel maps where the canopy is dense and where it is open
• The fragmentation index classifies each voxel by its neighbors: interior, edge, perforated

Voxel cross sections in 3D

• Slices through the 3D raster, placed where you want to look
• The lab’s optional part builds this voxel model from the Mid Pines lidar and slices it in the 3D view

Full waveform lidar

• A discrete return lidar records a few peaks of the echo; a full waveform system records the entire echo
• Mammoth Cave National Park: classified point cloud in LAS, raw waveform, and imagery
• Vertical structure of the canopy from the waveform, not only its top

• The returns we have been using are peaks picked from this curve by the instrument
• Keeping the curve gives the vegetation profile below the crown and bathymetry in clear water
• The Mammoth Cave data set is the example for the next two slides

Mammoth Cave: canopy

• Canopy height is DSM minus DEM, both from the same cloud
• Curvature of the first return surface: every crown is a bump, which is how individual trees are found
• The cross section shows what the DEM hides: 30 m of forest on top of the terrain

Mammoth Cave: bare earth

• Ground points only, under a full canopy; this is the lidar’s whole advantage over imagery
• The curvature map picks the plateau edge and the drainage; the sinkholes of the karst show as closed rings
• Everything here was computed from the classified LAS with the tools from this lecture

Voxel models of canopy density

Point count per voxel as a percentage, sliced through two forest stands with different structure (Petras et al. 2016).

• Same voxel approach on the Mid Pines lidar: green sparse, red dense canopy
• The two stands differ in where the dense layer sits, which no DSM would show
• The paper is the lab’s reference for all of the GRASS point cloud tools

Lidar today

• Section 7: what is changing in lidar, and where to get public data

Airborne lidar today

• Manned aircraft carry the newest sensors: full waveform, multispectral, and single photon lidar
• Single photon lidar detects one photon per return, so it maps terrain and forest structure over whole counties from high altitude in a few flights (Swatantran et al. 2016)
• Against full waveform lidar it needs fewer passes for the same density, but the cloud is noisier and needs more filtering (Mandlburger et al. 2019)
• Geiger mode lidar is the other photon counting design: each detector pixel records one photon per pulse, so the cloud is built from thousands of overlapping looks (Stoker et al. 2016)

• None of these are UAS sensors; they are large-area mapping systems on manned aircraft
• Single photon (Leica SPL100) and Geiger mode (L3Harris) are both photon counting, which is what lets them fly high and fast; they differ in the detector
• Geiger mode: avalanche photodiode array, one range per pixel per pulse, noise from sunlight, near infrared so no bathymetry; USGS tested both for 3DEP and found them usable for QL2 with more filtering than linear lidar
• Swatantran et al. mapped Garrett County, Maryland, at 12 points per square meter; canopy height and terrain matched conventional lidar
• Mandlburger et al. flew both systems over Vienna in 2018: single photon precision was mostly under 10 cm, with more noise on slopes and grass
• Full waveform was the Mammoth Cave example a few slides back

UAS lidar today

• Small spinning and solid state scanners made UAS lidar routine
  • a DJI Zenmuse L2 on a Matrice 350 maps a field at several hundred points per square meter
• Discrete return, a few returns per pulse, like the 2013 cloud in the lab
• Accuracy still rests on the GNSS and IMU of the platform

• Ten years ago UAS lidar was experimental; now it is a camera-sized payload
• It sees under the canopy like airborne lidar and at SfM densities; the GNSS/IMU quality decides the accuracy, not the scanner
• An attitude error of 0.1 degree at 100 m range moves a point 17 cm; that is why UAS lidar payloads carry a better IMU than a camera needs
• Same trade as always: active sensor, no color, more expensive than a camera

Public lidar data sources

Source	What you get
USGS 3DEP	Nationwide point clouds streamed as Entwine tiles, plus DEMs through The National Map
NOAA Digital Coast	Coastal point clouds with on the fly binning and download
NC OneMap	The statewide QL2 lidar flown 2014 to 2017: points, DEMs, and derived products
OpenTopography	Research lidar (NCALM and others) with on demand DEM generation

More on lidar in GRASS at the GRASS wiki, and a full workshop at Processing lidar and UAV point clouds in GRASS.

• The 2013 Mid Pines data in the lab is the Wake County part of the NC statewide program
• 3DEP through Entwine is the same COPC streaming idea as the lab’s browser viewer
• The workshop page is the long form of today’s lecture, with every command

Key terms

Return: one echo of a lidar pulse; a pulse through a canopy gives several, the last one most likely from the ground

Class: the label a provider or an algorithm assigned to a point: ground, unclassified, noise, water

Binning: a per cell statistic of the points, with no interpolation

Interpolation: a surface fitted to the points and sampled on a raster, at a resolution the point spacing can support

DEM vs DSM: ground points give the terrain; first returns or all points give the surface, canopy and roofs included

Voxel model: a 3D raster of point statistics, keeping the vertical structure a DSM throws away

• Midterm vocabulary, and the words the lab report has to use correctly
• The grading rubric asks for conclusions in these terms: returns, classes, binning versus interpolation

Wrap-up

Assignment 4B: Processing lidar and UAS SfM point clouds: bin and compare the densities of the two Mid Pines clouds, classify by returns, interpolate a DEM and DSM with their derivatives, and validate the lidar DEM against the surveyed GCPs

• Recap: lidar measures through the canopy with returns; SfM reconstructs the first surface in color
• Look at the density before choosing a resolution; bin for a quick look, interpolate for analysis
• Derivatives from the spline, not from a 3 x 3 window, when the DEM is noisy
• The lab validates the lidar DEM with the twelve surveyed GCPs the UAS flights used

References

Evans, Jeffrey S., and Andrew T. Hudak. 2007. “A Multiscale Curvature Algorithm for Classifying Discrete Return LiDAR in Forested Environments.” IEEE Transactions on Geoscience and Remote Sensing 45 (4): 1029–38. https://doi.org/10.1109/TGRS.2006.890412.

Mandlburger, Gottfried, Harald Lehner, and Norbert Pfeifer. 2019. “A Comparison of Single Photon and Full Waveform Lidar.” ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences IV-2/W5: 397–404. https://doi.org/10.5194/isprs-annals-IV-2-W5-397-2019.

Mitasova, Helena, Lubos Mitas, and Russell S. Harmon. 2005. “Simultaneous Spline Approximation and Topographic Analysis for Lidar Elevation Data in Open-Source GIS.” IEEE Geoscience and Remote Sensing Letters 2 (4): 375–79. https://doi.org/10.1109/LGRS.2005.848533.

Petras, Vaclav, Anna Petrasova, Justyna Jeziorska, and Helena Mitasova. 2016. “Processing UAV and Lidar Point Clouds in GRASS GIS.” The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences XLI-B7: 945–52. https://doi.org/10.5194/isprs-archives-XLI-B7-945-2016.

Stoker, Jason, Qassim Abdullah, Amar Nayegandhi, and Jayna Winehouse. 2016. “Evaluation of Single Photon and Geiger Mode Lidar for the 3D Elevation Program.” Remote Sensing 8 (9): 767. https://doi.org/10.3390/rs8090767.

Swatantran, Anu, Hao Tang, Terence Barrett, Phil DeCola, and Ralph Dubayah. 2016. “Rapid, High-Resolution Forest Structure and Terrain Mapping over Large Areas Using Single Photon Lidar.” Scientific Reports 6: 28277. https://doi.org/10.1038/srep28277.

Zhang, Keqi, Shu-Ching Chen, Dean Whitman, Mei-Ling Shyu, Jianhua Yan, and Chengcui Zhang. 2003. “A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LIDAR Data.” IEEE Transactions on Geoscience and Remote Sensing 41 (4): 872–82. https://doi.org/10.1109/TGRS.2003.810682.