使用 sapp() 将 terra::focal 应用到 SpatRaster

如何解决使用 sapp() 将 terra::focal 应用到 SpatRaster

我最近开始使用 terra 并且必须赞扬开发人员，它使我在 R 中处理许多大型栅格时的生活变得更加轻松。但是，我偶然发现了一个小问题，我正在尝试使用 sapp 将 focal 函数应用于 SpatRater 中的每一层，因为 focal 一次只能应用于一层。

使用一个小的可重现的 RasterSpat，我可以运行以下作为所需输出的示例：

library(terra)
packageVersion("terra")
>[1] ‘1.2.10’
    
s <- rast(system.file("ex/logo.tif",package="terra"))
s <- ifel(s == 255,1,NA)

r1 <- terra::focal(s[[1]],w=3,fun = "any",na.only=TRUE)
r2 <- terra::focal(s[[2]],na.only=TRUE)
r3 <- terra::focal(s[[3]],na.only=TRUE)
r <- c(r1,r2,r3)
r

#class       : SpatRaster 
#dimensions  : 77,101,3  (nrow,ncol,nlyr)
#resolution  : 1,1  (x,y)
#extent      : 0,77  (xmin,xmax,ymin,ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
#sources     : memory  
#              memory  
#              memory  
#names       : red,green,blue 
#min values  :   0,0 
#max values  :   1,1 

当我使用与 sapp 相同的语法为上述可重现数据运行 sapply 时：

f1 <- sapp(s,fun = function(x){terra::focal(x = x,w = 3,na.only = TRUE)})

我明白了：

h(simpleError(msg,call)) 中的错误：在为函数 'rast' 选择方法时评估参数 'x' 时出错：未使用的参数 (wopt = wopt)

如果我尝试：

f <- sapp(s,terra::focal,c(w= 3,na.only = TRUE))
f

我得到以下信息：

#class       : SpatRaster 
#dimensions  : 77,ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
#source      : memory 
#names       : red,0 
#max values  :   9,9,9 

注意图层的最大值。如何调整我的代码以使 sapp 按需要工作？

我可以使用 sapply 获得想要的结果，但我认为如果我可以让它工作，sapp 对于更大的数据会更有效。

f2 <- sapply(as.list(s),function(x){terra::focal(x = x,w= 3,na.only = TRUE)})
f2 <- rast(f2)
f2
#class       : SpatRaster 
#dimensions  : 77,ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
#sources     : memory  
               memory  
               memory  
#names       : red,1 

如果有人可以提供一些帮助，将不胜感激:)

解决方法

我已经解决了我的问题。

事实证明，您需要将 ... 参数包含到发送给 fun 的函数中，如下所示：

s <- sapp(s,fun = function(x,...) {focal(x,fun = "any",w = 3)})

然后就可以了

事实证明，sapp 内部调用 sapply 并且尚未实现并行化并且无法分配多核，因此与实现的 app 相比没有太大的速度提升并行化 vs apply。

如果有任何机构感兴趣，这里有一些基准......

terra::terraOptions(todisk = TRUE)
s <- rast(system.file("ex/logo.tif",package="terra"))
s <- terra::extend(s,c(500,500))
s <- disaggregate(s,10)
s <- terra::ifel(s == 255,1,0)

a <- function() {
    s <- sapply(as.list(s),function(x){terra::focal(x = x,w= 3,fun = "any")})
    s <- rast(s)}

b <- function(){
    s <- sapp(s,w = 3)})}


race <- microbenchmark::microbenchmark(
    a(),b(),times = 5)

结果：

race
Unit: seconds
 expr      min       lq     mean   median       uq      max neval cld
  a() 63.32553 63.93399 65.03783 65.61699 65.73256 66.58011     5   a
  b() 62.88114 63.85961 64.34571 64.16861 65.22703 65.59215     5   a