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"))