Monday, you stepped through a simple analysis in ArcMap. You created a normalized surface model (I called it the DZM in the lab document) by differencing the surface model (DSM) and the terrain model (DTM). The resulting raster layer contains cells whose values are 0 or near-zero. These cells represent the bare ground. The remaining values reflect an object’s height above the bare ground. I want to demonstrate an alternate workflow where we process the point cloud directly using R. I
If you want to try to follow along, here is the data I am using (data download link).
###lidR help: https://www.rdocumentation.org/packages/lidR/versions/2.0.0
library("lidR")
library("rgdal")
library("EBImage")
#####must change data paths
mypath<- "U:/FANR5640/Spring2019/labs/lab06Data/SEasternGrowers"
shapepath<- "U:/FANR5640/Spring2019/labs/lab06Data/SEasternGrowers"
shapename<- "Mask"
shapef<- file.path(shapepath,shapename)
shape_spdf<- readOGR(dsn = shapepath, layer = shapename)
shape_p<- shape_spdf@polygons[[1]]@Polygons[[1]]
myfile<- "clpSoutheasternGrowers_PointCloud.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)
las = lasclip(lastmp, shape_p)
las_original = las
plot(las)
summary(las)
#####Ground Classification
plot(las, color = "Classification")
las<- las_original
las<- lasground(las,csf(cloth_resolution = 0.5, class_threshold = 0.25), last_returns = FALSE)
las_ground <- las
plot(las, color = "Classification")
#####OR
##las<- las_original
##ws<- seq(0.75,0.75,0.75)
##th<- seq(0.5,1, length.out=length(ws))
##las<- lasground(las,pmf(ws,th), last_returns = FALSE)
#####NORMALIZE
las <- lasnormalize(las, knnidw(k = 8, p = 2))
las_normal<- las
summary(las)
plot(las_original)
plot(las_normal)
plot(las_normal, color = "Classification")
summary(las_Norm)
summary(las_original)
#####Canopy Height model
chm <- grid_canopy(las, res = 0.5, pitfree(c(0,2,5,10,15), c(0, 1.5)))
chm <- grid_canopy(las, res = 0.15, p2r(0.05))
plot(chm)
writeRaster(chm,filename = file.path(mypath,"chm4a.tif"), format="GTiff",overwrite=TRUE)
##writeLAS(las_15Norm,file.path(mypath,"i2_06Jan19_las_15norm.las"))
##writeLAS(las_1Norm,file.path(mypath,"i2_06Jan19_las_1norm.las"))
##writeLAS(las_3Norm,file.path(mypath,"i2_06Jan19_las_3norm.las"))
##writeLAS(las_original,file.path(mypath,"i2_06Jan19_las_original.las"))
##f = function(x) { x * 0.07 + 3}
ttops = tree_detection(las_original, lmf(ws=1, hmin=0.5, shape = "circular"))
x=plot(las_original)
add_treetops3d(x,ttops)
###or
plot(ttops)
writeOGR(obj=ttops,dsn=shapepath, layer="ttops2a", driver="ESRI Shapefile", overwrite=TRUE)