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

Question

## 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()循环版本，效果不错，但耗时太长。

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

现在我正在努力使用apply() family 转换此代码。

有没有可能通过for() 循环提高速度？或者请给我一些提示以在apply()中编写此代码

任何cmets都会帮助我，谢谢。

Answer 1

你的方法

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)