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?

  • 2
    What have you tried so far and what are you specifically having trouble with? We don't write scripts here, but we are happy to help troubleshoot.
    – MaryBeth
    Mar 8 '16 at 17:53
  • FUSION is coded in C++ and the only scripting available is calling existing exe functions through a DOS bat file. The cylinder subset is available in the GUI as a visualization tool but you will have to do some digging through the users manual to see if there is a standalone exe that would facilitate your needs. Mar 8 '16 at 18:39
  • 2
    Welcome to gis.stackexchange! Please note that a good question on this site is expected to show some degree of research on your part, i.e. what you have tried and - if applicable - code so far. For more info, you can check our faq
    – underdark
    Mar 8 '16 at 19:31

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:

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,",
             "] c:\\input\\file1.las c:\\output\\file1.csv",sep=""))

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,",
               "] 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.


+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)


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(cylinders, col="red", add= TRUE)

Assign cylinder ID's to subset lidar data

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

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)
  • 1
    +1 Very cool. How efficient is this process from a computational standpoint compared to running FUSION command line tools?
    – Aaron
    Mar 12 '16 at 21:04
  • 2
    @Aaron, honestly I don't know but I have a function that bins a point cloud and calculates several liar height and intensity metrics, returning raster output. I have processed 180 tiles with ~20M points each and it took about two hours so, not speedy. The advantage is that you can customize the output and stay in a modeling environment for additional analysis. Mar 12 '16 at 21:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.