Finding a ‘better’ crown delineation

If you don’t already have the input data, HERE it is.  The clip layer is here: CLIPLAYER.

Now that you have a good idea of how the various parameters affect the final output, I have modified the code so that we input the parameters near the beginning of the code (starting at line 30) and run the entire process.  Output images, shapefiles, and informational text files are differentiated using the myRunNumber variable.

Insert your better parameters and run the script.  Compare the tree top and crown delineation shapefiles (output on lines 73 and 84, repsectively) to the orthophoto.  Is your classification better?  How?  Can you make it more better?

#####https://rdrr.io/cran/lidR/
#####install libraries if not already
#####install.packages("lidR")
#####install.packages("raster")
#####install.packages("sp")
#####install.packages("gstat")
#####install.packages("dplyr")
library(lidR)
library(raster)
library(sp)
library(rgdal)

myfile<- "clpSG_PointCloud_th001.las"
inpath<- "c:/projects/week8"#

#####create folder for outputs
outfolder<- "outputs2"
if (file.exists(file.path(inpath,outfolder))){
} else {
  dir.create(file.path(inpath,outfolder))
}
outpath<- file.path(inpath,outfolder)
mylas<- readLAS(file.path(inpath,myfile))

# Extract a polygon from a shapefile
lakes = rgdal::readOGR(file.path(inpath), "clplayer")
lake = lakes@polygons[[1]]@Polygons[[1]]
subset_mylas = lasclip(mylas, lake)

#####PARAMETER LIST
myRunNumber<- "2"
myThreshhold<- 0.05
myResolution<- 0.5
myRigidness<- 3
myK<- 100
myPits<- c(0,.1, .2, 0.3, 0.5, 0.75, 1, 2, 5, 10)
myC<- c(0, 1, 1.5, 2)
myLMFMovingWindow<- 1
myLMFHmin<- 0.1

#####GROUND CLASSIFICATION OF POINTS; NEW OBJECT 'lasgrnd'
lasgrnd<- classify_ground(subset_mylas, csf(
  sloop_smooth = FALSE,
  class_threshold = myThreshhold,
  cloth_resolution = myResolution,
  rigidness = myRigidness
  )
)

#####OUTPUT DTM RASTER
dtm<- grid_terrain(lasgrnd,res = myResolution, algorithm = knnidw(k=myK, p=2), keep_lowest = TRUE)
outName<- paste("run",myRunNumber,"_myDTM_res",gsub("\\.","",toString(myResolution)),"_theK",myK,".tif",sep="")
writeRaster(dtm,filename=file.path(outpath,outName), overwrite=TRUE)

#####OUTPUT DSM RASTER
dsm<- grid_canopy(lasgrnd, res=myResolution, algorithm=pitfree(myPits, myC))
outName<- paste("run",myRunNumber,"myDSM_res",gsub("\\.","",toString(myResolution)),".tif",sep="")
writeRaster(dsm, filename=file.path(outpath,outName), overwrite=TRUE)

#####NORMALIZE POINT CLOUD; CREATE NEW OBJECT 'lasnorm'
lasnorm<- normalize_height(lasgrnd, knnidw(k=myK, p=2))  
#####filter points less than 0
lasnorm<- filter_poi(lasnorm, !(lasnorm$Z < 0))

#####CREATE CANOPY HEIGHT MODEL; SAME AS DSM BUT WITH NORMALIZED POINT CLOUD
chm <- grid_canopy(lasnorm, res = myResolution, pitfree(myPits,myC))
kernel<- matrix(1,3,3)
chm3a = raster::focal(chm, w=kernel, fun=mean)
chm3a = raster::focal(chm, w=kernel, fun=mean)

#####FIND TREE TOPS; EXPORT POINT SHAPEFILE
ttops<- find_trees(lasnorm,lmf(ws=myLMFMovingWindow, hmin = myLMFHmin))
myOutputShapefileName<- paste("run",myRunNumber,"_ttops",sep="")
writeOGR(ttops,dsn=outpath,layer=myOutputShapefileName, driver="ESRI Shapefile", overwrite_layer=TRUE)

#####CROWN SEGMENTATION
crowns<- segment_trees(lasnorm,dalponte2016(chm3a,ttops,th_tree=myLMFHmin))

#####DELINEATE TREE CROWNS
hulls<- delineate_crowns(crowns, func = ~list(Zmean = mean(Z)))

#####OUTPUT CROWNS; EXPORT SHAPEFILE
myOutputShapefileName<- paste("run",myRunNumber,"_crowns",sep="")
writeOGR(hulls,dsn=outpath,layer=myOutputShapefileName, driver="ESRI Shapefile", overwrite_layer=TRUE)

#####OUTPUT TEXT FILE WITH RUN # AND PARAMETERS
myOutputInformaitonfileName<- paste("run",myRunNumber,"_information",sep="")
myParameters<- NULL
myParameters<- data.frame(myRunNumber,
                    myThreshhold,
                    myResolution,
                    myRigidness,
                    myK,
                    gsub("\\,","",toString(myPits)) ,
                    gsub("\\,","",toString(myC)),
                    myLMFMovingWindow,
                    myLMFHmin
                    )
names(myParameters)<- c("myRunNumber","myThreshhold","myResolution","myRigidness","myK", "myPits","myC","myLMFMovingWindow","myLMFHmin")
write.csv(myParameters,file.path(outpath,paste(myOutputInformaitonfileName,".txt",sep="")))

Check out your tree crown polygon attribute table.  It is already populated with the tree apex height, ZTop, and  mean height, Zmean, from line 85.  Where are the blunders occurring?  Can you use this information for one final refinement?

Now that you have generated a better crown delineation.  I want you to do the following:

• Use screenshots and text to demonstrate and describe the effects of the various parameters on the tree crown delineation.  Present output from a run using a bad parameter (or set of parameters) and then using your better parameters.
• Demonstrate for:
  • grid_terrain (line 51)
  • normalize_height (line 61)
  • grid_canopy (line 66)
  • find_trees (line 72)
  • segment_trees (line 77)
• A Power Point slide with screenshots and text are fine.