#####https://www.rdocumentation.org/packages/lidR/versions/2.0.1
library(“lidR”)
library(“rgdal”)
library(“EBImage”)
#####Load Files
mypath<- “/media/biomata/SSD1/Dickens_07Dec18/out”
myfile<- “Dickens_07Dec18_llp.las”
infile<- file.path(mypath,myfile)
lastmp = readLAS(infile, select = “xyzirc”)##, Intensity = F, ReturnNumber = F, NumberOfReturn = F, ScanAngle = F, EdgeOfFlightline = F, ScanDirectionFlag = F, UserData = F, PointSourceId = F, color = T)
summary(lastmp)
lasd<- lasfilterdecimate(lastmp,random(100))
plot(lasd)
summary(lasd)
#####Classify Ground Points. This is used in the normalization process
lasd<- lasground(lasd,csf(cloth_resolution = 0.5, class_threshold = 1, rigidness = 1), last_returns = FALSE)
plot(lasd,color=”Classification”)
#####normalize, might take a little while to run
lasd_n<- lasnormalize(lasd, knnidw(k=10,p=2))
plot(lasd_n,color=”Classification”)
summary(lasd_n)
lasd_n<- lasfilter(lasd_n,(Z >= 0))
plot(lasd_n)
summary(lasd_n)
lasd_nbak<- lasd_n
#####create terrain and surface model
col<- pastel.colors(250)
#####dtm1 = terrain model “bare ground model”
##dtm1<-grid_terrain(lasd_n, res = 0.25, algorithm = knnidw(k=25,p = 0.5))
##plot(dtm1)
#####chm = surface model
chm<- grid_canopy(lasd_n, res=0.5,p2r(0.3))
plot(chm)
ker <- matrix(1,3,3)
chm <- raster::focal(chm, w = ker, fun = mean, na.rm = TRUE)
plot(chm)
#####detect individual tree tops and crowns
ttops <- tree_detection(chm, lmf(4, 2))
lasd_n<- lastrees(lasd_n,mcwatershed(chm,ttops))
lasd_n_trees<- lasfilter(lasd_n,!is.na(treeID))
plot(lasd_n_trees,color=”treeID”,colorPalette=random.colors(1000),backend=”pcv”)
metric<- tree_metrics(lasd_n_trees, .stdtreemetrics)
hulls<- tree_hulls(lasd_n_trees)
plot(hulls)
hulls@data<- dplyr::left_join(hulls@data, metric@data)
#####Rasterway
##crowns<- watershed(chm)()
##plot(lasd_n, col = random.colors(length(crowns)))
##x = plot(lasd_n, color = “treeID”, colorPalette = col)
##add_treetops3d(x, ttops)
#####Export DTM, DSM, Crowns, and tree tops; view these in ArcMap. They do not have a defined coordinate system but the data are to scale; meters are the unit of measure.
dtmOut<- “TerrainModel.tif”
chmOut<- “CanopyHeightModel.tif”
crnOut<- “TreeCrowns.tif”
ttopOut<- “TreeTops”
crnPolyOut<- “TreeCrownPoly”
##writeRaster(dtm1 ,filename = file.path(mypath,dtmOut), format=”GTiff”,overwrite=TRUE)
writeRaster(chm, file.path(mypath,chmOut), format=”GTiff”,overwrite=TRUE)
##writeRaster(crowns, file.path(mypath,crnOut), format=”GTiff”,overwrite=TRUE)
writeOGR(obj=ttops,dsn=mypath, layer=ttopOut, driver=”ESRI Shapefile”, overwrite=TRUE)
writeOGR(obj=hulls,dsn=mypath,layer=crnPolyOut,driver=”ESRI Shapefile”,overwrite=TRUE)
#####Generate statistics for each tree; Trees were identified with the “lastrees(lasd_n,mcwatershed(chm,ttops))” command
myMetrics = function(z, zr)
{
metrics = list(
zmean = mean(z),
zmax = max(z),
zmin = min(z),
zrmean = mean(zr),
zrmax = max(zr),
zrmin = min(zr)
)
return(metrics)
}
treestats<- tree_metrics(lasd_n, myMetrics(Z, Zref))
metricsOut<- “TreeMetrics.csv”
write.csv(treestats, file= file.path(mypath,metricsOut))
#####Statistics for groups of trees @ 30ft (30/3.28)m, 60ft (60/3.28)m, 90ft (90/3.28)m
#####First create filtered copies of the lasd_n point cloud
lasd_n30ft<- lasfilter(lasd_n,(Z <= (30/3.28)))
lasd_n60ft<- lasfilter(lasd_n,((Z > (30/3.28)) & (Z <= (60/3.28))))
lasd_n90ft<- lasfilter(lasd_n,(Z > (90/3.28)))
plot(lasd_n30ft, color=”treeID”, colorPalette = col)
plot(lasd_n60ft, color=”treeID”, colorPalette = col)
plot(lasd_n90ft, color=”treeID”, colorPalette = col)
####output 30ft stats
treestats30ft<-tree_metrics(lasd_n30ft, myMetrics(Z, Zref))
csvoutname<- “TreeMetrics_30ft.csv”
write.csv(treestats30ft, file = file.path(mypath, csvoutname))
####output 60ft stats
treestats60ft<-tree_metrics(lasd_n60ft, myMetrics(Z, Zref))
csvoutname<- “TreeMetrics_60ft.csv”
write.csv(treestats60ft, file = file.path(mypath, csvoutname))
####output 90ft stats
treestats90ft<-tree_metrics(lasd_n90ft, myMetrics(Z, Zref))
csvoutname<- “TreeMetrics_90ft.csv”
write.csv(treestats90ft, file = file.path(mypath, csvoutname))