GISNOTE
GISNOTE
FORS7690 Spring2019

Processing point cloud in R

by |Published

 

 

Lab 07:  Process point cloud in R!!!

In this lab, you will use the R programming language to process a point cloud.  As in previous labs, the goal is to automatically count trees across the landscape.  These methods are relevant even if you are not interested in tree-counting.  Think about how these trees are differentiated from the non-tree areas…  The ‘tree’ points have a higher Z value than their immediate surroundings.  Insert the word ‘mound’ if you are trying to locate fire ant hills or some other animal’s den on your site.

I have posted the R code below and linked the data (HERE) (or download Lab07Data.zip from the N:\FANR5640_7640\ folder) or .  Save and uncompress the file to your workspace and then open a new script in RStudio.  The general workflow is this:

  1. Read the LAS file (the point cloud)
  2. Filter the las file
  3. Classify the ground points (give you two approaches)
  4. Generate a canopy height model (CHM)
  5. Individual tree detection

Take note of the following:

  • LAS data are vector points.  Output these as a shapefile using the writeOGR command.
  • The output of the grid_metrics and the CHM are raster data.  Output these as a TIFF using the writeRaster command.
  • Use the summary command to view the LAS headers. Pay attention to…
    • Area: area covered by point cloud
    • Points: number of points in the point cloud
    • Density: points per square meter
    • Min X Y Z: The minimum X & Y coordinate and the minimum Z value
    • Max X Y Z:  The maximum X & Y coordinate and the maximum Z value
  • I create a backup of my working data before I get to a parameterized command so I can run several times with the original data. For instance:
    • …
    • las<- readLAS(infile, select = “xyzirc”)
    • las_bak<- las
    • …
  • I supply the links to the lidR documentation inline. Use them to figure out what the command is doing and how the parameters affect the output.

These are the R/lidR commands you will use:

  • plot()
  • summary()
  • las_ground()
  • lasnormalize()
  • lasfilter()
  • grid_metrics()
  • grid_canopy()
  • writeRaster() – Make sure you provide a unique output name
  • writeOGR() – Make sure you provide a unique output name

There is not a deliverable for this lab, yet…  I do, however, want you to familiarize yourself with the commands you use today.  Once you’ve run and re-run the code, you will have created several raster TIFFs and vector shapefiles.  Load them in ArcMap and compare your results with what you see in the original imagery.  Bring these outputs to class Wednesday.  Also, there are points within the code that ask you to record summary results for several datasets.  Create a table like the one below.  Bring this to class, too.

Dataset # Points Density Min Z Max Z Description
…
…

Before you begin, you probably will have to install the lidR and the raster packages.  From RStudio, look for “Install Packages” under the Tools pull-down-menu.  Type lidR and raster into the middle box and hit OK.  Take the defaults.

library("lidR")
library("rgdal")
######
######1. Read the LAS file
######Change the path below to your working directory
mypath<- "G:/UAV/clp_Sparta"
myfile<- "clp_26Oct2017.las"
infile<- file.path(mypath,myfile)
las<- readLAS(infile, select = "xyzirc")
las_orig<- las
plot(las)
summary(las)
######
######Record Summary Results las
######
######
######2. reduce the number of points you will be working with
#####Reduce the density of the point cloud from 201.33 points/m3 to 10 points/m3
#####lasfilterdecimate documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/lasfilterdecimate
las<- las_orig
######the random parameter specifies the desired points-per-square-meter 
lasd<- lasfilterdecimate(las,random(10))
summary(las)
summary(lasd)
plot(lasd)
######
######Record Summary Results for lasd
######
#####I already had written the following code using the variable 'las'. I am just setting the las variable (used later) to the decimated las file.
las<- lasd
las_orig<- las

######
######3. 
#####Classify 'ground' points, will use for terrain model (ground model) for topographic analysis
#####Link to lasground/csf documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/csf
plot(las, color = "Classification")
summary(las)
######reset las to the original (decimated) point cloud
las<- las_orig  
#####adjusting cloth_resolution, class_threshold, and rigidness parameters affects the ground point classification 
las<- lasground(las,csf(cloth_resolution = 0.5, class_threshold = 0.15, rigidness = 1), last_returns = FALSE)
las_ground <- las
plot(las_ground, color = "Classification")
summary(las_ground)
las_ground_bak<- las_ground
######change the k and p parameters to adjust the K-Nearest Neighbor search window
las_ground_normalized<- lasnormalize(las_ground,knnidw(k=20,p=2))
plot(las_ground_normalized)
summary(las_ground_normalized)
######
######Record Summary Results for las_ground_normalized, before removal of outliers
######
######Remove below 0m and above 2m points - using normalized values
las_ground_normalized<- lasfilter(las_ground_normalized,(Z >= 0))
las_ground_normalized<- lasfilter(las_ground_normalized,(Z < 2))
######
######Record Summary Results for las_ground_normalized, after removal of outliers
######
######Create Terrain model, save it out as a TIFF
######ZRef is carried over from the lasnormalize command (see help documentation)
######lasnormalize documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/lasnormalize
lasunngrd<- grid_metrics(las_ground_normalized, func=min(Zref), 2)
writeRaster(lasunngrd, filename = file.path(mypath,"treeainZ4.tif"), format="GTiff", overwrite=TRUE)
plot(las_ground_normalized,color="Z")
plot(las_ground_normalized,color="Zref")

######
######3.
#####Generate the terrain model, approach 2 
#####grid_terrain documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/grid_terrain
las<- las_orig
plot(las)
######again, adjust the k and p parameters to change the KNN search window
dtm1<-grid_terrain(las, res = 0.25, algorithm = knnidw(k=5,p = 0.5), keep_lowest = TRUE)
plot(dtm1)
writeRaster(dtm1 ,filename = file.path(mypath,"dtm5_05.tif"), format="GTiff",overwrite=TRUE)


######
######4.
#####grid_canopy documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/grid_canopy
#####CHM, canopy height model, should be near identical to the surface model exported from Agisoft
######review the grid_canopy documentation to figure out how the parameters will effect the putput
chm <- grid_canopy(las, res = 1, pitfree(c(0,2,5,10,15), c(0, 1.5)))
plot(chm)
writeRaster(chm,filename = file.path(mypath,"chm_1m.tif"), format="GTiff",overwrite=TRUE)

######
######5.
#####Individual tree detection documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/tree_detection
#####LMF algorithm documentation: https://www.rdocumentation.org/packages/lidR/versions/2.0.0/topics/lmf
######another opportunity to adjust parameters (see the documentation for assistance)
ttops = tree_detection(las_ground_normalized, lmf(ws=2, hmin=1, shape = "circular"))
######Write tree tops out as a point shapefile
writeOGR(obj=ttops,dsn=mypath, layer="ttops2", driver="ESRI Shapefile", overwrite=TRUE)

##################################################################
######EOF
##################################################################

 

######Bonus Tree Crown Segmentation

library(dplyr)
las<- lastrees(las,li2012(R=3,speed_up = 5))
col <- pastel.colors(200)
plot(las, color = "treeID", colorPalette = col)
metric<- tree_metrics(las,.stdtreemetrics)
hulls<- tree_hulls(las)
hulls@data<- dplyr::left_join(hulls@data, metric@data)
spplot(hulls, "Z")