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

编程问答 2022-05-15

问题描述

我有一个 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,95157  (nrow,ncell)
resolution : 0.25,ymax)
crs        : +proj=longlat +datum=wgs84 +no_defs 
source     : memory
names      : layer 
values     : 1,365  (min,max)

> EOS
class      : RasterLayer 
dimensions : 291,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：计算 R1998 年第 81 天和第 110 天之间的平均降雨量。
  • ...对所有像素依此类推。

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

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

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

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

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

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

非常感谢

解决方法

terra 包有一个名为 rapp（范围应用）的方法来解决这个问题。你可以这样做

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

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

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

示例数据

library(terra)
r <- rast(ncol=9,nrow=9)
values(r) <- 1:ncell(r)
rain <- c(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,fun="mean")

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

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

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

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