我有一个365层的RasterStack(R1998(。每一层代表一天的降雨量数据(即从01月1日到12月31日(。然后，我有两个光栅层，为每个像素报告生长季节开始的当天(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=第97天，EOS=第128天：计算1998年第97天和第128天之间的平均降雨量
- pixel[1，2]的SOS=第95天，EOS=第131天：计算1998年第95天和第131天之间的平均降雨量
- pixel[1,3]的SOS=第81天，EOS=第110天：计算1998年第81天和第110天之间的平均降雨量
- 。。。等等
我无法在stackApply中真正实现此选项。我花了大部分时间尝试stackSelect，但这只是提取R1998中SOS/EOS位置的像素值。我试图将rasterStack转换为数据帧(rasterToPoint(，但它变得一团糟。
我相信这很简单，但我就是找不到解决方案。如有任何帮助，我们将不胜感激。
(这是利伯曼气候指标计算的一部分；利伯曼等人，2012年-气候杂志(
后续问题：对于每个像素，我现在必须计算SOS和EOS之间的日降雨量减去年平均降雨量(MAP；每个像素的长期平均值(。例如：
- 像素[1,1]SOS=第97天，EOS=第128天：计算rain_day97–MAP；rain_day98–MAP；[…]；rain_day128–MAP
- Pixel[1,2]SOS=第95天，EOS=第131天：计算rain_day95–MAP；rain_day96–MAP；[…]；rain_day131–MAP
- …依此类推所有像素
rapp和app都为每个像素返回一个值(例如，fun=mean返回SOS和EOS之间的平均值(，但这里我需要SOS和EOS之间的每一天都返回一个数值(rain_day–MAP((即文件列表(tapp似乎很合适，但我发现很难回收SOS/EOS索引。
非常感谢
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中修复了github。在早期版本中，您需要在eos中添加1
rr <- rapp(rain-MAP, sos, eos+1, allyrs=TRUE)