基于另一个栅格堆栈的栅格堆栈中的总像素值

Question

我有一个代表396层蒸散发（ET）的栅格堆栈（共3个栅格层，共11年-2009年至2019年）。对于每个月，栅格图层始终代表该月的第1天，第11天和第21天，称为dekads。这是样本数据集

library(raster)

#create a raster with random numbers
r <- raster(ncol=5, nrow=5, xmx=-80, xmn=-150, ymn=20, ymx=60)
values(r) <- runif(ncell(r))

#create a random raster stack for 3 raster a month for 11 years
n <- 396 #number of raster
s <- stack(replicate(n, r)) # convert to raster stack

#rename raster layers to reflect date
d =rep(c(1,11,21),132)
m =rep(1:12, 11, each =3)
y = rep (2009:2019, each =36)
df.date <- as.Date(paste(y, m, d,sep="-"), "%Y-%m-%d")
names(s) = df.date

我还有另外两个栅格堆栈，其像素值表示2009年至2019年的季节开始（ss）11层和季节结束（se）11层。

#create a raster stack representing season start (ss) and season end (se)
# The pixel value represents dekad number. Each raster layer covers exactly three calendar years with the target year in the middle.
# (1-36 for the first year, 37-72 for the target year, 73-108 for the next year). 
ss.1 = r # season start raster
values(ss.1)= as.integer(runif(ncell(ss.1), min=1, max=72))
se.1 = ss.1+10 # season end raster
yr = 11
ss <- stack(replicate(yr, ss.1)) # season start raster stack
se <- stack(replicate(yr, se.1)) #season end rasterstack

现在，我需要从“ s”栅格堆栈中估算出每年的季节总和，以使每个像素求和的时间段应通过考虑3年的移动窗口来对应于“ ss”和“ se”的像素值。

这里是我需要一个时间步长（3年窗口），一个季节开始（ss）栅格和一个季节结束（se）栅格的输出示例。但是真正让我们感到震惊的是遍历三个栅格堆栈（s代表数据集，ss代表季节开始日期，se代表季节结束日期）。感谢您的帮助。

# Example to calculate pixel based sum for 1 time step
#subset first 3 years - equal to 108 dekads
s.sub = subset(s, 1:108)
# sum each grid cells of "s" raster stack using "ss.1" and "se.1" as an indicator for the three year subset.
for (i in 1:ncell(s.sub)) {
  x[i] <- sum(s[[ss.1[i]:se.1[i]]][i], na.rm = T)
}

Answer 1

您的示例的最后一部分没有用，我将其更改为此（一年）

x <- rep(NA, ncell(s))
for (i in 1:ncell(s)) {
  x[i] <- sum(s[i][ss.1[i]:se.1[i]], na.rm = T)
}
x <- setValues(ss.1, x)
x
#class      : RasterLayer 
#dimensions : 5, 5, 25  (nrow, ncol, ncell)
#resolution : 14, 8  (x, y)
#extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : 0.6505058, 10.69957  (min, max)

您可以这样获得结果

idx <- stack(ss.1, se.1)
thefun <- function(x, y){
    apply(cbind(y, x), 1, function(i) sum(i[(i[1]:i[2])+2], na.rm = T))
}  
z <- overlay(s, idx, fun=thefun)

有一个类似问题的more examples here。

鉴于这是一个普遍的问题，我在rapp（terra的替代品）中为其添加了一个函数raster（适用范围）--- available here;应该在7月初发布在CRAN上。

library(terra)
r <- rast(ncol=5, nrow=5, xmx=-80, xmn=-150, ymn=20, ymx=60)
values(r) <- 1:ncell(r)
s <- rast(replicate(396, r))

ss.1 <- r 
values(ss.1) <- as.integer(runif(ncell(ss.1), min=1, max=72))
se.1 <- ss.1+10 

x <- rapp(s, c(ss.1, se.1), sum)