R如何使用stackApply计算每年的季节性平均值？

Question

我想计算每年的季节性平均值，而不是整个期间的平均值。我将我的季节定义如下：DJF（12 月至 2 月）、MAM（3 月至 5 月）、JJA（6 月至 8 月）和 SON（9 月至 11 月）。

受Fredrick 问题解决方案的启发，我创建了一个代表季节的索引“组”，然后我应用了命令“stackApply”，但这个命令计算了整个期间的平均季节值。我解释说，获得的最终图层仅包含 4 个栅格，但就我而言，我想计算“每年的季节性平均值，因此每年必须有 4 个栅格，并且栅格堆栈总共应该有 136 个栅格。

在我的代码下面

感谢您的帮助

library(raster)
set.seed(123)
r <- raster(ncol=10, nrow=10)
r_brick <- brick(sapply(1:408, function(i) setValues(r, rnorm(ncell(r), i, 3))))
dim(r_brick)

dates <- seq(as.Date("1982-01-01"), as.Date("2015-12-31"), by="month")  
months <- format(dates, "%Y-%m")

groups <- function(x) {
  d <- as.POSIXlt(x)

  ans <- character(length(x))
  ans[d$mon %in%  c(11,0:1)] <- "DJF"
  ans[d$mon %in%  2:4] <- "MAM"
  ans[d$mon %in%  5:7] <- "JJA"
  ans[d$mon %in% 8:10] <- "SON"
  ans
}

data.frame(dates, groups(dates))

r_brick.s <- stackApply(r_brick, indices=groups(dates), fun=mean,na.rm=TRUE) 
nlayers(r_brick.s)

Answer 1

您的脚本一切正常，除了您正在生成的组（即 DJF，...）。您计划用于组的索引对于每个目标组应该是唯一的，而在这个阶段，1982 年的 JJA 和 1983 年的 JJA 之间没有区别，以此类推！

#data.frame(dates, paste(substr(dates, 1,4), groups(dates), sep="_"))
idx <- paste(substr(dates, 1,4), groups(dates), sep="_")

r_brick.s <- stackApply(r_brick, indices=idx, fun=mean, na.rm=TRUE) 

nlayers(r_brick.s)
#136 that is the number of seasons

Answer 2

您的示例数据

library(raster)
r <- raster(ncol=10, nrow=10)
b <- brick(sapply(1:408, function(i) setValues(r, rnorm(ncell(r), i, 3))))    
dates <- seq(as.Date("1982-01-01"), as.Date("2015-12-31"), by="month")  

正如 Majid 指出的，如果你想按年份分组，你需要使用这些

years <- as.integer(format(dates, "%Y"))
months <- as.integer(format(dates, "%m"))

现在需要对月份进行分组。请注意，当您从 12 月开始时，您必须确保不要将同一年的 1 月和 12 月组合在一起。相反，您想将 i 年的 12 月与 i+1 年的 1 月和 2 月结合起来。这是实现这一目标的一种方法（让这一年从 12 月开始！）

n <- length(months)

# move all months back one month   
mnt <- c(months[-1], ifelse(months[n] < 12, months[n]+1, 1)) 

# move the years along
yrs <- c(years[-1],  ifelse(months[n] < 12, years[n], years[n]+1)) 

# group by trimesters using integer division (or do: floor((mnt-1) / 3))
trims <- (mnt-1) %/% 3  

# get names instead of 0, 1, 2, 3
trimnms <- c("DJF", "MAM", "JJA", "SON")[trims + 1]

结合年份和名字

yt <- paste(yrs, trimnms, sep="_")

使用索引

s <- stackApply(b, indices=yt, fun=mean, na.rm=TRUE) 

如果上述将月份向后移动的业务难以遵循，请尝试仅使用几个日期（前 15 天左右）