.
Processing a terrestrial laser scan of a single tree… (TLS data)
#####New libraries to install: concaveman, shapefiles
library("lidR")
library("rgdal")
library("concaveman")
library("shapefiles")
#####Load Files
mypath<- "U:/FANR5640/Spring2019/labs/lab09Data"
myfile<- "SingleTree.laz"
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<- lastmp
plot(lasd)
summary(lasd)
#####use for loop later on
valMinZ<- lasd@header@PHB$'Min Z'
valMaxZ<- lasd@header@PHB$'Max Z'
valDifZ<- valMaxZ - valMinZ
print(valMinZ)
print(valMaxZ)
#####manually set ground to valMinZ -
#####the next line creates a 'ground' object whose raster values are not actually the ground
grnd<- grid_metrics(lasd,min(Z),0.05)
#####reset all cell values tothe minimum valMinZ value
grnd[grnd > valMinZ]<- valMinZ
#####Normalize lasd based on the grnd object created above
lasd_n<- lasnormalize(lasd,grnd)
summary(lasd_n)
plot(lasd_n)
#####maximum heigt
valMaxZNorm<- lasd_n@header@PHB$'Max Z'
col<- pastel.colors(250)
#####chm = surface model
chm<- grid_canopy(lasd_n, res=0.5,p2r(0.3))
ker <- matrix(1,3,3)
#####smooth the canopy height mocel
chm <- raster::focal(chm, w = ker, fun = mean, na.rm = TRUE)
plot(chm)
####tree detection required for convex_hull procedures below
#tree_detection(lasd,lmf(ws=25,hmin=(65)))
ttops <- tree_detection(chm, lmf(4, 2))
lasd_n<- lastrees(lasd_n,mcwatershed(chm,ttops))
plot(lasd_n,color="treeID")
#####SLICE TREE
seqStep<- 3
for (i in seq(0,valMaxZNorm,seqStep)){
charI<- as.character(i)
print(charI)
charI<- gsub("\\.","_",charI)
lowval<- i
highval<- i + seqStep
#####the point cloud below stores the sliced point cloud and is overwritten each loop
lasd_n_xxft<- lasfilter(lasd_n,((Z>=lowval) & (Z<highval)))
#####https://rdrr.io/cran/lidR/man/tree_hulls.html
convex_hulls<- tree_hulls(lasd_n_xxft)
sp::plot(convex_hulls)
bbox_hulls<- tree_hulls(lasd_n_xxft,"bbox")
sp::plot(bbox_hulls)
concave_hulls<- tree_hulls(lasd_n_xxft,"concave")
sp::plot(concave_hulls)
#####Export hulls to shapefile
#####Populate new shapefile with area, height range, bottom height, top height, volume
#######
data.frame(convex_hulls)
convex_hulls$myarea<- (convex_hulls@polygons[[1]]@area)
convex_hulls$htrange<- paste(as.character(i) , " - " , as.character(i + seqStep))
convex_hulls$bottomheight<- lowval
convex_hulls$topheight<- highval
convex_hulls$volume<- convex_hulls$myarea * seqStep
#####pattern to create output file name
outshape<-(paste("convexHull_",charI))
writeOGR(obj=convex_hulls,dsn=mypath, layer=outshape, driver="ESRI Shapefile", overwrite=TRUE)
print("written convex hull")
data.frame(bbox_hulls)
bbox_hulls$myarea<- (bbox_hulls@polygons[[1]]@area)
bbox_hulls$htrange<- paste(as.character(i) , " - " , as.character(i + seqStep))
bbox_hulls$topheight<- highval
bbox_hulls$bottomheight<- lowval
bbox_hulls$volume<- bbox_hulls$myarea * seqStep
#####pattern to create output file name
outshape<-(paste("boundingBox_",charI))
writeOGR(obj=bbox_hulls,dsn=mypath, layer=outshape, driver="ESRI Shapefile", overwrite=TRUE)
print("written bounding box")
data.frame(concave_hulls)
concave_hulls$myarea<- (concave_hulls@polygons[[1]]@area)
concave_hulls$htrange<- paste(as.character(i) , " - " , as.character(i + seqStep))
concave_hulls$topheight<- highval
concave_hulls$bottomheight<- lowval
concave_hulls$volume<- concave_hulls$myarea * seqStep
#####pattern to create output file name
outshape<-(paste("concaveHull",charI))
writeOGR(obj=concave_hulls,dsn=mypath, layer=outshape, driver="ESRI Shapefile", overwrite=TRUE)
print("written concave hull")
}
.