Retrieving LiDAR metrics from relative vertical layers of point cloud?

Question

I need to obtain the proportion of all returns within 5 equally spaced height bins between the cutoff (1 meter) and the maximum height value. I need to get this proportion of all returns for 1000 cylinders with diameter of 28 m (I have the center coordinates of each cylinder).

I am planning to use FUSION to do this task. It seems that PolyClipData will work to clip the cylinders if I generate a shapefile file buffering my points. Then, CloudMetrics (strata:[#,#,#,...) would compute the returns in various height strata. However, the example in the user's manual shows that strata will be fixed (from the default strata) i.e. 0.15, 1.37, 5, 10, 20, 30.

What I want is to have 5 equally spaced height bins between the cutoff and the maximum height value, so the height of the strata will vary according to the highest return. How can specify this function in the switch strata?

Answer 1

It is possible to run CloudMetrics as a first round to get the maximum height, and then, run it again specifying strata with the desired intervals. The trick is to automate such process.

I did it using the software R, but you can do it using other software (I guess even directly from the command prompt window).

I assume you already have one cylinder per lidar file, and also that the point clouds are normalized. If not, you will find help in:

• How to clip LAS data using shapefile polygons and open source software?
• Normalizing point cloud data

First, let's exemplify how it would be considering just one lidar file. Here is the code, properly commented.

#Call Cloudmetrics .exe file from R. Change file paths and file names accordingly.
system("c:\\fusion\\cloudmetrics c:\\input\\file1.las c:\\output\\file1_temp.csv")

#Import CloudMetrics output in R.
cloudmetrics1 = read.csv("c:\\output\\file1_temp.csv",sep=",")

#Get maximum height from CloudMetrics output in 15th column. Perhaps you'll need to tweak the column number. 
max_height = cloudmetrics1[1,15]
#Set the vertical layers range.
layer_range = max_height/5
#Get strata values to be used in the second CloudMetrics run.
strata_value = c(layer_range,2*layer_range,3*layer_range,4*layer_range)

#Call Cloudmetrics again, but now using the switch 
#'strata' to specify which breaks for the vertical layers.
system(paste("c:\\fusion\\cloudmetrics /strata:[0,1,",
             strata_value[1],",",strata_value[2],",",strata_value[3],",",strata_value[4],
             "] c:\\input\\file1.las c:\\output\\file1.csv",sep=""))
#Done.

Now, the same code, but with a for loop that will automate the process. Note that the name of the input files need to have some kind of ID to be identified inside the loop.

#Loop through the las files with names 'file#.las' where # is equal some integer number.
#The example below loop through two las files named file1.las and file2.las, respectively.

for (i in c(1:2)){

  system(paste("c:\\fusion\\cloudmetrics c:\\input\\file",i,".las c:\\output\\file",i,"_temp.csv", sep=""))

  cloudmetrics1 = read.csv(paste("c:\\output\\file",i,"_temp.csv",sep=""),sep=",")
  max_height = cloudmetrics1[1,15]
  layer_range = max_height/5
  strata_value = c(layer_range,2*layer_range,3*layer_range,4*layer_range)

  system(paste("c:\\fusion\\cloudmetrics /strata:[0,1,",
               strata_value[1],",",strata_value[2],",",strata_value[3],",",strata_value[4],
               "] c:\\input\\file",i,".las c:\\output\\file",i,".csv",sep=""))
}

The csv outputs will contain the number and proportion of all returns in each layer specified.

You may continue using R to organize the generated statistics in one table, and continue the processing workflow from there.

Answer 2

+1 to @Andre Silva for providing a solid workflow using Fusion. Here is a solution that cuts out Fusion and allows you to keep the analysis in R.

Add required libraries and simulate some lidar data with height (0-40)

library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(spatialEco)

e <- extent(553462.6, 554070.4, 5225868, 5226352)
  e <- as(e, "SpatialPolygons")
lidar <- spsample(e, 500000, type="random")  
  lidar <- SpatialPointsDataFrame(lidar, data.frame(height = 
                            runif(length(lidar), 0, 40)))  

Create some sampling cylinders with 10m radius

samp.pts <- spsample(e, 10, type="regular")
  samp.pts <- SpatialPointsDataFrame(samp.pts, data.frame(ID=1:length(samp.pts)))   
cylinders <- gBuffer(samp.pts, byid=TRUE, id=samp.pts$ID, width = 28, quadsegs=10)
  plot(e)
    plot(cylinders, col="red", add= TRUE)

Assign cylinder ID's to subset lidar data

lidar.sub <- point.in.poly(lidar, cylinders) 
 head(lidar.sub@data)

Write a function to stratify heights into 5 bins

strata <- function(x) { 
  bks <- max(x) / 5
  ht.strata <- c(   
    length( x[x > 1 & x <= bks] ) / length( x[x > 1] ), 
    length( x[x > bks & x <= 2 * bks] ) / length( x[x > 1] ), 
    length( x[x > 2 * bks & x <= 3 * bks] ) / length( x[x > 1] ), 
    length( x[x > 3 * bks & x <= 4 * bks] ) / length( x[x > 1] ), 
    length( x[x > 4 * bks] ) / length( x[x > 1] ) )
    names(ht.strata) <- c("strata1", "strata2", "strata3", "strata4", "strata5")    
    return( ht.strata )
}     

We can then use tapply to pass the strata function to each cylinder

tapply(lidar.sub$height, lidar.sub$ID, FUN = strata)