如何将多个 Rasterstack 聚合为一个
How to aggregate multiple Rasterstacks into one
我有几个 Rasterstacks 从几个时间序列 Netcdf 文件创建的。我想将这些汇总到 mean/median 和相关的 95% 置信区间或标准差统计数据。输出将是相同维度的单个 Rasterstack，表示跨所有 Rasterstacks 的 mean/median/stdev。
我尝试使用 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
编辑：post 提供了三种解决问题的方法。大型 RasterStacks 最快的是第三种方法，它将堆栈强制转换为数组并对其执行计算。
方法一：叠加
我假设你想要 layer-wise 统计数据，即你希望你的结果是 RasterStack 三层，第一层是四个堆栈的第一层的中位数（即中位数栅格 r1、r3、r4 和 r1），第二个是四个堆栈第二层的中值（r2、r4[ 的中值=24=]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。
方法二：uberlay
对于较大的栅格，上述方法似乎会慢很多。也许我做错了......另一种方法是更直接地调用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 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))))
首先，将堆栈强制转换为矩阵并绑定到 three-dimensional 数组。
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（如果需要，可以用 stack() 强制转换为 RasterStack），例如：
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, ...
我有几个 Rasterstacks 从几个时间序列 Netcdf 文件创建的。我想将这些汇总到 mean/median 和相关的 95% 置信区间或标准差统计数据。输出将是相同维度的单个 Rasterstack，表示跨所有 Rasterstacks 的 mean/median/stdev。
我尝试使用 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
编辑：post 提供了三种解决问题的方法。大型 RasterStacks 最快的是第三种方法，它将堆栈强制转换为数组并对其执行计算。
方法一：叠加
我假设你想要 layer-wise 统计数据，即你希望你的结果是 RasterStack 三层，第一层是四个堆栈的第一层的中位数（即中位数栅格 r1、r3、r4 和 r1），第二个是四个堆栈第二层的中值（r2、r4[ 的中值=24=]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。
方法二：uberlay
对于较大的栅格，上述方法似乎会慢很多。也许我做错了......另一种方法是更直接地调用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 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))))
首先，将堆栈强制转换为矩阵并绑定到 three-dimensional 数组。
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（如果需要，可以用 stack() 强制转换为 RasterStack），例如：
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, ...