Custom voxel metrics catalog processing

[RCBlackburn]

Hello,

I have been struggling with applying custom voxel functions to a catalog. When I apply them (only on a catalog), differences appear in some of the chunks (see attached image). This does not happen when applying .std_metrics. I'm not looking for anyone to debug all the code (shown below), but I am throwing it out there in case anyone can tell right off the bat why applying a voxel function to a catalog creates different ranges of those metrics in each chunk.

Any thoughts are much appreciated!!!

## write custom function

## function to get points per voxel
vox_mt <- function(z)
{
  SVi = length(z)
  
  metrics =list(
    SVi = SVi
  )
  return(metrics)}

## function to get unique voxel metrics per voxel
vox_mt2 <- function(vox, res)
{
  fullvox <- data.table(vox, 
                        id_xyz = paste0(vox$X,"-", vox$Y,"-",vox$Z))
  fullvox$SVi[is.na(fullvox$SVi)] <- 0
  
  fullvox$FR_SVi <- fullvox$SVi/sum(fullvox$SVi)

  unique_Z <- unique(fullvox$Z)
  
  
  point_blow <- fullvox[,as.list(lapply(unique(Z), function(x){
    sum(SVi[Z < x])})),
    by = list(X,Y)]
  colnames(point_blow)[3:ncol(point_blow)] <- paste0(unique(fullvox$Z))
  point_blow <- point_blow %>%
    tidyr::pivot_longer(cols = 3:ncol(point_blow), names_to = "Z", values_to = "npoints_below")
  
  point_blow$id_xyz = paste0(point_blow$X,"-",point_blow$Y,"-",point_blow$Z)
  fullvox <- fullvox %>% left_join(point_blow[,c(4,5)], by = "id_xyz")
  
  point_above <- fullvox[,as.list(lapply(unique(Z), function(x){
    sum(SVi[Z > x])})),
    by = list(X,Y)]
  
  colnames(point_above)[3:ncol(point_above)] <- paste0(unique(fullvox$Z))
  point_above <- point_above %>%
    tidyr::pivot_longer(cols = 3:ncol(point_above), names_to = "Z", values_to = "npoints_above")
  
  point_above$id_xyz = paste0(point_above$X,"-",point_above$Y,"-",point_above$Z)
  fullvox <- fullvox %>% left_join(point_above[,c(4,5)], by = "id_xyz")
  
  
  fullvox$FR_Di = fullvox$npoints_above/sum(fullvox$SVi) 
  
  metrics =data.frame(
    X = fullvox$X,
    Y = fullvox$Y,
    Z = fullvox$Z,
    SVi = fullvox$SVi,
    FRDi = fullvox$FR_Di,
    PDi = fullvox$npoints_below,
    PDi_above = fullvox$npoints_above
  )
  names(metrics)[4:ncol(metrics)] <- paste0(names(metrics)[4:ncol(metrics)])
  return(metrics)
}

## wrap up above functions to get a voxelized pint cloud to use in pizel metrics
ind_vox <- function(las, res = res){
  vox <- lidR::voxel_metrics(las, func = vox_mt(Z), res = res, all_voxels = TRUE)
  vox_2 <- vox_mt2(vox, res = res)
  out <- LAS(vox_2, header = las@header, crs = las@crs, index = las@index )
}

## write summary function for pixel_metrics() function in lidR
vox_sum_raster <- function(SVi, FRDi, PDi, PDi_above, vox_res){
  las_vox <- data.frame(SVi,FRDi, PDi, PDi_above)
  means <- apply(las_vox, 2, mean, na.rm = TRUE)
  names(means)<- paste0(names(means), "_", vox_res, "_", "mean")
  meds <- apply(las_vox, 2, median, na.rm = TRUE)
  names(meds)<- paste0(names(meds), "_", vox_res, "_", "med")
  var <- apply(las_vox, 2, var, na.rm = TRUE)
  names(var)<- paste0(names(var), "_", vox_res, "_", "var")
  sd <- apply(las_vox, 2, sd, na.rm = TRUE)
  names(sd)<- paste0(names(sd), "_", vox_res, "_", "sd")
  cv <- apply(las_vox, 2, function(x) {sd(x, na.rm = T)/mean(x, na.rm = T)})
  names(cv)<- paste0(names(cv), "_", vox_res, "_",  "cv")
  IQR <- apply(las_vox, 2,IQR, na.rm = T)
  names(IQR)<- paste0(names(IQR), "_", vox_res, "_",  "IQR")
  skew <- apply(las_vox, 2, function(x) {(sum((x - mean(x, na.rm = T))^3,na.rm = T)
                                          /length(x))/(sum((x - mean(x, na.rm = T))^2,na.rm = T)
                                                       /length(x))^(3/2)})
  names(skew)<- paste0(names(skew), "_", vox_res, "_", "skew")
  kurt <- apply(las_vox, 2, function(x) {length(x)*sum((x- mean(x, na.rm = T))^4, na.rm = T)/
      (sum((x - mean(x, na.rm = T))^2, na.rm = T)^2)})
  names(kurt)<- paste0(names(kurt),  "_", vox_res, "_",  "kurt")
  pct_fill_vox <- data.frame(pct_fill_vox = as.numeric(sum(las_vox[,1] > 0 ))/ length(las_vox[,1]))
  names(pct_fill_vox) <- paste0(names(pct_fill_vox),"_",vox_res)
  
  data <- data.frame(c(means, meds, var, sd, cv, IQR, skew, kurt, pct_fill_vox))
  
  return(data) 
}

all_layers_catalog = function(cluster)
{ 
  #pre-processing
  las = readLAS(cluster)
  las = las_update(las) 
  if (lidR::is.empty(las)) return(NULL) 
  las = filter_poi(las, Classification != LASNOISE)
  if (lidR::is.empty(las)) return(NULL)
  las = normalize_height(las, knnidw(k = 8, p = 2)) 
  las = filter_poi(las, Z < 50 & Z >= 1.37 ) 
  if (lidR::is.empty(las)) return(NULL)
  las <- decimate_points(las, random_per_voxel(res = 1, n = 8))
  if (lidR::is.empty(las)) return(NULL)
  
  # cloud metrics
  gmetrics <- pixel_metrics(las, .stdmetrics_z, res = 10)
  
 
  # voxel metrics
  las_vox3 <- ind_vox(las, res = 2) 
  vmetrics_3 <- pixel_metrics(las_vox3,
                              func = vox_sum_raster(SVi, FRDi, PDi, PDi_above, vox_res = 2),
                              res = 10)
  
  vmetrics_3 <- resample(vmetrics_3, gmetrics)
  tile0_ext = lidR::ext(las)
  gmetrics <- crop(gmetrics,tile0_ext)
  vmetrics_3 <- crop(vmetrics_3, tile0_ext)
  all_metrics <- c(vmetrics_3, gmetrics)
  
  return(wrap(all_metrics))
}

## load in catalog and apply functions to it
ctg <- readLAScatalog("path-to-ctg")
opt_chunk_buffer(ctg) <- 10 
opt_chunk_size(ctg) <- 250    
opt_select(ctg) <- "xyzicr"      
opt_wall_to_wall(ctg) <- TRUE
opt_chunk_alignment(ctg) <- c(0,0)
opt_progress(ctg) <- TRUE
opt_stop_early(ctg) <- TRUE

layer_list <- catalog_apply(ctg, all_layers_catalog)

The result regardless of what data I use looks something like below. Note the clear chunk outlines in the yellow oval:

[Jean-Romain]

Please provide a minimal reproducible example including data (data shipped with the package if possible) and removing as much code as possible to get a meaningful output but avoid sharing hundreds of lines of code.

At first glance, I see one more or less obvious problem. Your voxel resolution if not a multiple of your pixel resolution. Thus, the voxel grid is not aligned with the pixel grid, and you missed some data by computing a sub grid inside a grid that does not contain the adequate data for the sub grid. But, as I said, the code is too complex to grasp the problem just reading the code.

[RCBlackburn]

Thanks for the input! I think you are correct that alignment is at least one of the issues. However, even changing the voxel resolution to 2 produces some differences in chunks. These differences are even more apparent however when I play with the opt_chunk_alignment() option (see bottom image). Note where chunks overlap only a portion of the catalog, the voxel values are generally higher. Most of the catalogs I work with are not square if that means anything.

I have run an example on the catalog provided in the lidR book. using the functions provided above. However, I changed the resolution to 2 in the ind_vox() and vox_sum_raster() functions. I also edited the original post so the functions can be copied and pasted to produce the results from the following code:

## Load in catalog example from https://r-lidar.github.io/lidRbook/engine.html
ctg <- readLAScatalog("ctg_example/")
opt_chunk_buffer(ctg) <- 10 
opt_chunk_size(ctg) <- 250    
opt_select(ctg) <- "xyzicr"      
opt_wall_to_wall(ctg) <- TRUE
opt_chunk_alignment(ctg) <- c(0,0)
opt_progress(ctg) <- TRUE
opt_stop_early(ctg) <- TRUE

plot(ctg, chunk = TRUE)

layer_list <- catalog_apply(ctg, all_layers_catalog)

rasters <- lapply(layer_list, terra::unwrap)
rasters <- do.call(terra::merge, rasters)

plot(rasters[[2]])

Here is the output of the above code with some chunk borders circled in yellow (voxel res = 2; chunk alignment is c(0,0)):

Here is the output of the above code with a different alignment. Some chunk borders are circled in yellow (voxel res = 2; chunk alignment is c(70,0)):

[Jean-Romain]

As I said earlier, please share a minimal reproducible example, including data. I'm keen to help, but I won't spend hours digging and trying to reproduce myself with my own data. I need something simple I can copy paste to run and reproduce.

[bi0m3trics]

@RCBlackburn - are you saying that the above code you provided produces the issue (pictured) with the example .laz files from the book (which has all the .laz files to build the catalog)? Specifically, the nine files in section 14.3 read as a catalog produce what you circled?

[RCBlackburn]

@bi0m3trics That is correct. I wanted to provide a sample using the same chunk size that I found the issue with.

@Jean-Romain here is some code with the Megaplot.laz in the library. It is harder to see with a smaller catalog with a pixel_metrics() resolution of 10m. However, if you set the pixel_metrics() resolution to 4m it is more visible. I've also included the GitHub repository install for the voxel functions. There are minor differences with the original voxel function posted in this thread and the lidRattR package, but not regarding this issue.

devtools::install_github("RCBlackburn/lidRattR")
library(lidRattR)
library(lidR)
library(terra)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile)

opt_chunk_buffer(ctg) <- 10 
opt_chunk_size(ctg) <- 100    
opt_select(ctg) <- "xyzicr"      
opt_wall_to_wall(ctg) <- TRUE
opt_chunk_alignment(ctg) <- c(0,0)
opt_progress(ctg) <- TRUE
opt_stop_early(ctg) <- TRUE

plot(ctg, chunk = TRUE)
all_layers_catalog = function(cluster)
{ 
  las = readLAS(cluster)
  las = las_update(las) 
  if (lidR::is.empty(las)) return(NULL) 
  las = filter_poi(las, Classification != LASNOISE)
  if (lidR::is.empty(las)) return(NULL)
  las = normalize_height(las, knnidw(k = 8, p = 2)) 
  las = filter_poi(las, Z < 50 & Z >= 1.37 ) 
  if (lidR::is.empty(las)) return(NULL)
  las <- decimate_points(las, random_per_voxel(res = 1, n = 8))
  if (lidR::is.empty(las)) return(NULL)
  
  ### voxel metrics
  las_vox <- ind_vox(las, res = 2) 
  vmetrics <- pixel_metrics(las_vox,
                              func = vox_sum_raster(SVi, FRDi, PDi, PDi_above, vox_res = 2),
                              res = 10)
  
  tile0_ext = lidR::ext(las)
  vmetrics <- crop(vmetrics, tile0_ext)
  
  return(wrap(vmetrics))
}

layer_list <- catalog_apply(ctg, all_layers_catalog)

raster <- lapply(layer_list, terra::unwrap)
raster <- do.call(terra::merge, raster)

plot(raster[[2]])

[Jean-Romain]

Is is what I'm supposed to see ? Indeed it looks weird.

[RCBlackburn]

Yep. Notice the areas where chunks are smaller (i.e., edges) the pixel values have a range that is greater than the same size chunks. Also, even when chunks are the same size, you can still notice slight differences from the linear patterns in the pixels (3/4 to the top, just right of center). If you change the chunk alignment, these patterns will become more or less apparent, but they never fully disappear.

[Jean-Romain]

Let me have a guess. I did not dig into the code of your vox stuff and it is late tonight. So just a quick guess.

My guess is that ind_vox relies on the maximum height (or similar) of the point cloud to generate a voxel-full point cloud. Then vox_sum_raster use that somehow and the output depends on the number of layers of voxels. But the maximum height, and thus number of layers, is not absolute. It depends on the size and location of the chunk. Use larger chunks and you will have locally different zmax and thus different results.

I predict that the values depend on the size and alignment of the chunks. I tested it and I'm correct.

[Jean-Romain]

So, my guess was correct. First let me rewrite a bit your code

all_layers_catalog = function(las)
{ 
  las = filter_poi(las, Classification != LASNOISE)
  if (lidR::is.empty(las)) return(NULL)
  las = normalize_height(las, knnidw(k = 8, p = 2)) 
  las = filter_poi(las, Z < 50 & Z >= 1.37 ) 
  if (lidR::is.empty(las)) return(NULL)
  las <- decimate_points(las, random_per_voxel(res = 1, n = 8))
  if (lidR::is.empty(las)) return(NULL)
  
  ### voxel metrics
  las_vox <- ind_vox(las, res = 2) 
  vmetrics <- pixel_metrics(las_vox,
                            func = vox_sum_raster(SVi, FRDi, PDi, PDi_above, vox_res = 2),
                            res = 10)
  
  return(vmetrics)
}

The function can be applied

raster <- catalog_map(ctg, all_layers_catalog)

Now let apply all_layers_catalog on a point cloud where I will edit a single point to change the global maximum

las = readLAS(LASfile)
res1 = all_layers_catalog(las)
las$Z[1] = 45
res2 = all_layers_catalog(las)

plot(res1)
plot(res2)

You can observe that changing the global maximum affects all the pixels including those that are very far from the global maximum point.

To make it simple, your function does not compute a metric per pixel that is a function of the local points in the pixels, but that is also a function of the global point maximum in the point cloud. Your function does not compute valid metrics and the pixel values are potentially affected by an outlier 1 km further.

When processing by chunk, you reduced the global view, get a different global maximum for each chunk and thus your values are irrelevant and not related. This is why you can see the edges.