lidR::catalog_apply()中自定义函数的下标越界错误

5cg8jx4n  于 2023-02-10  发布在  其他
关注(0)|答案(1)|浏览(66)

我正尝试使用lidR::包为大量激光雷达数据创建像素度量栅格。我想首先删除点云中的任何异常点,将点云归一化为数字地形模型,最后,在20米X 20米网格上计算标准z像素度量。我遵循了lidR::包的指导S的书和它的小插曲使用catalog_apply()引擎。我已经创建了一个“低级API”函数,它首先有一个条件来检查输入是否是LAScatalog,然后通过catalog_apply运行该函数,然后检查输入是否是LAScluster,然后直接运行该函数并从输出中剪切块缓冲区,最后检查输入是否为LAS,然后显式运行该函数。我正在努力使该函数在LAScatalog上正确运行。当我在LAS文件上运行该函数时,它可以在out error的情况下工作,但是当我在LAScatalog上运行它时,所有的块在绘图上都显示了一个错误,当例程完成时,它抛出了这个错误:

Error in any_list[[1]] : subscript out of bounds
In addition: There were 15 warnings (use warnings() to see them)

这个错误让我觉得我错过了一些catalog_apply引擎或SpatRaster驱动程序选项,这些选项告诉函数如何将输出块合并在一起以形成最终输出,但我不确定是哪个选项,而且我还没有能够在lidR:: wiki页面、小插曲或书籍上找到任何答案。我在Stackoverflow上也找不到类似的问题。2任何建议都将不胜感激。3下面是我可复制的例子:

##Loading Necessary Packages##
library(lidR)
library(future)

#Reading in the data##
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
ctg <- readLAScatalog(LASfile) # As LAScatalog
las<-readLAS(LASfile) # As LAS

####Custom function####
raster_metrics<-function(las, dtm_ras,  grid_size, sensitivity){#start function
  if(is(las, "LAScatalog")){#Start first conditional for LASCatalog
    options <-list(automerge=TRUE, need_buffer=TRUE)
    output<-catalog_apply(las, raster_metrics,  grid_size=grid_size, sensitivity=sensitivity, .options=options)
    return(output)
  } else { #end first condition start first else
    if (is(las, "LAScluster")){  #start second conditional for LAScluster
      las<-readLAS(las)
      if (is.empty(las)){return(NULL)}#Conditional for empty chunk (self contained)
      output_tmp<-raster_metrics(las, dtm_ras,  grid_size, sensitivity)
      bbox<-sf::st_bbox(las)
      output<-st_crop(output_tmp, bbox)
      return(output)
    } else {# End second conditional begin second else
      if (is(las, "LAS")){
        p95 <- pixel_metrics(las, ~quantile(Z, probs = 0.95), grid_size)
        las <- merge_spatial(las, p95, "p95")
        las <- filter_poi(las, Z < p95*sensitivity)
        las$p95 <- NULL
        norm<-las - dtm_ras
        output<-pixel_metrics(norm, .stdmetrics_z, grid_size)
        return(output)
      }else { #end final conditional begin final else
        stop("This type is not supported.")
      }#end final else
    } #end second else
  } #end first else
}  #end function

##Creating a rasterized dtm to feed to the function##
dtm_ras<-rasterize_terrain(ctg, algorithm = knnidw())

##Defining some function and engine option setttings##
grid_size<-20.0
sensitivity<-1.2
chunk_size<-grid_size*50
chunk_buffer<-grid_size*2

##Setting driver and engine option parameters##
opt_output_files(ctg)<-paste0(tempdir(), "/{XCENTER}_{YCENTER}_{ID}_Norm_Height")
opt_chunk_size(ctg)<-chunk_size
opt_chunk_buffer(ctg)<-chunk_buffer
opt_wall_to_wall(ctg)<-TRUE
opt_stop_early(ctg)<-FALSE
opt_filter(ctg)<-"-drop_withheld"
opt_select(ctg)<-"xyz"
ctg@output_options$drivers$SpatRaster$param$overwrite<-TRUE

##Setting up parallel processing##
plan(multisession, workers = nbrOfWorkers()-1)
set_lidr_threads(nbrOfWorkers()-1)

##Running the function##
example1<-raster_metrics(las=ctg, dtm_ras = dtm_ras, grid_size = grid_size, sensitivity = sensitivity)#Throws error
example2<-raster_metrics(las=las, dtm_ras = dtm_ras, grid_size = grid_size, sensitivity = sensitivity)#Works without error

2023年2月3日更新

我自己做了一些调查,发现这个错误是由内部lidR:::函数engine_merge()引发的,该函数有一个参数any_list=,这让我认为我的函数违反了catalog_apply()的模板规则之一,但我从小插图中逐字复制了模板,希望这能阐明错误的来源。

8zzbczxx

8zzbczxx1#

1.您错过了传播dtm_ras

output<-catalog_apply(las, raster_metrics,  dtm_ras = dtm_ras, grid_size=grid_size, sensitivity=sensitivity, .options=options)

1.您使用了错误的 Package 进行裁剪

bbox <-terra::ext(las)
output<-terra::crop(output_tmp, bbox)

通过以下函数,它可在sequential模式下工作

raster_metrics<-function(las, dtm_ras,  grid_size, sensitivity)
{
  if(is(las, "LAScatalog")) 
  {
    options <-list(automerge=FALSE, need_buffer=TRUE)
    output<-catalog_apply(las, raster_metrics,  dtm_ras = dtm_ras, grid_size=grid_size, sensitivity=sensitivity, .options=options)
    return(output)
  } 
  else if (is(las, "LAScluster"))
  {
    las<-readLAS(las)
    if (is.empty(las)){return(NULL)}
    output_tmp <- raster_metrics(las, dtm_ras,  grid_size, sensitivity)
    bbox <-terra::ext(las)
    output<-terra::crop(output_tmp, bbox)
    return(output)
  } 
  else if (is(las, "LAS"))
  {
    p95 <- pixel_metrics(las, ~quantile(Z, probs = 0.95), grid_size)
    las <- merge_spatial(las, p95, "p95")
    las <- filter_poi(las, Z < p95*sensitivity)
    las$p95 <- NULL
    norm <- las - dtm_ras
    output<-pixel_metrics(norm, .stdmetrics_z, grid_size)
    return(output)
  }
  else 
  {
    stop("This type is not supported.")
  }
}

然而,它不能并行工作,因为terraSpatRaster不可串行化。简单地说,当dtm_ras被发送到每个工作者时,它就不再存在。这不是lidR的问题,而是terra的问题。在lidR函数中,我使用一个内部变通方法来处理SpatRaster,将其转换为raster
在您这边,最简单的选择是使用raster中的RasterLayer

相关问题