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

Question

我有一个代表蒸发蒸腾 (ET) 的栅格堆栈，有 396 个图层（一个月 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" 栅格堆栈中估算每年的季节性总和，以便每个像素求和的时间段应对应于 "ss" 和 [=23= 中的像素值] 通过考虑 3 年移动 window。

这是我需要一个时间步长（3 年 window）的输出示例，其中包含一个季节开始 (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个类似问题。

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

library(terra)
r <- rast(ncols=5, nrows=5, xmin=-150, xmax=-80, ymin=20, ymax=60)
values(r) <- 1:ncell(r)
s <- rast(replicate(36, 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, ss.1, se.1, sum)

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

Sum pixel values in a Raster stack based on another raster stack

r

time-series

r-raster