我正尝试使用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()
的模板规则之一,但我从小插图中逐字复制了模板,希望这能阐明错误的来源。
1条答案
按热度按时间8zzbczxx1#
1.您错过了传播
dtm_ras
1.您使用了错误的 Package 进行裁剪
通过以下函数,它可在
sequential
模式下工作然而,它不能并行工作,因为
terra
的SpatRaster
不可串行化。简单地说,当dtm_ras
被发送到每个工作者时,它就不再存在。这不是lidR
的问题,而是terra
的问题。在lidR
函数中,我使用一个内部变通方法来处理SpatRaster
,将其转换为raster
。在您这边,最简单的选择是使用
raster
中的RasterLayer
。