如何将多个 Rasterstacks 聚合为一个

Question

我有几个 Rasterstacks 从几个时间序列 Netcdf 文件创建。我想将这些汇总为均值/中位数以及相关的 95% 置信区间或标准偏差统计数据。输出将是一个相同维度的Rasterstack，代表所有Rasterstacks 的均值/中值/标准差。

我尝试使用overlay 函数，但它似乎不起作用。这是一个可重现的示例：

library(raster)
library(rgdal)
library(ncdf4)

r <- raster(ncol=10, nrow=10)
r1 <- init(r, fun=runif)
r2 <- init(r, fun=runif)
r3 <- overlay(r1, r2, fun=function(x,y){return(x+y)})
r4 <- overlay(r1, r2, fun=function(x,y){(x*y)} )
r5 <- overlay(r1, fun=sqrt)

#create rasterstacks
s1 <- stack(r1, r2,r3)
s2 <- stack(r3, r4,r5)
s3 <- stack(r4, r5, r2)
s4 <- stack(r1, r4, r3)

z<-overlay(s1, s2, s3, s4, fun=function(a,b,c,d){return(median(a,b,c,d))} )
Error in (function (x, fun, filename = "", recycle = TRUE, ...)  : 
cannot use this formula, probably because it is not vectorized

Answer 1

编辑：这篇文章提供了三种解决问题的方法。对于大型 RasterStacks 最快的是第三种方法，它将堆栈强制转换为数组并对其执行计算。

---

方法一：叠加

我假设您需要逐层统计，即您希望您的结果是具有三层的RasterStack，第一层是四叠第一层的中位数（即栅格的中位数r1， r3、r4和r1），第二个是四个堆栈第二层的中位数（r2的中位数，r4,r5, andr4`），以此类推。

您可以通过Vectorize 函数mean、median 和sd 来实现这一点：

overlay(s1, s2, s3, s4, fun=function(...) Vectorize(median, 'x')(list(...)))

## class       : RasterBrick 
## dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       :    layer.1,    layer.2,    layer.3 
## min values  : 0.01763912, 0.01018932, 0.24531431 
## max values  :  0.9933407,  0.9050321,  1.4268951

根据需要将median 替换为mean 或sd。

---

方法二：优步

对于较大的栅格，上述方法似乎会减慢很多。也许我做错了......另一种方法是更直接地调用mapply：

uberlay <- function(..., fun) {
  fun <- match.fun(fun)
  L <- lapply(list(...), unstack)
  stack(do.call(mapply, c(FUN=function(...) calc(stack(...), fun), L)))
}

将 RasterStacks 传递给 ...，将函数传递给 fun。

uberlay(s1, s2, s3, s4, fun='median')

## class       : RasterStack 
## dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names       :    layer.1,    layer.2,    layer.3 
## min values  : 0.01763912, 0.01018932, 0.24531431 
## max values  :  0.9933407,  0.9050321,  1.4268951

---

方法 3：superduperlay

@Joe mentioned uberlay 方法需要大约一个小时才能处理他的数据。对于大堆栈，将堆栈强制转换为数组（或者，例如，data.table）并对其执行计算会更快。

让我们使用@Joe 的维度创建一些假数据：

library(raster)
library(abind)

nc <- nr <- 17
nl <- 5829

s1 <- stack(replicate(nl, raster(matrix(runif(nr*nc), nr))))
s2 <- stack(replicate(nl, raster(matrix(runif(nr*nc), nr))))
s3 <- stack(replicate(nl, raster(matrix(runif(nr*nc), nr))))
s4 <- stack(replicate(nl, raster(matrix(runif(nr*nc), nr))))
s5 <- stack(replicate(nl, raster(matrix(runif(nr*nc), nr))))

首先，将堆栈强制转换为矩阵并绑定到一个三维数组。

A <- abind(as.matrix(s1), as.matrix(s2), as.matrix(s3), as.matrix(s4), as.matrix(s5), 
           along=3)

现在将您的函数应用到边距1:2，调整尺寸并转置，然后堆叠回RasterBrick：

z <- apply(A, c(1:2), median) # substitute median with desired function
dim(z) <- c(nr, nc, nl)
z <- apply(z, c(1, 3), t)
b <- brick(z)

对于median 和sd，整个过程（包括创建阵列）在我的系统上只需要 30 多秒。对于mean，您可以利用colMeans，将速度提高到 3 秒以下。为方便起见，我们可以将它们全部封装成一个函数：

superduperlay <- function(..., fun) {
  require(abind)
  require(raster)
  fun <- match.fun(fun)
  L <- list(...)
  A <- do.call(abind, c(lapply(L, as.matrix), along=3))
  if(as.character(match.call()['fun'])=='mean') {
    A <- aperm(A, c(3, 1, 2))
    z <- colMeans(A)
  } else {
    z <- apply(A, c(1:2), fun)
  }
  dim(z) <- c(nr, nc, nl)
  z <- apply(z, c(1, 3), t)
  b <- brick(z)
}

system.time(my_mean <- superduperlay(s1, s2, s3, s4, s5, fun='mean'))
##    user  system elapsed 
##    2.68    0.04    2.72 

system.time(my_median <- superduperlay(s1, s2, s3, s4, s5, fun='median'))
##    user  system elapsed 
##   31.75    0.06   31.92 

每个对象都是一个RasterBrick（如果需要，可以强制转换为RasterStack，使用stack()），例如：

my_mean

## class       : RasterBrick 
## dimensions  : 17, 17, 289, 5829  (nrow, ncol, ncell, nlayers)
## resolution  : 0.05882353, 0.05882353  (x, y)
## extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       :    layer.1,    layer.2,    layer.3,    layer.4, ... 
## min values  : 0.19478752, 0.14775996, 0.15108237, 0.14281812, ... 
## max values  :  0.8388662,  0.8577153,  0.8396123,  0.7781535, ...