光栅层序列表上光栅统计量的计算
假设我有一个由5个层组成的rasterstack：
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"))如何计算统计数据(例如平均数)：
s中的第1:2层，即my.seq1
s中的第三层，即my.seq2
s中的第4层:5层，即my.seq3
我觉得这最好用循环来实现，但我不知道它是如何工作的。
任何帮助都将受到感激。
编辑
这就是我想做的IRL。我想我应该把它包括在这里，以便为这个问题提供更大的背景。
我试着从每日温度数据中计算月平均值。我正在使用的栅格堆栈有数百个层，每个月的每一天都有一层。所以层1:31合在一起，1970年1月。在上面的示例中，my.seq是一个包含每个月天数的列表。所以我试着计算1:31 (1月)，然后32:60 (2月)等等的平均值。
回答 3
Stack Overflow用户
发布于 2018-06-29 01:30:15
你可以使用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')Stack Overflow用户
发布于 2018-06-28 19:05:51
我不理解你的my.seq的例子，关于你在想要不同光栅子集时所解释的内容。下面是一个有用的示例，其中我使用您定义的子集创建了一个list对象。您所需要做的就是将所需的栅格子集数字索引传递给堆栈对象的双括号索引。列表对象也必须是使用双括号的子集，因此，它看起来是一个混乱的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中的date类创建索引，更方便的是，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) )Stack Overflow用户
发布于 2018-06-28 19:07:20
该函数使用my.seq提取感兴趣的层(包括可能重复的层，如示例所示)。
您可以将它们放入一个新的stack中，或者在下面代码中提取它们的代码中的层级别上计算/组合统计数据。
new_stack <- function(x){
tmp <- stack()
for(i in x){
new <- stk[[i]]
tmp <- stack(tmp, new)
}
return(tmp)
}
tmp2 <- new_stack(my.seq) 现在，我还不能百分之百地确定你想要怎样计算平均值，所以我用了以下几种方法：
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“在第2层中，层平均为-0.0981706786020096”1“在第1层中为-0.0490853393010048”1“在第2层中为-0.0490853393010048”1“而在第2层中为-0.0981706786020096”。
然后，这三个层次的平均数：
1“均值：-0.081808898835008”
最后，您了解了stack对象的平均值，这显然是一件事情，尽管我对此不太了解，也不知道这一切到底意味着什么：
mean(tmp2) # equivalently: mean(stack_means(my.seq))
类别: RasterLayer尺寸: 10，10,100 (nrow，ncol，ncell)分辨率: 36，18 (x，y)范围：-180,180，- 90，90 (xmin，xmax，ymin，ymax)。参考：+proj=longlat +datum=WGS84数据源:内存名称:图层值：-5.742679，4.206359 (最小，最大)
https://stackoverflow.com/questions/51089335
复制相似问题
腾讯云开发者
