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)