1

I am trying to create subsets of a point cloud based on the ScanAngleRank and compute some metrics for each subset. Then return a stack of values(rasters) for the metric. I modified the rumple index example as shown below:

met_calc = function(cluster, res)
{
  las = readLAS(cluster)
  if (is.empty(las)) return(NULL)

  outlist <- list()
  for(i in y){
    ls <- lasfilter(las, abs(ScanAngleRank) >= min(i) & abs(ScanAngleRank) < max(i))
    if (is.empty(ls)) return(NULL)
    out <- grid_metrics(ls,~mean(Z),res)
    outlist <- c(outlist,out)
  }
  #bbox   <- raster::extent(cluster)
  #outlist <- raster::crop(outlist, bbox)
  return(outlist) 
}

opt_chunk_buffer(ldr) <- 1
opt_chunk_size(ldr) <- 300
opt_select(ldr) <- "*"
opt <- list(raster_alignment = 20, automerge = TRUE)
output <- catalog_apply(ldr, met_calc, res = 20, .options = opt)

y <- list(0:10, 10:20, 20:30, 30:40, 40:50)

In the first iteration, for example, I tried to subset las for minimum and maximum scan angle of 0 and 10 degerees respectively, and so on. In this case, that means 5 rasters.

bbox and outlist are comments because they were generating errors.

I get the following warning:

The list returned by 'catalog_apply' contains unsupported objects. Merging is impossible. A list has been returned. 

Output is a complicated list of RasterLayers for small areas.

I tried using grid_metrics() to return a list of metrics as shown in the example in the documentation. Have not been able to figure that.

I understand how to work with a single las file. My confusion is with the usage of LAScatalog and related options for this purpose, efficiently.

1

outlist is a list. The output of catalog_apply() is thus a list with one element per chunk and each element is a list. The engine is not able to merge such complex output with nested lists. So automerge option fails and the output is the regular, unmerged, list + a warning to tell you that automerge did nothing.

You are trying compute mean elevation for different subsets with more or less points with large angles. Instead of returning a list of RasterLayer you could return a RasterStack or a RasterBrick.

met_calc = function(cluster, res)
{
  las = readLAS(cluster)
  if (is.empty(las)) return(NULL)

  # Creation of a layout so each raster will be the same
  # no matter the content of the point cloud
  layout <- grid_metrics(las, ~length(Z), res)
  layout[] <- NA # see issue #318

  y <- list(0:50, 5:20, 10:20, 5:8, 2:6)

  # Initialize a first layer
  i <- y[[1]]
  ls <- lasfilter(las, abs(ScanAngleRank) >= min(i) & abs(ScanAngleRank) < max(i))
  out <- grid_metrics(ls,~mean(Z), layout)

  # Loop on other layers
  for(i in y[-1])
  {
    ls <- lasfilter(las, abs(ScanAngleRank) >= min(i) & abs(ScanAngleRank) < max(i))

    if (!is.empty(ls))
      r <- grid_metrics(ls, ~mean(Z), layout)
    else # if empty we have to return something anyway
      r <- layout

    out <- raster::addLayer(out, r)
  }

  # crop the buffer
  bbox   <- raster::extent(cluster)
  out <- raster::crop(out, bbox)
  names(out) <- c("0:50", "5:20", "10:20", "5:8", "2:6")
  return(out) 
}

library(lidR)

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

opt_chunk_buffer(ldr) <- 1
opt_chunk_size(ldr)   <- 200
opt_chunk_alignment(ldr) <- c(60,60)
opt_select(ldr)       <- "*"
opt <- list(raster_alignment = 10, automerge = TRUE)
output  <- catalog_apply(ldr, met_calc, res = 10, .options = opt)

plot(output)

enter image description here

5
  • Thankyou! Processing using LAScatalog is something I haven't been able to grasp very well yet! I was ready to subset and execute manually! Would you say this (the answer) is the most efficient way of creating such subsets? It would be nice to be able to subset the data in a custom function for grid_metrics(), for example. – K_D Jan 29 '20 at 17:42
  • I acknowledge that it is not straightforward. We are currently writing an extensive vignette that document in depth the engine. I hope it will help. For feature request ask on github or by email. – JRR Jan 29 '20 at 17:46
  • Update: It worked well for a small area which I was testing upon to save time. On applying it to a larger area (more las tiles), the function returns the same complicated list structure. Everything is the same except for the data. Was wondering if it is a random error? Didnt think a new question would be very helpful as everything is the same. – K_D Feb 3 '20 at 8:58
  • Followup: So what was happening with the larger area is that there was an error in reading some chunk and the program stopped working. I used the 'opt_stop_early(ldr) <- FALSE' option. Didn't understand why there was an error in the first place. Hope this information is helpful. I guess this is a difference I didn't mention in the previous comment. – K_D Feb 3 '20 at 23:59
  • You can report bug on github or by email if you think you found a bug – JRR Feb 4 '20 at 0:57

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.