使用其他栅格作为指标对栅格堆栈的每个网格单元求和的函数

## input raster
s <- stack(list.files("~/dailyraster", full.names=TRUE)) # daily raster stack
r_start <- raster("~/stackSumSTART.asc") # this raster contain starting Julian day
r_end <- raster("~/stackSumEND.asc") # this raster contain ending Julian day
noNAcells <- which(!is.na(r[])) # cell numbers which contain values

## dummy raster
x <- r
x[] <- NA 

## loop
for (i in noNAcells) {      
  x[i] <- sum(s[[r_start[i]:r_end[i]]][i])
}

我想创建一个像这样的函数stackApply()，但我希望它能够在单元格的基础上工作。

上面是一个for()循环版本，效果很好，但需要太多时间。

关键是每个单元格的范围sum()从两个栅格层，r_start, r_end在上面的脚本中。

现在我正在努力使用以下代码来转换这段代码apply() family.

有没有可能提高速度for()环形？或者请给我一些编写此代码的提示apply()

任何评论都会对我有帮助，谢谢。

你的方法

x <- s$layer.1

system.time(
for (i in 1:ncell(x)) {
     x[i] <- sum(s[[r_start[i]:r_end[i]]][i], na.rm = T)
   }
)

   user  system elapsed 
  0.708   0.000   0.710

我的建议

您可以在堆栈末尾添加用作索引的栅格，然后使用calc高度加速该过程（~30-50x）。

s2 <- stack(s, r_start, r_end)
sum_time <- function(x) {sum(x[x[6]:x[7]], na.rm = T)}

system.time(
   output <- calc(s2, fun = sum_time)
)

   user  system elapsed 
  0.016   0.000   0.015

all.equal(x, output)

[1] TRUE

样本数据

library(raster)

# Generate rasters of random values
r1 <- r2 <- r3 <- r4 <- r5 <- r_start <- r_end <- raster(ncol=10, nrow=10)

r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
r4[] <- rnorm(ncell(r4), 1, 0.2)
r5[] <- rnorm(ncell(r5), 1, 0.2)
s <- stack(r1,r2,r3,r4,r5)

r_start[] <- sample(1:2, ncell(r_start),replace = T)
r_end[] <- sample(3:5, ncell(r_end),replace = T)