计算栅格图层顺序列表的栅格统计数据

Question

假设我有一个 rasterstack 由 5 层组成：

library(raster)
r <- raster(nrows=10,ncols=10)
r[] <- rnorm(100)
stk <- stack(r,r*2,r+10,r-7, r*6)

我想根据不同的数字列表计算一组顺序层的统计数据。

my.seq<-as.numeric(c("2", "1", "2"))

我如何计算统计数据（例如平均值）：

层 1:2 在 s 即 my.seq[1]

s 中的第 3 层，即 my.seq[2]

层 4:5 在 s 即 my.seq[3]

我觉得最好使用循环来实现，但我不知道它是如何工作的。

如有任何帮助，我们将不胜感激。

编辑

这就是我正在尝试做的 IRL。我想我会把它包括在这里以帮助给这个问题更多的背景。

我正在尝试根据每日温度数据计算月均值。我正在使用的 rasterstack 有数百层，每个月的每一天都有一层。因此，层 1:31 放在一起时是 1970 年 1 月。在上面的示例中，my.seq 是一个列表，其中包含每个月的天数。所以我正在尝试计算1:31（一月）然后32:60（二月）等天数的平均值。

Answer 1

我不理解你的 my.seq 示例关于你在想要不同的栅格子集时所解释的内容。这是一个有效的示例，其中我使用您定义的子集创建了一个列表对象。您需要做的就是将所需的栅格子集的数字索引传递给堆栈对象的双括号索引。列表对象也必须使用双括号作为子集，所以看起来一团糟 s[[idx[[i]]]]。但是，如果将其分解，对于第一个栅格子集，它仅仅是：raster::calc(s[[1:2]], mean)

library(raster)
file <- system.file("external/test.grd", package="raster")
s <- stack(file, file, file)
s <- addLayer(s, raster(file)/2, raster(file)*2)

idx <- list( c(1:2), c(3), c(4:5) )
s.mean <- stack()
  for(i in 1:length(idx)) {
    s.mean <- addLayer(s.mean, calc(s[[idx[[i]]]], mean) )
  }

s.mean
plot(s.mean)

关于您的扩展问题，您可以在 R 中使用日期 class 创建索引，更方便的是，rts 包允许这些类型的时间序列摘要。

在这里，让我们创建一个包含 365 个图层的堆栈来表示一年中的几天。

library(raster)
f <- raster(nrows=50, ncols=50, xmn=0, xmx=10)
s <- stack()
  for(i in 1:365) {
    x <- f
    x[] <- runif(ncell(x),0,10)
    s <- addLayer(s, x)
}

现在，我们可以创建相应的日期向量。

( d <- seq(as.Date("1970/1/1"), as.Date("1970/12/31"), "days") )

可以查询日期向量以提供栅格索引。假设我们想要 12 月的平均值，我们可以使用 which 来生成索引。

( dec.idx <- which( months(d) == "December") )
( dec.mean <- calc(s[[dec.idx]], mean) )

您可以轻松创建一个列表，其中包含每个月的栅格堆栈索引，这相当于您在 my.seq 对象中描述的内容。

months.idx <- list()   
  for(m in unique( months(d) ) ) {
    months.idx[[m]] <- which( months(d) == m)  
  }  
months.idx

但是，这是在 rts 包中执行这些类型的时间摘要的内置功能，因此，我们可以通过将堆栈强制为 rts 对象然后使用其中一个应用函数来缩短任何 for 循环。

library(rts)
s.date <- rts::rts(s, d)
( s.month <- apply.monthly(s.date, mean) )

Answer 2

此函数使用 my.seq 提取感兴趣的层（包括可能重复的层，如示例中所示）。

您可以在下面提取它们的代码部分将它们放入新的 stack 或 compute/combine 层级统计数据中。

new_stack <- function(x){
  tmp <- stack()
  for(i in x){
    new <- stk[[i]]
    tmp <- stack(tmp, new)  
  }
  return(tmp)
}

tmp2 <- new_stack(my.seq)

现在，我不是 100% 确定你想要计算平均值的确切方式，所以我用以下几种方法进行了计算：

stack_mean <- function(x){
  tmp         <- stack()
  layer_means <- numeric()

  for(i in x){
    new <- stk[[i]]

    layer_mean <- mean(stk[[i]]@data@values) 
    print(paste("In layer ", i, "the layer mean is ", layer_mean))

    layer_means <- c(layer_means, layer_mean)

    tmp <- stack(tmp, new)  

  }
  print(paste("Mean of means:", mean(layer_means)))
  return(tmp)
}

stack_mean(my.seq)

这计算了几种方法。首先，层的方法（我知道不是你想要的 - 我只是为了后代展示如何）：

[1] "In layer  2 the layer mean is  -0.0981706786020096"
[1] "In layer  1 the layer mean is  -0.0490853393010048"
[1] "In layer  2 the layer mean is  -0.0981706786020096"

然后 3 层的平均值：

[1] "Mean of means: -0.081808898835008"

最后，你得到了 stack 对象的平均值，这显然是一回事，尽管我对此了解不多，也不知道这一切的真正含义：

mean(tmp2) # equivalently: mean(stack_means(my.seq))

class       : RasterLayer 
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer 
values      : -5.742679, 4.206359  (min, max)

Answer 3

你可以使用 stackApply

示例数据

library(raster)
r <- raster(nrows=10, ncols=10, vals=rnorm(100))
stk <- stack(r,r*2,r+10,r-7, r*6)
my.seq <- c(2, 1, 2)

获取索引并使用 stackApply

i <- rep(1:length(my.seq), my.seq)
x <- stackApply(stk, i, fun='mean')

计算栅格图层顺序列表的栅格统计数据

Compute raster statistics on sequential list of raster layers

r

raster

r-raster