r 由两个不同栅格确定的单元格中的栅格砖总和值，如何加快计算速度

Question

我正在处理包含每日数据的气候数据文件，因此在大多数年份中，我们将 365 个栅格数据整合在一起。我想总结几天子集的文件中的值 - 比如说第 x 天到第 y 天。这可以通过 stackApply 来完成。我在下面创建了一些代码，用于生成一些栅格、创建砖块并使用 x 和 y、1 和 3 的特定值应用 stackApply。

我需要的是从两个栅格图层中获取 x 和 y。在下面的代码中，它们被称为 raster.start 和 raster.end。在第一组代码下面我有第二组可以工作但速度很慢。

library(raster)
r <- raster(nrows=100, ncols=100)
s <- stack(lapply(1:5, function(i) setValues(r, runif(ncell(r), min = -10*i, max = 10))))
raster.start <- setValues(r, sample(2, ncell(r), replace=TRUE))
raster.end <- raster.start + 3
rasterb <- brick(s)

indices <- format(as.Date(names(rasterb), format = "layer.%d"), format = "%d")
indices <- c(1,1,1,1,1)

datasum.all <- stackApply(rasterb, indices, fun = sum)
datasum.sub1 <- stackApply(rasterb[[c(1:3)]], indices, fun = sum)

这个想法是通过开始和结束栅格的行和列来对砖进行子集化并对其进行操作。这是我为此开发的代码。

raster.out <- r
for (i in 1:nrow(r)){
  for (j in 1:ncol(r)){
    start <- raster.start[[1]][i,j] # get the starting day
    end <- raster.end[[1]][i,j] # get the ending day
    raster.out[i,j] <- sum(rasterb[[start:end]][i,j])
  }
}

但是，即使对于这个玩具示例，计算时间也很慢。完成大约需要 1.3 分钟。我尝试用函数替换一些代码，如下所示，但它对完成时间没有影响。非常感谢任何有关如何加快此过程的建议。

startEnd <- function(raster.start, raster.end, i,j) {
  start <- raster.start[i,j] # get the starting day
  end <- raster.end[i,j] # get the ending day
  return(c(start,end))
}

rasterOutValue <- function(rasterb, i, j, startEnd){
  return(sum(rasterb[[startEnd]][i,j]))
}

for (i in 1:nrow(raster.in1)){
  for (j in 1:ncol(raster.in1)){
    raster.out[i,j] <-rasterOutValue(rasterb, i, j, startEnd(raster.start, raster.end, i,j))
  }
}

Answer 1

您的示例数据

library(raster)
r <- raster(nrows=100, ncols=100)
set.seed(88)
b <- stack(lapply(1:5, function(i) setValues(r, runif(ncell(r), min = -10*i, max = 10))))
r.start <- setValues(r, sample(2, ncell(r), replace=TRUE))
r.end <- raster.start + 3

首先是您的示例的改进版本，该版本有效，但速度太慢。下面的速度要快得多，但仍然相当慢。

raster.out <- r
for (i in 1:ncell(r)){
    start <- raster.start[i] # get the starting day
    end <- raster.end[i] # get the ending day
    raster.out[i] <- sum(rasterb[i][start:end])
}

这将我的时间从 74 秒缩短到 5 秒。但是你永远不应该遍历单元格，这总是太慢了。相反，你可以这样做（对我来说在 0.04 秒内）：

s <- stack(r.start, r.end, b)
x <- calc(s, fun=function(x) sum(x[(x[1]:x[2])+2]))
#class      : RasterLayer 
#dimensions : 100, 100, 10000  (nrow, ncol, ncell)
#resolution : 3.6, 1.8  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : -129.5758, 30.31813  (min, max)

这似乎是正确的

a <- s[1]
a
#     layer.1.1 layer.2.1 layer.1.2 layer.2.2  layer.3   layer.4   layer.5
#[1,]         1         4 -1.789974  2.640807 4.431439 -23.09203 -5.688119    

fun <- function(x) sum(x[(x[1]:x[2])+2])
fun(a)
#[1] -17.80976
x[1]
#[1] -17.80976

calc 是光栅对象，apply 是矩阵。 （这就是为什么它在terra 中被称为app。

开始的地方是首先编写一个函数，用一个向量来做你想做的事情。

x <- 1:10
test1 <- function(start, end, values) {
    mean(values[start:end]) 
}
test1(2, 5, x)
test1(5, 8, x)

calc 只接受一个参数，所以像这样的函数

test2 <- function(values) {
    # the +2 to skip the first two elements in the computation
    start <- values[1] + 2
    end <- values[2] + 2
    mean(values[start:end]) 
}

test2(c(2, 5, x))
test2(c(5, 8, x))

还有一个更简洁的版本

test3 <- function(v) {
    mean(v[ (v[1]:v[2])+2 ] ) 
}
 test3(c(2, 5, x))
 #[1] 3.5
 test3(c(5, 8, x))
 #[1] 6.5

第二次添加（并提醒您始终检查 NA 值！）。当索引之一（开始和结束）为NA 时，test3 中断（如果其他索引为NA，则可以）

test3(c(NA, 5, x))
#Error in v[1]:v[2] : NA/NaN argument

所以我们需要一个函数来捕捉这些

test4 <- function(v) {
    if (any(is.na(v[1:2]))) {
        NA
    } else {
        mean(v[ (v[1]:v[2])+2 ] ) 
    }
}

test4(c(NA, 5, x))
#[1] NA
test4(c(1, 5, x))
#[1] 3

通常“开始”和“结束”都将同时为 NA，因此也可以使用更简单的版本

test5 <- function(v) {
    if (is.na(v[1])) {
        NA
    } else {
        mean(v[ (v[1]:v[2])+2 ] ) 
    }
}

这种calc 的方法可能会很慢，因为它将 RasterBrick 变成了具有 365 + 2 层的 RasterStack。这会显着减慢读取数据的速度。因此，您可以改用overlay 尝试这种方法（这里再次使用sum）

f <- function(i, v) {
    j <- !is.na(i[,1])
    r <- rep(NA, nrow(i))
    x <- cbind(i[j,,drop=FALSE], v[j,,drop=FALSE])
    r[j] <- apply(x, 1, function(y) sum(y[ (y[1]:y[2])+2 ] )) 
    r
}
cal <-stack(r.start, r.end)
x <- overlay(cal, b, fun= f, recycle=FALSE)
x
#class      : RasterLayer 
# ...
#values     : -129.5758, 30.31813  (min, max)

您可以通过用 Rcpp/C++ 编写算法来加快算法速度

library(Rcpp)
cppFunction('std::vector<double> gtemp(NumericMatrix cal, NumericMatrix wth) {
    std::vector<double> out(cal.nrow(), NAN);
    for (int i=0; i<cal.nrow(); i++) {
      if (!std::isnan(cal(i,0))){
         NumericVector v = wth(i,_);
         size_t start = cal(i,0)-1;
         size_t end = cal(i,1);
         out[i] = std::accumulate(v.begin()+start, v.begin()+end, 0.0);
      }  
    }
    return out;
}')

x <- overlay(cal, b, fun=gtemp, recycle=FALSE)

---

下面是使用terra（版本>= 0.6-14）和rapp（范围应用）方法的方法。

示例数据

library(terra)
d <- rast(nrows=100, ncols=100, nl=5)
rstart <- rast(d, nlyr=1)
nc <- ncell(d) 
set.seed(88)
values(d) <- t(sapply(1:5, function(i) runif(nc, min = -10*i, max = 10)))
values(rstart) <- sample(2, nc, replace=TRUE)
rend <- rstart + 3

解决方案

idx <- c(rstart, rend)
z <- rapp(d, idx, "sum")
z  
#class       : SpatRaster 
#dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
#resolution  : 3.6, 1.8  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#data source : memory 
#names       :      lyr1 
#min values  : -184.6918 
#max values  :  34.93876