LiDAR in R
Sabrina Desjardins, PhD
2024-05-22
Using lidR to process point clouds
By using the lidR package, we can read and validate point clouds, add classifications (i.e., ground vs vegetation), generate a digital terrain model (DTM), generate a digital surface model (DSM), normalize heights, and produce single-tree detection, among other applications. Below, we will look into reading and validating point cloud data, generate a DTM, look into normalization, and produce single-tree detection. All of the examples provided in this section are from The lidR package book
Reading and Validating Point Clouds
The function readLAS() reads and creates an object that contains both the header and the payload. To save memory, readLAS() can take an optional parameter select which enables the user to selectively load the attributes of interest. Meanwhile, the function las_check() checks if a LAS object meets the ASPRS LAS specifications and whether it is valid for processing, giving warnings if otherwise
las <- readLAS("FileName.las")
las <- readLAS("FileName.las", select = "xyz") # load XYZ only
las <- readLAS("FileName.las", select = "xyzi") # load XYZ and intensity only
las <- readLAS("FileName.las", filter = "-keep_first") # Read only first return
las_check(las)Classifying ground data
The lidR package currently provides three algorithms for classifying
ground points using the function classify_ground():
1. Progressive morphological filter (PMF) described by Zhang et
al. 2003
2. Cloth simulation function (CSF) described by Zhang et al. 2016
3. Multi-scale curvature classification (MCC) described by Evans and
Hudak 2016
las <- classify_ground(las, algorithm = pmf(ws = 5, th = 3)) #PMF
las <- classify_ground(las, algorithm = csf()) #CSF
las <- classify_ground(las, mcc(1.5,0.3)) #MCCGenerating DTM
The 3 main ways of generating a DTM are:
- Triangular irregular network (TIN): this is the most simple solution because it involves no parameters. The drawbacks of the method are that it creates a non-smooth DTM and that it cannot extrapolate the terrain outside the convex hull delimited by the ground points since there are no triangle facets outside the convex hull. Moreover, the interpolation is weak at the edges because large irrelevant triangles are often created.
dtm_tin <- rasterize_terrain(las, res = 1, algorithm = tin())- Invert distance weighing (IDW): one of the simplest and most readily available methods that can be applied to create DTMs. It is based on an assumption that the value at an unsampled point can be approximated as a weighted average of values at points within a certain cut-off distance d, or from a given number k of closest neighbours. Compared to tin() this method is more robust to edge artifacts because it uses a more relevant neighbourhood but generates terrains that are “bumpy” and probably not as realistic as those generated using TINs.
dtm_idw <- rasterize_terrain(las, algorithm = knnidw(k = 10L, p = 2))- Krigging: most advanced approach and utilizes advanced geostatistical interpolation methods that take into account the relationships between the returns and their respective distances from each other. This method is very advanced, difficult to manipulate, and extremely slow to compute, but probably provides the best results with minimal edge artifacts.
dtm_kriging <- rasterize_terrain(las, algorithm = kriging(k = 40))A hillshade layer adds shadows to the DTM and is done using functions from the terra package. The terrain() and hillShade() functions can be combined to take the DTM raster layers as input and return a hillshade raster:
dtm <- rasterize_terrain(las, algorithm = tin(), pkg ="terra")
dtm_prod <- terrain(dtm, v = c("slope", "aspect"), unit = "radians")
dtm_hillshade <- shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
plot(dtm_hillshade, col =gray(0:30/30), legend = FALSE)
Point Cloud Normalization
Point cloud normalization interpolates the elevation of every single point locations using ground points. To compute the continuous normalization, we can feed normalize_height() with an algorithm for spatial interpolation instead of a raster.
nlas <- normalize_height(las, knnidw())
#All the ground points should be exactly 0. Lets check it:
hist(filter_ground(nlas)$Z, breaks = seq(-0.6, 0.6, 0.01), main = "", xlab = "Elevation")
*Bonus: to reverse normalization, we can use the function unnormalize_height()
Individual Tree Detection
Individual tree detection (ITD) is the process of spatially locating trees and extracting height information. Individual tree segmentation (ITS) is the process of individually delineating detected trees. In lidR, detecting and segmenting functions are decoupled to maximize flexibility. Tree tops are first detected using the locate_trees() function, followed crown delineation using segment_trees().
During ITD, tree tops can be detected by applying a Local Maximum Filter (LMF) on the loaded data set. The LMF can be applied with either a constant or variable window size. The following example is a window size of 5 meters (meaning that for a given point the algorithm looks to the neigbourhood points within a 2.5 radius circle to figure out if the point is the local highest). While the algorithm does not need any CHM to work, one was provided for better visualization
ttops <- locate_trees(las, lmf(ws = 5))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)
A large window is suitable for large scattered trees while a small windows size if preferable for small and close trees. In reality trees of variable sizes may be present in a single scene leading to sub-optimal outputs. To overcome this issue, a windows size that adapts to the height of pixels or height of the points in our case becomes necessary.
f <- function(x) {x * 0.1 + 3}
heights <- seq(0,30,5)
ws <- f(heights)
#When applied within lmf(), the function yields the following result:
ttops <- locate_trees(las, lmf(f))
plot(chm, col = height.colors(50))
plot(sf::st_geometry(ttops), add = TRUE, pch = 3)
Geospatial Raster and Vector Data
Raster and vector data in R are typically manipulated using the terra and sf packages, respectively. The following examples provided in this section are from the “Introduction to Geospatial Raster and Vector Data with R” tutorial by datacarpentry.org
DSM data
Before reading in the data, it is best to check the metadata to better understand what you are working with (and to reduce computational load if the data is not what you need)
describe(DSMfile)## [1] "Driver: GTiff/GeoTIFF"
## [2] "Files: NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif"
## [3] "Size is 1697, 1367"
## [4] "Coordinate System is:"
## [5] "PROJCRS[\"WGS 84 / UTM zone 18N\","
## [6] " BASEGEOGCRS[\"WGS 84\","
## [7] " ENSEMBLE[\"World Geodetic System 1984 ensemble\","
## [8] " MEMBER[\"World Geodetic System 1984 (Transit)\"],"
## [9] " MEMBER[\"World Geodetic System 1984 (G730)\"],"
## [10] " MEMBER[\"World Geodetic System 1984 (G873)\"],"
## [11] " MEMBER[\"World Geodetic System 1984 (G1150)\"],"
## [12] " MEMBER[\"World Geodetic System 1984 (G1674)\"],"
## [13] " MEMBER[\"World Geodetic System 1984 (G1762)\"],"
## [14] " MEMBER[\"World Geodetic System 1984 (G2139)\"],"
## [15] " ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
## [16] " LENGTHUNIT[\"metre\",1]],"
## [17] " ENSEMBLEACCURACY[2.0]],"
## [18] " PRIMEM[\"Greenwich\",0,"
## [19] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [20] " ID[\"EPSG\",4326]],"
## [21] " CONVERSION[\"UTM zone 18N\","
## [22] " METHOD[\"Transverse Mercator\","
## [23] " ID[\"EPSG\",9807]],"
## [24] " PARAMETER[\"Latitude of natural origin\",0,"
## [25] " ANGLEUNIT[\"degree\",0.0174532925199433],"
## [26] " ID[\"EPSG\",8801]],"
## [27] " PARAMETER[\"Longitude of natural origin\",-75,"
## [28] " ANGLEUNIT[\"degree\",0.0174532925199433],"
## [29] " ID[\"EPSG\",8802]],"
## [30] " PARAMETER[\"Scale factor at natural origin\",0.9996,"
## [31] " SCALEUNIT[\"unity\",1],"
## [32] " ID[\"EPSG\",8805]],"
## [33] " PARAMETER[\"False easting\",500000,"
## [34] " LENGTHUNIT[\"metre\",1],"
## [35] " ID[\"EPSG\",8806]],"
## [36] " PARAMETER[\"False northing\",0,"
## [37] " LENGTHUNIT[\"metre\",1],"
## [38] " ID[\"EPSG\",8807]]],"
## [39] " CS[Cartesian,2],"
## [40] " AXIS[\"(E)\",east,"
## [41] " ORDER[1],"
## [42] " LENGTHUNIT[\"metre\",1]],"
## [43] " AXIS[\"(N)\",north,"
## [44] " ORDER[2],"
## [45] " LENGTHUNIT[\"metre\",1]],"
## [46] " USAGE["
## [47] " SCOPE[\"Engineering survey, topographic mapping.\"],"
## [48] " AREA[\"Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.\"],"
## [49] " BBOX[0,-78,84,-72]],"
## [50] " ID[\"EPSG\",32618]]"
## [51] "Data axis to CRS axis mapping: 1,2"
## [52] "Origin = (731453.000000000000000,4713838.000000000000000)"
## [53] "Pixel Size = (1.000000000000000,-1.000000000000000)"
## [54] "Metadata:"
## [55] " AREA_OR_POINT=Area"
## [56] "Image Structure Metadata:"
## [57] " COMPRESSION=LZW"
## [58] " INTERLEAVE=BAND"
## [59] "Corner Coordinates:"
## [60] "Upper Left ( 731453.000, 4713838.000) ( 72d10'52.71\"W, 42d32'32.18\"N)"
## [61] "Lower Left ( 731453.000, 4712471.000) ( 72d10'54.71\"W, 42d31'47.92\"N)"
## [62] "Upper Right ( 733150.000, 4713838.000) ( 72d 9'38.40\"W, 42d32'30.35\"N)"
## [63] "Lower Right ( 733150.000, 4712471.000) ( 72d 9'40.41\"W, 42d31'46.08\"N)"
## [64] "Center ( 732301.500, 4713154.500) ( 72d10'16.56\"W, 42d32' 9.13\"N)"
## [65] "Band 1 Block=1697x1 Type=Float64, ColorInterp=Gray"
## [66] " Min=305.070 Max=416.070 "
## [67] " Minimum=305.070, Maximum=416.070, Mean=359.853, StdDev=17.832"
## [68] " NoData Value=-9999"
## [69] " Metadata:"
## [70] " STATISTICS_MAXIMUM=416.06997680664"
## [71] " STATISTICS_MEAN=359.85311802914"
## [72] " STATISTICS_MINIMUM=305.07000732422"
## [73] " STATISTICS_STDDEV=17.83169335933"If there are no issues with the metadata, we can read in the file using rast(), and can find simplified summary statistics using summary(). Important to note, rast() can only import one single band at a time - to see the number of bands in the raster file we can use nlyr()
DSM_HARV <- rast(DSMfile) #read in the DSM raster file
summary(values(DSM_HARV)) #summary stats## HARV_dsmCrop
## Min. :305.1
## 1st Qu.:345.6
## Median :359.7
## Mean :359.9
## 3rd Qu.:374.3
## Max. :416.1nlyr(DSM_HARV) #number of layers## [1] 1To visualize the raster data in ggplot2, we need to use the terra package to convert the file into a dataframe before plotting
DSM_HARV_df <- as.data.frame(DSM_HARV, xy = TRUE) #convert to dataframe
ggplot() +
geom_raster(data = DSM_HARV_df , aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_viridis_c() +
coord_quickmap()
We can also add hillshade raster data to the plot. First we load in the hillshade data usng rast(), convert it to a dataframe, then overlay it on the DSM plot:
DSM_hill_HARV <- rast(DSMhillfile)
#convert to df so we can use it in ggplot2
DSM_hill_HARV_df <- as.data.frame(DSM_hill_HARV, xy = TRUE)
#now we can overlay the DSM data ontop of the hillshade data
ggplot() +
geom_raster(data = DSM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dsmCrop)) +
geom_raster(data = DSM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DSMhill)) +
scale_fill_viridis_c() +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
ggtitle("Elevation with hillshade") +
coord_quickmap()
DTM Data
DTM data can be manipulated much like DSM data, so in this example we look at when the metadata for the hillshade and DTM rasters are different:
#load data
DTM_HARV <- rast(DTMfile)
DTM_hill_HARV <- rast(DTMhillfile)
#check summarized metadata
terra::crs(DTM_HARV, proj = TRUE)## [1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"terra::crs(DTM_hill_HARV, proj = TRUE)## [1] "+proj=longlat +datum=WGS84 +no_defs"In the above example, DTM_HARV is in the UTM projection, with units of meters. DTM_hill_HARV is in Geographic WGS84 - which is represented by latitude and longitude values. So lets reproject DTM_hill_HARV to UTM:
DTM_hill_UTMZ18N_HARV <- project(DTM_hill_HARV,
crs(DTM_HARV))
#double check to make sure it worked
terra::crs(DTM_hill_UTMZ18N_HARV, proj = TRUE)## [1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"Now that we reprojected the hillshade data we need to check the resolution of it vs the DTM data
res(DTM_hill_UTMZ18N_HARV)## [1] 1.001061 1.001061res(DTM_HARV)## [1] 1 1They don’t match - so we need to add a resolution line when coding the reprojection:
DTM_hill_UTMZ18N_HARV <- project(DTM_hill_HARV,
crs(DTM_HARV),
res = res(DTM_HARV))
#double check:
res(DTM_hill_UTMZ18N_HARV)## [1] 1 1Now the reprojected hillshade data and DTM data need to be converted into a dataframe so we can plot using ggplot2
DTM_hill_HARV_df <- as.data.frame(DTM_hill_UTMZ18N_HARV, xy = TRUE)
DTM_HARV_df <- as.data.frame(DTM_HARV, xy = TRUE)
ggplot() +
geom_raster(data = DTM_HARV_df ,
aes(x = x, y = y,
fill = HARV_dtmCrop)) +
geom_raster(data = DTM_hill_HARV_df,
aes(x = x, y = y,
alpha = HARV_DTMhill_WGS84)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
coord_quickmap()
Canopy Height Model (CHM)
CHM subtracts the DTM from the DSM to find the actual height of trees. Using the lapp() function we can subtract DTM from DSM:
CHM_ov_HARV <- lapp(sds(list(DSM_HARV, DTM_HARV)),
fun = function(r1, r2) { return( r1 - r2) })
#convert to a df so we can plot using ggplot
CHM_ov_HARV_df <- as.data.frame(CHM_ov_HARV, xy = TRUE)
#plot:
ggplot() +
geom_raster(data = CHM_ov_HARV_df,
aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_quickmap()
If we want to export as a GEOtiff we can use the writeRaster function:
writeRaster(CHM_ov_HARV, "CHM_HARV.tiff",
filetype="GTiff",
overwrite=TRUE,
NAflag=-9999)Adding a vector-based polygon on a CHM raster
We can use the crop() function to crop a raster to the extent of another spatial object. Plotting the CHM with the AOI boundary polygon:
## Reading layer `HarClip_UTMZ18' from data source
## `/Users/sabrinad/Desktop/r_projects/LiDAR Intro/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
## Projected CRS: WGS 84 / UTM zone 18Nggplot() +
geom_raster(data = CHM_ov_HARV_df, aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
geom_sf(data = aoi_boundary_HARV, color = "blue", fill = NA) +
coord_sf()
Now using crop(), we can create a new object with only the portion of the CHM data that falls within the boundaries of the AOI:
CHM_HARV_Cropped <- crop(x = CHM_ov_HARV, y = aoi_boundary_HARV)Now we can plot the cropped CHM data, along with a boundary box showing the full CHM extent. However, remember, since this is raster data, we need to convert to a data frame in order to plot using ggplot. To get the boundary box from CHM, the st_bbox() will extract the 4 corners of the rectangle that encompass all the features contained in this object. The st_as_sfc() converts these 4 coordinates into a polygon that we can plot
CHM_HARV_Cropped_df <- as.data.frame(CHM_HARV_Cropped, xy = TRUE)
ggplot() +
geom_sf(data = st_as_sfc(st_bbox(CHM_ov_HARV)), fill = "green",
color = "green", alpha = .2) +
geom_raster(data = CHM_HARV_Cropped_df,
aes(x = x, y = y, fill = HARV_dsmCrop)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_sf()
The plot above shows that the full CHM extent (plotted in green) is much larger than the resulting cropped raster. Our new cropped CHM now has the same extent as the aoi_boundary_HARV object that was used as a crop extent (blue border below).
ggplot() +
geom_raster(data = CHM_HARV_Cropped_df,
aes(x = x, y = y, fill = HARV_dsmCrop)) +
geom_sf(data = aoi_boundary_HARV, color = "blue", fill = NA) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_sf()
There are plenty of other data manipulation possibilities within the terra, sf, and lidR packages that were not highlighted in this document due to time constraints. A full R script with other code not highlighted in this document is available upon request.