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

Posted 2023-04-18

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

【英文标题】：r raster brick sum values in the cells determined by two different rasters, how to speed up calculations

【发布时间】：2020-05-03 17:04:58

【问题描述】：

我正在处理包含每日数据的气候数据文件，因此在大多数年份中，我们将 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))
  

【问题讨论】：

您是否尝试过使用*apply-function？

我不太熟悉这些功能。我确实对堆栈溢出进行了一些研究，但没有发现任何与我想要做的相似的东西。

这段代码使用 stackApply 并且快了大约 3 倍 - for (i in 1:nrow(raster.in1)) for (j in 1:ncol(raster.in1)) datasum.sub1 <- stackApply(rasterb[[start:end]], indices, fun = sum)

【参考方案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 

【讨论】：

更简洁一点是轻描淡写。我不确定这是否符合我的要求。我需要以下内容。对于每个 1/2 度像素，获取 r.start 和 r.end 值。使用这些来确定在 b 中求和的日期数。我真正的 b 有 365*10 层（加上闰年更多）。我想创建一个新的栅格，对同一像素的 b 值求和。然后除以 10 得到 10 年期间的平均值。我意识到这个细节不在原始示例中。我试图估计的真实世界解释是作物日历中的生长期。

我相信这就是它的作用。除了用“平均值”替换“总和”——我扩展了我的答案，试图展示如何到达那里（或其他地方）

在研究您的补充（这非常有帮助！）时，我尝试使用我的数据运行您的代码。这是堆栈代码 - s <- stack(croppingCalendar_plantstart, croppingCalendar_plantend, gdd);以及来自x <- calc(s, fun = function(x) sum(x[(x[1]:x[2]) + 2])) Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function 的结果知道有什么问题吗？

我敢打赌这是因为 NA 值。查看解决方案的更新答案

还有很多 NA。栅格仅在有土地的地方才有值。我可以以某种方式掩盖并让该功能仅在带土地的像素上运行吗？同时我会试试你的新功能！