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

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

  • Steve

    Thanks for sharing this Andrew, really nice to call the FUSION exicutables from R.

    Did you mean to add the rgeos package to the script? I got a dissolve error somewhere after line 55. I was also having trouble getting rLiDAR to install in v. 3.2.2, but all seems to be running well in v 3.1.2. It may have been a dependency issue, but added rgeos anyway.

  • bi0m3trics

    rgeos is a “Suggests” on raster and sp, so I always have it installed. Guess I must have missed it. I added it after line 54.

    Did you install rLiDAR using the source? i.e., install.packages(“rLiDAR”, type=”source”)? I had a little trouble getting it to run on 3.2.2 but installing from the source on r-forge fixed it for me… install.packages(“rLiDAR”, repos=”http://R-Forge.R-project.org”) Also, the dissolve isue you had is likely with CHMsmoothing (line 9 of the second chunk). Everything before that is just plotting. Were you running it on your data or the sample I provided?

    Really good seeing you this last week and I really appreciate your insights on using WorldView2. It helped alot…

  • wuodasembo

    Hallo,
    Am trying to run this script but am experiencing a problem. the error seems to be the directories. I am using the example data here and I have placed it in the inDirectory, and am having an error result and a warning message below:

    GroundFilter v1.20 (FUSION v3.42) (Built on Mar 28 2014 09:44:05) DEBUG
    Using LASzip.dll (c) Martin Isenburg–Rapidlasso for LAS/LAZ file access
    LASzip.dll V2.3 r0 (build 150402)
    Command line: C:\Fusion\GROUND~1.EXE /gparam:0 /wparam:1 /tolerance:1 /iterations:10 D:/LAS in R with FUSION/outputDirectory/Example_GroundPts.las 1 D:/LAS in R with FUSION/inputDirectory/Example.las
    Run started: Fri Feb 19 15:39:46 2016
    Error for R: Can’t open file
    Error for with: File does not exist
    Error for FUSION/outputDirectory/Example_GroundPts.las: File does not exist
    Error for 1: File does not exist
    Error for D:/LAS: File does not exist
    Error for in: File does not exist
    Error for R: Can’t open file
    Error for with: File does not exist
    Error for FUSION/inputDirectory/Example.las: File does not exist
    Run finished: Fri Feb 19 15:39:46 2016 (elapsed time: 0 seconds)
    ***There were errors during the run
    Done
    Warning message:
    running command ‘C:/Fusion/groundfilter.exe /gparam:0 /wparam:1 /tolerance:1 /iterations:10 D:/LAS in R with FUSION/outputDirectory/Example_GroundPts.las 1 D:/LAS in R with FUSION/inputDirectory/Example.las’ had status 2

    Any assistance on what i might have done wrong?

    • bi0m3trics

      The problem seems to be with the directory name you’ve chosen… Fusion doesn’t like spaces in folder or file names. Try changing your folder from “LAS in R with FUSION” to “LASInRWithFUSION” and let me know if that fixes your problem… basically, the error “Error for R: Can’t open file: says it can’t find the file “Example_GroundPts.las” the way you have it set up.

      When it runs the output in R will look like this:
      GroundFilter v1.20 (FUSION v3.42) (Built on Mar 28 2014 09:44:05) DEBUG
      Using LASzip.dll (c) Martin Isenburg--Rapidlasso for LAS/LAZ file access
      LASzip.dll V2.2 r0 (build 140907)
      Command line: C:\Fusion\GROUND~1.EXE /gparam:0 /wparam:1 /tolerance:1 /iterations:10 C:/Dropbox/LiDAR_Prj/outputDirectory/Example_GroundPts.las 1 C:/Dropbox/LiDAR_Prj/inputDirectory/Example.las
      Run started: Sat Feb 20 10:39:10 2016
      Data files:
      C:/Dropbox/LiDAR_Prj/inputDirectory/Example.las: 194244 points
      Ground point file produced:
      C:/Dropbox/LiDAR_Prj/outputDirectory/Example_GroundPts.las Feb 20, 2016 @ 10:39 AM
      96129 ground points
      Run finished: Sat Feb 20 10:39:11 2016 (elapsed time: 1 seconds)
      Done

  • Stevo

    >>loc<-FindTreesCHM(sCHM, fws, minht)

    Error in `colnames<-`(`*tmp*`, value = c("x", "y", "height")) :
    length of 'dimnames' [2] not equal to array extent

    • Stevo

      Why does this error persist?

      • bi0m3trics

        You didn’t provide enough information for me to say really, but if I had to guess it’s because of how the release code for the FindTreesCHM was written. At this point, I’ll say I’m not the maintainer of this package, but I do know a thing or two and I’ve encountered this before… Specifically, the SpatialPoints call in that function needs the projection string defined otherwise the next line’s call to over may have issues. Replace the FindTreesCHM function with the following version and it should work….

        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))
            htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
            treeList <- cbind(XYmax, htExtract)
            colnames(treeList) <- c("x", "y", "height")
            return(treeList)
        }

        Basically, adding the “proj4string = crs(chm)” to the call to SpatialPoints may fix your problem. Also, this fix is commented out in the Github repo for the package and I find it’s only an issue when using it on large datasets. Let me know if it works….

        • Stevo

          Sorry for not explaining at length. These are the steps I took:

          1. I used the Lidar.R script and the Example.las file downloaded from this webpage.
          2. I installed all the necessary library packages required.
          3. The script ran fine till when it reached at this line:


          loc<-FindTreesCHM(sCHM, fws, minht)

          where this (below) error is called irregardless of the R studio installation
          I used or the number of trial runs I attempted


          Error : Error in `colnames<-`(`*tmp*`, value = c("x", "y", "height")) :
          length of 'dimnames' [2] not equal to array extent

          4. I looked at your solution and ran a trace on the function above i.e.


          trace(FindTreesCHM(sCHM, fws, minht))

          which described my ORIGINAL FindTreesCHM function as:


          function(chm, fws, minht)
          {
          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))
          htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
          treeList <- cbind(XYmax, htExtract)
          colnames(treeList) <- c("x", "y", "height")
          return(treeList)
          }

          This is the same as the solution you posted above. I tried one more solution…. I changed the function name and embedded the function posted by you in the script i.e.


          trialfunct <- 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))
          htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
          treeList <- cbind(XYmax, htExtract)
          colnames(treeList) <- c("x", "y", "height")
          return(treeList)
          }

          loc <- trialfunct(sCHM, fws, minht)

          but still the error persisted. Is there another solution?

          • bi0m3trics

            Interesting… for the life of me I can’t reproduce your errors using the code I’ve provided here on my website. However, I was able to get it using a larger dataset… Also, I’m not sure how the return from TRACE on the original FindTreesCHM function could include the SpatialPoints(xyFromCell(setNull, Which(setNull == 1, cells = TRUE)), proj4string = crs(chm)) as the versions (0.1) on CRAN, rforge, and in gitHub repository all read XYmax < - SpatialPoints(xyFromCell(setNull, Which(setNull==1, cells=TRUE))). However, that's clearly not the problem...

            Based on what you've told me, I can't help but wonder if the problem lies in the call to cbind when creating treeList... Does it work if you replace treeList < - cbind(XYmax, htExtract) with treeList< -cbind(as.matrix(XYmax@coords),htExtract) in the modified version of FindTreesCHM?

            Another option would be to pass sCHM to chm and see if the code internal to the FindTreesCHM function (just running everything inside the function locally) is building treeList correctly.

  • Solution:

    I realised what the mishap was….the segment:

    treeList <- cbind(XYmax, htExtract)
    colnames(treeList) <- c("x", "y", "height")

    The XYmax is an object of the SpatialPoints class….meaning “cbind” only creates TWO columns up here –
    1. A spatial data column
    2. htExtract column

    but colnames is set to change THREE columns…..

    Therefore for all who will have a problem with this segment use :

    treeList <- cbind(coordinates(XYmax), htExtract)

    instead of

    treeList <- cbind(XYmax, htExtract)

    and the whole code will run smoothly. Haven’t attempted the “as.matrix” part yet, but maybe it will do also although “coordinates” only takes the X and Y columns directly.

    Do the following steps as discussed previously :

    1. write the “FindTreesCHM()” in the code script itself i.e. copy paste this in your script before calling it in the loc object and everything will run smoothly:

    #first

    function(chm, fws, minht)
    {
    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))
    htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
    treeList <- cbind(coordinates(XYmax), htExtract)
    colnames(treeList) <- c("x", "y", "height")
    return(treeList)
    }

    #second

    loc<-FindTreesCHM(sCHM, fws, minht)

    add the first section above the second segment in your code.

    Cheers.
    #interesting piece of code for the CHM

  • instead of


    function(sCHM, fws, minht){....}

    change “function” to


    FindTreesCHM(sCHM, fws, minht){....}

    or at least whichever function name you used.

    • bi0m3trics

      Glad you got it working! I edited the original post to reflect these changes and once I get around to it I’ll shoot the package maintainer a note. I honestly didn’t think anyone had seen this post…

      The xxx< -function(){...} part was my bad - maybe I wasn't clear... the FindTreesCHM< -function(chm, fws = 5, minht = 1.37){...} creates a custom function (slightly edited in this case) that will replace the original FindTreesCHM(), only if you run it after the library has been loaded (i.e., it will override the original FindTreesCHM()). So calling it using funciton() won’t work. Standard custom function creation stuff…

      Lastly, my solution with the “as.matrix” part basically does the same thing as your solution for cbind (without an addition call for the coordinates() function from the sp package). Old habits… So given this, there’s no need to proceed as your solution is fine and is actually easier to understand when reading the code.

      Cheers!

  • la

    thanks a lot for discussing and solving this problem!!!

  • Rodrigo

    Dear,

    Thanks for sharing this code. I am not an expert and it has worked with some data, but now it has failed to:

    canopy <-ForestCAS (schm, loc, maxcrown, exclusion)

    When I run it gives the following error:

    ……………………………………………………………………………………………suggested minimum tolerance: 1.86265e-08

    Error in points2grid (points, tolerance, round):
       Dimension 2: coordinate intervals are not constant

    Please appreciate your help.

    Regards,

    • bi0m3trics

      I’m sorry, but I just don’t follow… If I have to guess, I would say your problem is the call to SpatialPixelsDataFrame in DF2raster. So have a look at the DF2raster function within ForestCAS and maybe uncomment the call to coordinates(spP). I think it’s the only place where points2grid is called…

      • jimichili

        I had the same error. I changed the tolerance value in the call to SpatialPixelsDataFrame in function DF2Raster and now it is working fine
        ———
        m <- suppressWarnings(SpatialPixelsDataFrame(points = spP[c("x", "y")], data = spP,tolerance=0.0001))
        ———
        I temporarily removed the suppressWarnings() and noticed that the point2grid error occurred in a case where the following warning message is issued :
        ———
        Warning message:
        In points2grid(points, tolerance, round) :
        grid has empty column/rows in dimension 2
        ———
        But I do not know if this has any relationship with the point2grid error.

        • bi0m3trics

          Good to know and thanks for posting. According to rforge rLiDAR hasn’t had a commit in 17 months. In the mean time, I’ve been using the hell out of Jean-Romain Roussel’s lidR package. That and I bit the bullet and implemented my own version of Li et al’s (2012) algorithim in c++…

  • stukrause

    Thanks for the FindTreesCHM fix! Any Ideas in how I should citate this for my master thesis?

    • bi0m3trics

      I’d just cite the rLiDAR package and if you really need to you can cite the specific fix for FindTreesCHM (i.e., you think it is crucial for replication) I think a personal communication with me, Andrew Sanchez Meador, would be sufficient… I’d be interested in seeing what you did for your thesis.
      Andrew

  • Peter

    Hi Andrew,thanks for the impressive blog!

    In this step:
    # 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=” “))

    the DTM2ASCII took a very long time to finish, e.g., my las file has 160 MB(5 million points) but it has been one day and it is not done. Are there some other tools (or ways) can make this process faster?

    Look forward to hearing from you.

    Regards,

    Peter(pingyang.whu@gmail.com)

    • bi0m3trics

      I typically chunk the data using buffered polygons (e.g., stands) and process it using the foreach and doParallel packages. I then clip the results to the original polygon(s) before merging them into one output file. Maybe I’ll get around to posting on that in the next couple of weeks… once things slow down.

  • Hi Andrew. I am the author of the rLiDAR package.

    First all, thanks for have this group to discuss about LiDAR, FUSION and rLiDAR package.

    It is very nice to know that someone is using the rLiDAR. Also, thanks for fixing the bugs! I don’t have sure, but maybe I will upload a new version of it in the rforge by the end of the year. If you are interested to be part of the project, I can send you a rforge invitation!

    This is a example of how to cite the package: 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, .

    Also, this is on publication that came out about it: rLiDAR:http://www.tandfonline.com/doi/full/10.1080/07038992.2016.1196582?scroll=top&needAccess=true

    I am interested to know more about your implementation of Li et al’s (2012) algorithm for individual tree detection. Are you able to share it? Thanks!

    • bi0m3trics

      Wow Carlos! It’s nice to finally meet you and I see that you’ve been busy. Last time I looked at your web application list, there was only one or two applications on there. Now there’s ten!

      In case you didn’t know, Nick pointed me to your work last time he was in town, and I’d say I’ve made good use of it in the past year. rLiDAR is a really nice piece of work and I’d be more than happy to contribute bug fixes and some minor documentation corrections. I’ve used it on a few projects, a couple of which I hope will be submitted in the next few months and I’m really pleased that this post has received so much attention. The comment section has especially received quite a bit of attention lately and I’ve got a couple of follow-ups I need to post – but life and work keeps getting in the way.

      My implementation of the Li et al’s algorithm is all in c++ and I’ve only compiled it for windows (it was for one of the above projects I mentioned). I’d be happy to share it and maybe this will prompt me to finally getting around to wrapping it in rcpp or exporting the functions properly (so it can linked). I’m afraid it might require some rewriting to facilitate its use in R. The project I was working on had a quick timeline, so I just built executables and called it via commandline using system and doparallel… I know. I’m a hack!

      Thanks for stopping by and thanks for the link – I’ll have to check it out tomorrow once I’m back in the office. Also, I added a footnote providing readers with your suggested citation…

      • JR

        Hi Andrew. I’m the author of lidR package.

        I began to implement Li et al’s algo in R. Very inefficient as R is not design for this kind of stuff. I plan to write it in C++ within Rcpp. But if you already did it I’m interested in your version of the algorithm. I can try to compile it for GNU/linux and wrap it with Rcpp. Would you like to share it?

        Thanks

        Jean-Romain

  • Leo A. Yanez

    if you want to export the thiessen polys and tree points to shape files this is how to do it :
    #rgdal is superb for this kind of thing. http://www.gdal.org/ has most of the information you could every want on what formats are supported.

    #In this case you want the writeOGR function
    # this will create a shapefile for your thiessen poly within the working directory
    library(rgdal)
    #test convert point shape file
    coord = loc
    coorddf <- SpatialPointsDataFrame(coord, data = data.frame(loc))
    writeOGR(boundaryTrees, dsn = 'c://fusion//testoutput', layer = 'stand_thiessen3', driver = "ESRI Shapefile")
    writeOGR(coorddf, dsn = 'c://fusion//testoutput', layer = 'stand_loc_xy_heights', driver = "ESRI Shapefile")

  • Steve Sesnie

    Hi Andrew,

    Working with this script a little bit today, and thought others may want to batch process multiple tiles at once :

    #Create a list of .las files that can be read by FUSION
    startingDir <- "/pathtoyourlasfiles/"
    #Make a list
    lasFiles <- list.files(startingDir, pattern = ".las", full.names = TRUE)
    head(lasFiles) #check the list
    #Create a proper list for FUSION
    lasFiles <- gsub("/", "\\", lasFiles, fixed=TRUE)
    head(lasFiles) #This list will run in FUSION
    #Write to a .txt file
    write(lasFiles, "/pathtoyourlasfiles/lasFiles.txt")

    #For running FUSION .exe
    #Canopy height model example for multiple tiles using a list (lasFiles.txt)
    system(paste(file.path("C:","Fusion", "canopymodel.exe"),
    paste("/class:3,4,5, /ground:",file.path(mainDir, outDir, "Example_DTM.dtm"), sep=""),
    file.path(mainDir, outDir, "Example_CHM.dtm"),
    "3 F F 2 0 2 2",
    file.path(mainDir, inDir, "lasFiles.txt"),
    sep=" "))

    #Note that this is example uses a State Plane projection in feet
    #I've also added /class: if working with classified returns
    #This works great for building extraction if /class:6 is available
    #Memory will be a limiting factor if a large number of tiles are selected

    Hope you are well!!

Leave a Reply

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>