R Lidar Voxel Calculation – Using lidR to Calculate the Number of Null Voxels

laslidarlidrrvoxel

I have a number of LAS files (750+). I am using the package lidR to extract different metrics. I am using the following function to calculate the number of voxels in each files that contain a LiDAR return.

files <- list.files(path= "path/to/files")

O = lapply(files, function(x) { 
  las <- readLAS(x, select = "xyzc", filter = "-drop_z_below 0 -drop_class 1 2 6 7 18")
  las_voxel <- voxelize_points(las, 1.0) # Creating 1.0^3 voxels
  coverage <- npoints(las_voxel) # number of voxels containing >= 1 LiDAR return
  return(data.frame(file = x, coverage)
})

I wish to calculate the percentage of voxels that contain vegetation. I am wondering if it is possible to calculate the number of null voxels – i.e., those with 0 LiDAR returns – within this function? I have come stuck when trying to use voxel_metrics: Using lidR to calculate null voxels

Best Answer

You must define better what you mean by number of voxels with 0 points. Because there are virtually an infinite of them above the point cloud. So I guess you want to set a threshold. But what is this thresholds?

The highest point in your 750 files?
The highest point of each file?
An arbitrary value?

Iĺl suppose you want the solution number 2. You can adapt for your needs. First you must know that voxel_metrics() has an argument all_voxels (the question you linked is older than this option)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)
vox_size = 8
vm <- voxel_metrics(las, ~length(Z), vox_size, all_voxels = TRUE) 
sum(is.na(vm$V1))
#> 1622

But I don't think it is an optimal idea. With a resolution of 1m it will likely generate more voxels than points. You can simply compute the number of voxels

las_voxel <- voxelize_points(las, vox_size)
coverage <- npoints(las_voxel)
nX = diff(range(las_voxel$X))/vox_size + 1
nY = diff(range(las_voxel$Y))/vox_size + 1
nZ = diff(range(las_voxel$Z))/vox_size + 1
nV = nX*nY*nZ
nV - coverage
#> 1622

Related Solutions

R LiDAR – How to Filter Above Ground LiDAR Points Using Catalog in lidR

In your case the filter you are using is a simple one: Classification != 2. And you don't need the ground points at all. You are better to use a streaming filter and a streamed processing.

ctg <- catalog("/...2015TestGroup")
opt_chunk_size(ctg) <- 0
opt_chunk_buffer(ctg) <- 0
opt_output_files(ctg) <- ".../Outputs/2015nonground/{XLEFT}_{YBOTTOM}_2015"
opt_filter(ctg) <- "-drop_class 2"
new_ctg <- catalog_retile(ctg)

catalog_retile is a streamed function. It means that not a single point is loaded at the R level. Each point is read, filtered and written immediately without calling any R code and without loading any memory. It is thus way faster than using the R function lasfilter.

That being said, you are better to use las2las from LAStools for this kind of job. It does exactly the same actually but is more efficient and simple than calling more or less equivalent code from R.

lidR R Package – Fixing Error “$ Operator Not Defined for S4 Class”

You cannot use chunk@data$Number.of.1st.return. Chunk is a LAscluster. This class is not documented and you know nothing about it. You only know that you can write las = readLAS(chunk) and use st_bbox(chunk). I also fixed many typos, error, missing variables and so on.

calculate_las_metrics <- function(chunk) 
{
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  
  
  # Calculate point density for las file
  total_point_density <- density(las)
  
  # extract unique classifications for las file
  classes <- unique(las$Classification)
  
  # extract other information/descriptive data
  num_total_points <- npoints(las)
  #num_returns_total <- sum(las$NumberOfReturns) does not make sense. 
  num_of_1st_return <- las$`Number of points by return`[1]
  
  # Create data frame of stats
  metrics_df_total <- data.frame(
    Classification = classes,
    Num_Points = num_total_points,
    #Num_Returns = num_returns_total,
    #Points_Per_Return = points_per_return_total,
    num_of_1st_return = num_of_1st_return,
    Point_Density = total_point_density
  )
  
  return(metrics_df_total)
}

Best Answer

Related Solutions

R LiDAR – How to Filter Above Ground LiDAR Points Using Catalog in lidR

lidR R Package – Fixing Error “$ Operator Not Defined for S4 Class”

Related Question