Using rLiDAR and FUSION to delineate individual trees through canopy height model segmentation

example3Let’s face it. Light detection and ranging (LiDAR) is awesome and it’s completely arrived. Currently, the best sources for free nationwide LiDAR datasets are the United States Interagency Elevation Inventory, USGS Center for LIDAR Information Coordination and Knowledge, and NSF’s OpenTopography. Just as there are quite a few sources for datasets, your choices are equally diverse when it comes to tools for processing LiDAR data and for the detection and delineation of individual trees from LiDAR (see Jakubowski et al 2013 for a nice review). Personally, I use a combination of FUSION, Global Mapper’s LiDAR Module, LiForest’s implementation of Li et al’s 2012 point cloud segmentation method, Swetnam and Falk’s 2014 variable area local maxima algorithm (implemented in MatLab), and the local maximum with a fixed window size algorithm implemented in rLiDAR by Carlos Alberto Silva*.

The following is a worked example using a forested LiDAR dataset, FUSION to build a CHM (called form R) and processing the resulting output using rLiDAR*. You need to install FUSION from here, and I also assume it’s install in “C:/FUSION/”. I hope you find the script useful…

# First things first, we need to set up some directories to keep the raw data 
# separate from that produced by the analysis. I basically have an input (.las) 
# and output (.las, .dtm, etc) directory in a dropbox folder.     
mainDir <- "C:/Dropbox/LiDAR_Prj"
inDir <- "inputDirectory"
outDir<- "outputDirectory"
dir.create(file.path(mainDir, outDir), showWarnings = FALSE)

# First, we process the data in FUSION

# Read in the .las file and use FUSION to produce a .las file of points that 
# approximate the ground's surface (bare-earth points). 
# http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=94&zoom=auto,70,720
system(paste(file.path("C:","Fusion", "groundfilter.exe"),
       "/gparam:0 /wparam:1 /tolerance:1 /iterations:10",
       file.path(mainDir, outDir, "Example_GroundPts.las"),
       1,
       file.path(mainDir, inDir, "Example.las"),
       sep=" "))

# Next we use ridSurfaceCreate to compute the elevation of each grid cell using the 
# average elevation of all points within the cell. Check the manual for arguments and uasge 
# http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=88&zoom=auto,70,720 
system(paste(file.path("C:","Fusion", "gridsurfacecreate.exe"),
       file.path(mainDir, outDir, "Example.dtm"),
       "1 M M 1 12 2 2",
       file.path(mainDir, outDir, "Example_GroundPts.las"),
       sep=" "))
	   
# Next we use CanopyModel to create a canopy surface model using a LIDAR point cloud. 
# By default, the algorithm used by CanopyModel assigns the elevation of the highest return within 
# each grid cell to the grid cell center.	   
#http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=32&zoom=auto,70,720
system(paste(file.path("C:","Fusion", "canopymodel.exe"),
       paste("/ground:",file.path(mainDir, outDir, "Example.dtm"), sep=""),
       file.path(mainDir, outDir, "Example.dtm"),
       "1 M M 1 12 2 2",
       file.path(mainDir, inDir, "Example.las"),
       sep=" "))

# Lastly, we use DTM2ASCII to convert the data stored in the PLANS DTM format into ASCII raster
# an file. Such files can be imported into GIS software such as ArcGIS or QGIS.
# http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=88&zoom=auto,70,720
system(paste(file.path("C:","Fusion", "dtm2ascii.exe"),
       file.path(mainDir, outDir, "Example.dtm"),
       file.path(mainDir, outDir, "Example.asc"),
       sep=" "))

# Second, we process the resulting CHM in rLiDAR

#install.packages("rLiDAR", type="source")
#install.packages("raster", dependencies = TRUE)
library(rLiDAR)
library(raster)
library(rgeos)

# Import the LiDAR-derived CHM file that we just made in the above section and plot it
chm<-raster(file.path(mainDir, outDir, "Example.asc"))
plot(chm)

Here’s the resulting canopy height model
Rplot02

# If we want, rLiDAR can smooth the CHM using a Gaussian filter
# Set the window size (ws)
ws<-3 # dimension 3x3
# Set the filter type
filter<-"Gaussian"
# filter<-"mean"
# Set the sigma value for the Gaussian filter
sigma<-0.5
sCHM<-CHMsmoothing(chm, filter, ws, sigma) 
plot(sCHM)

#########################
# Edited - Have a look at the comments section in this post to see the discussion that triggered this fix... basically, before we can call FindTreesCHM(), it needs to be
# replaced with an edited, custom version. Why? - well as of 04/06/2016 the FindTreesCHM funciton has a couple of bugs - namely one form a call to SpatialPoints() with specifying the projection
# and one in the call to colnames() (really, it's cbind(), but who's counting) 
FindTreesCHM<-function(chm, fws = 5, minht = 1.37) 
{
    if (class(chm)[1] != "RasterLayer") {
        chm <- raster(chm)
    }
    if (class(fws) != "numeric") {
        stop("The fws parameter is invalid. It is not a numeric input")
    }
    if (class(minht) != "numeric") {
        stop("The minht parameter is invalid. It is not a numeric input")
    }
    w <- matrix(c(rep(1, fws * fws)), nrow = fws, ncol = fws)
    chm[chm < minht] <- NA
    f <- function(chm) max(chm)
    rlocalmax <- focal(chm, fun = f, w = w, pad = TRUE, padValue = NA)
    setNull <- chm == rlocalmax
    XYmax <- SpatialPoints(xyFromCell(setNull, Which(setNull == 
        1, cells = TRUE)), proj4string = crs(chm))                # Edited
    htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
    treeList <- cbind(coordinates(XYmax), htExtract)              # Edited
    colnames(treeList) <- c("x", "y", "height")
    return(treeList)
}
#########################

# Setting the fixed windo size (fws)
fws<-3 # dimention 3x3
# Set the specified height above ground for the detection break
minht<-2.0
# Create the individual tree detection list and summarize it
loc<-FindTreesCHM(sCHM, fws, minht)
summary(loc)

# Set the maxcrown parameter - maximum individual tree crown radius expected
maxcrown=10.0
# Set the exclusion parameter - A single value from 0 to 1 that represents the % of pixel
# exclusion. E.g. a value of 0.5 will exclude all of the pixels for a single tree that has 
#a height value of less than 50% of the maximum height from the same tree. Default value is 0.3.
exclusion=0.1
# Compute canopy areas for the individual tree detections
canopy<-ForestCAS(sCHM, loc, maxcrown, exclusion)

# Retrieve the boundary for individual tree detection and canopy area calculation
boundaryTrees<-canopy[[1]]

# Retrieve the list of individual trees detected for canopy area calculation
canopyList<-canopy[[2]] # list of ground-projected areas of individual tree canopies
summary(canopyList)     # summary
canopyList$crad<-sqrt(canopyList$ca/pi) # Compute the corresponding crown radii

# Write the output to a CSV file. This will make bringing it into ArcGIS or QGIS by others easier
write.csv(canopyList, file.path(mainDir, outDir, "Example_Out.csv"))

# Plot the results in ggplot

# Let's convert the tree results into spatial points to be used in ggplot
library(sp)
XY<-SpatialPoints(canopyList[,1:2])    # Spatial points
XY<-data.frame(XY) # Converted to a dataframe
# We're not dealing with lat and long, so let's rename the columns
names(XY)[names(XY)=="longitude"]<-"x"
names(XY)[names(XY)=="latitude"]<-"y"

# Rasters are a little problematic, so convert the values into rows in a dataframe with 
# the corresponding information
CHMdf <- rasterToPoints(sCHM); CHMdf <- data.frame(CHMdf)
colnames(CHMdf) <- c("X","Y","Ht")
# Build the breaks for plotting
b.chm <- seq(0,50,10)

# Plotting the individual tree canopy boundary over the CHM
library(ggplot2)
ggplot(CHMdf) +
  geom_raster(data=CHMdf,aes(X,Y,fill=Ht)) +
  scale_fill_gradientn(name="Canopy Height",colours = terrain.colors(length(b.chm))[length(b.chm):1],breaks=b.chm) +
  geom_polygon(data = fortify(boundaryTrees), aes(x=long, y=lat,
                                                  group = group),colour='black', fill='transparent')+
  geom_point(data=XY, aes(x=x, y=y), color="black", shape=3, size=0.5)+
  scale_shape_discrete(name = "Tree Locations", labels=c("Tree Centroid","")) + 
  scale_alpha(range = c(0, 0.5)) +
  ggtitle("Location map showing individual trees and associated crown areas \n segmented from a LiDAR-derived canopy height model") +
  coord_equal() + theme_bw()

You can download the script from here and here’s the resulting figure.

Rplot01* – As per Carlos suggestion in the comments, here’s how one might cite rLiDAR – “Silva, C.A., Crookston, N.L., Hudak, A.T., and Vierling, L.A. 2015. rLiDAR: An R package for reading, processing
and visualizing LiDAR (Light Detection and Ranging) data, version 0.1, accessed Oct. 15 2015.”

 

Facebooktwittergoogle_plusredditpinterestlinkedinmailFacebooktwittergoogle_plusredditpinterestlinkedinmail

More Fun with LiDAR

LiDARPlotView Just a little animated GIF showing LiDAR from the North Kaibab. Notice the penetration achieved and the resolution it provides for picking out large trees, small bushes, and topography.
Facebooktwittergoogle_plusredditpinterestlinkedinmailFacebooktwittergoogle_plusredditpinterestlinkedinmail

Fun with LiDAR

LiDAR
So I finally got my hands on the North Kaibab LiDAR acquisition and it’s huge (~2TB). I was able to find a few minutes to play with it (in ENVI) to produce this simple figure of the landscape surrounding Matt Tuten’s research sites. Cool beans…

Facebooktwittergoogle_plusredditpinterestlinkedinmailFacebooktwittergoogle_plusredditpinterestlinkedinmail