This R script allows the user to process multiple LAS/LAZ files at once using the catalog functionality of the lidR library. If you look closely, this and the R code from Monday are almost identical. Line 12, here, we create a LAS catalog instead of reading in an individual LAS like in line 20 of Monday’s code.
#####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<- "U:/FORS7690/Fall2020/lectures/MaineProperty/lidar2015/las_utm19"
outfolder<- "out_catalog"
if (file.exists(file.path(mypath,outfolder))){
} else {
dir.create(file.path(mypath,outfolder))
}
lascat<- catalog(mypath)
spplot(lascat,"Max.Z")
###########################
#####Generate Digital Terrain Model "bare ground model"
###create subdirectory in 'mypath' called "tiles"
newdir<- "dtm"
if (file.exists(file.path(mypath,newdir))){
} else {
dir.create(file.path(mypath,newdir))
}
###opt_output_files specifies the output directory; files will be tagged with the left/bottom XY
###will use this option repetedly throughout
opt_output_files(lascat) <- file.path(mypath,newdir,"dtm_{XLEFT}_{YBOTTOM}")
###create digital terrain model, 2meter cell, using knn w/ 10 neighbors w/ a power of 2
dtm <- grid_terrain(lascat, res = 2, knnidw(k = 10, p = 2), keep_lowest = FALSE)
plot(dtm)
###saves new raster called "a_dtm.tif" in mypath
###after you load raster in ArcMap > Calculate Statistics to display
writeRaster(dtm, file.path(mypath,newdir,"/a_dtm.tif"))
###create slope, aspect, and hillshade
slope <- terrain(dtm, opt='slope')
aspect <- terrain(dtm, opt='aspect')
hs <- hillShade(slope, aspect, angle=45, direction=315)
###output slope, aspect, and hillshade rasters to the mypath folder
writeRaster(slope, file.path(mypath,newdir,"/a_slope.tif"))
writeRaster(aspect, file.path(mypath,newdir,"/a_aspect.tif"))
writeRaster(hs, file.path(mypath,newdir,"/a_hillshade.tif"))
############################
#####Normalize 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
newdir<- "norm_pointcloud"
if (file.exists(file.path(mypath,newdir))){
} else {
dir.create(file.path(mypath,newdir))
}
opt_output_files(lascat) <- file.path(mypath,newdir,"ndtm_{XLEFT}_{YBOTTOM}")
lasnorm<- lasnormalize(lascat,dtm, na.rm=TRUE)
plot(lasnorm)
############################
#####Create Digital Surface Model
###create and output DSM using the approach described here:
###https://github.com/Jean-Romain/lidR/wiki/Rasterizing-perfect-canopy-height-models
newdir<- "dsm"
if (file.exists(file.path(mypath,newdir))){
} else {
dir.create(file.path(mypath,newdir))
}
opt_output_files(lascat) <- file.path(mypath,newdir,"dsm_{XLEFT}_{YBOTTOM}")
###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(lascat, res=2, pitfree(c(0,2,5,10,15), c(0, 1)))
###Add to ArcMap and calculate statistics (catalog>right-click>calculate statistics
writeRaster(dsm,file.path(mypath,newdir,"/a_dsm.tif"))
#################################
#####Generate Canopy Height Model (CHM)
###generate normalized surface model, AKA CHM... Canopy Height Model
###we are working with rasters here
ndsm<- dsm - dtm
###remove all cells less than 0
ndsm[ndsm<0]=0
ndsm
plot(ndsm)
writeRaster(ndsm,file.path(mypath,newdir,"/a_chm.tif"), overwrite=TRUE)
##################################
#####Summarize across a 10 meter cell
###Calculate mean elevation of first return
newdir<- "mean_elev"
if (file.exists(file.path(mypath,newdir))){
} else {
dir.create(file.path(mypath,newdir))
}
opt_output_files(lascat) <- file.path(mypath,newdir,"meanelev10x_{XLEFT}_{YBOTTOM}")
opt_filter(lascat) <- "-keep_first"
###https://www.rdocumentation.org/packages/lidR/versions/2.2.3/topics/grid_metrics
metrics <- grid_metrics(lascat, ~mean(Z), 10)
###remove means less than zero
metrics[metrics<0]=0
plot(metrics)
writeRaster(metrics,file.path(mypath,newdir,"/a_mean10x.tif"))