Using LAZ point cloud linked on our eLC site
#####Processing 1 LAZ file. Will generate slope, aspect, and hillshade images; normalize the point cloud; create
#####canopy height model and a custom metric
#####Save the
#####might need to install libraries
#####install.packages(lidR)
library(lidR)
library(rgdal)
library(future)
plan(multisession, workers = 4L)
#####Setting my working directory, output, and my LAZ file name
mypath<- "E:/landowner"
outfolder<- "out"
if (file.exists(file.path(mypath,outfolder))){
} else {
dir.create(file.path(mypath,outfolder))
}
myfile<- "20150507_19tdk417957_utm19.las"
#####read the LAZ into memory, store as "mylas"
mylas<- readLAS(file.path(mypath,myfile))
#####Plot the point cloud; this might take some time so plot at your own risk
plot(mylas)
#####Report some high-level information about the point cloud
summary(mylas)
hist(mylas$Z)
quantile(mylas$Z, 0.99)
quantile(mylas$Z, 0.01)
##remove bottom and top 1% of points WRT elevation
mylas<- lasfilter(mylas,(Z<quantile(mylas$Z, 0.99)) & (quantile(mylas$Z, 0.01)))
hist(mylas$Z)
#####terrain model
###create digital terrain model, 2meter cell, using knn w/ 10 neighbors w/ a power of 2
###https://rdrr.io/cran/lidR/man/grid_terrain.html
dtm <- grid_terrain(mylas, res = 2, knnidw(k = 10, p = 2), keep_lowest = FALSE)
plot(dtm)
###saves new raster called "dtm.tif" in mypath
###after you load raster in ArcMap > Calculate Statistics to display
writeRaster(dtm, file.path(mypath,outfolder,"/dtm.tif"), overwrite=TRUE)
###create slope, aspect, and hillshade
slope <- terrain(dtm, opt='slope')
aspect <- terrain(dtm, opt='aspect')
hs <- hillShade(slope, aspect, angle=45, direction=315)
plot(hs)
###output slope, aspect, and hillshade rasters to the mypath folder
writeRaster(slope, file.path(mypath,outfolder,"/slope.tif"))
writeRaster(aspect, file.path(mypath,outfolder,"/aspect.tif"))
writeRaster(hs, file.path(mypath,outfolder,"/hillshade.tif"))
### Normalize the point cloud
###normalize point cloud by subtracting dtm elevation from original elevation; output is point cloud
###reset output directory for normalized point cloud
###Z values in the lasnorm catalog represent height above ground; ground located using the grid_terrain function
###here, we are working with the point cloud; exporting normalized point cloud
lasnorm<- normalize_height(mylas,dtm, na.rm = TRUE)
##uncomment below to plot point cloud
#plot(lasnorm)
summary(lasnorm)
hist(lasnorm$Z)
min(lasnorm$Z)
quantile(lasnorm$Z, 0.99)
quantile(lasnorm$Z, 0.01)
##remove <0 and very large normalized points) determined by trial and error
lasnorm<- lasfilter(lasnorm,Z>=0)
lasnorm<- lasfilter(lasnorm,Z< 28)
hist(lasnorm$Z)
min(lasnorm$Z)
### Create surface model (DSM) and another surface model based on the normalized point cloud
###create and output DSM using the approach described here:
###https://github.com/Jean-Romain/lidR/wiki/Rasterizing-perfect-canopy-height-models
###generate pit-free chm
###https://www.rdocumentation.org/packages/lidR/versions/2.2.3/topics/grid_canopy
###https://www.rdocumentation.org/packages/lidR/versions/2.2.3/topics/pitfree
###https://www.ingentaconnect.com/content/asprs/pers/2014/00000080/00000009/art00003?crawler=true
###in general pitfree(c(vertical thresholds), c(horizontal thresholds)
dsm<- grid_canopy(mylas, res=2, pitfree(c(0,2,5,10,15), c(0, 1)))
summary(dsm)
plot(dsm)
###Add to ArcMap and calculate statistics (catalog>right-click>calculate statistics
writeRaster(dsm,file.path(mypath,outfolder,"/dsm.tif"))
###canopy model based on normalized point cloud. values = tree height
dsmn<- grid_canopy(lasnorm, res=2, pitfree(c(0,2,5,10,15), c(0, 1)))
summary(dsmn)
plot(dsmn)
###Add to ArcMap and calculate statistics (catalog>right-click>calculate statistics
writeRaster(dsmn,file.path(mypath,outfolder,"/dsmn.tif"))
### Calculate the mean height, not-overlapping 10x10 ft cell size
###https://rdrr.io/cran/lidR/man/grid_metrics.html
###mean normalized chm
metrics<- grid_metrics(lasnorm,~mean(Z),10)
plot(metrics)
writeRaster(metrics,file.path(mypath,outfolder,"/chm_10ft.tif"),overwrite=TRUE)