根据季节开始和结束计算每像素平均降雨量值

Question

我有一个 365 层的 RasterStack (R1998)。每层代表一天的降雨数据（即从 01/01 到 31/12）。然后，我有两个 rasterLayers，它们针对每个像素报告一年中生长季节开始的日子 (SOS) 和一年中生长季节结束的日子 (EOS)：

> R1998
class      : RasterStack 
dimensions : 291, 327, 95157, 365  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -18.25, 63.5, -35, 37.75  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
values     : ...

> SOS
class      : RasterLayer 
dimensions : 291, 327, 95157  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : -18.25, 63.5, -35, 37.75  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 1, 365  (min, max)

> EOS
class      : RasterLayer 
dimensions : 291, 327, 95157  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : -18.25, 63.5, -35, 37.75  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 1, 365  (min, max)

我不知道如何计算每个像素仅在生长季节开始和结束之间的平均降雨量值（和其他指标）。例如：

• pixel[1,1] 的 SOS=day97 和 EOS=day128：计算 R1998 第 97 天和第 128 天之间的平均降雨量。
• pixel[1,2] 的 SOS=day95 和 EOS=day131：计算 R1998 第 95 天和第 131 天之间的平均降雨量。
• pixel[1,3] 的 SOS=day81 和 EOS=day110：计算 1998 年第 81 天和第 110 天之间的平均降雨量。
• ...所有像素以此类推。

我无法在 stackApply 中真正实现此选项。我大部分时间都在尝试使用 stackSelect，但这只是提取 R1998 中 SOS/EOS 位置的像素值。我尝试将 rasterStack 转换为数据框 (rasterToPoint)，但它变得混乱。

我确信这很简单，但我就是找不到解决方案。任何帮助将不胜感激。

（这是 Liebmann 气候指标计算的一部分；Liebmann 等人，2012 - Journal of Climate）

后续问题：对于每个像素，我现在必须计算 SOS 和 EOS 之间的日降雨量减去年平均降雨量（MAP；每像素长期平均值）。例如：

• Pixel[1,1] SOS=day97 and EOS=day128：计算rain_day97 - MAP； rain_day98 – 地图； […]； rain_day128 - 地图
• Pixel[1,2] SOS=day95 and EOS=day131：计算rain_day95 - MAP； rain_day96 – 地图； […]； rain_day131 - 地图
• ...所有像素以此类推

rapp 和 app 都为每个像素返回一个值（例如，fun=mean 返回 SOS 和 EOS 之间的平均值），但这里我需要一个值 ( rain_day – MAP）在 SOS 和 EOS 之间的每一天（即文件列表）。 tapp 似乎合适，但我发现很难回收 SOS/EOS 索引。

非常感谢

Answer 1

terra 包有一个名为rapp (range-apply) 的方法来解决这个问题。你可以这样做

x <- rapp(R1998, SOS, EOS, fun=mean)

前三个参数是SpatRaster 对象。您可以使用 Raster* 对象或文件中的 terra::rast 创建这些。

这是一个最小的、独立的、可重现的示例

示例数据

library(terra)
r <- rast(ncol=9, nrow=9)
values(r) <- 1:ncell(r)
rain <- c(r, r, r, r, r, r)
rain <- rain * 1:6
rain[1:2] <- NA
sos <- eos <- rast(r)
sos[] <- 1:3
eos[]   <- 4:6

解决方案

r <- rapp(rain, sos, eos, fun="mean")

后续问题的解决方案：

MAP <- rain+10
r <- rapp(rain-MAP, sos, eos, fun="mean")

如果你不想总结，但保留sos和eos之间的所有值，并将其他设置为NA，你可以这样做

r <- rapp(rain-MAP, sos, eos, allyrs=TRUE)

但我现在在这里看到了一个错误（terra v1.1-18 中的fixed on github。在早期版本中，您需要将 1 添加到 eos

rr <- rapp(rain-MAP, sos, eos+1, allyrs=TRUE)