如何用terra包对栅格数据按月份聚合计算95分位数？

我有每日气温栅格数据，想要聚合为月度95分位数，尝试terra::tapp和gdalcubes::aggregate_time后未成功，测试代码如下：

library(terra)
r <- rast( system.file("ex/elev.tif", package="terra") )
rast <- c(r,r,r,r)
time(rast) <- as.POSIXct(c("2023-01-16","2023-01-17","2023-01-18","2023-01-19"),tz="UTC")

# 正常运行的最大值聚合
rast2 <- tapp(rast, "months", max ) 

# 无法运行的95分位数聚合
rast2 <- tapp(rast, "months", fun = function(i) quantile(i, probs = 0.95, na.rm = T))

用terra::tapp解决的核心修改 #

tapp要求自定义聚合函数返回单个数值，但quantile()默认会返回带分位数名称的命名向量（比如95%作为名称），这导致tapp无法解析。只需修改函数提取纯数值：

# 修正后的95分位数聚合代码
rast2 <- tapp(rast, "months", fun = function(i) {
  # 提取quantile返回的数值部分，去掉命名属性
  quantile(i, probs = 0.95, na.rm = TRUE)[[1]]
})

运行后即可得到每个像元的月度95分位数栅格。

用gdalcubes实现的备选方案 #

如果偏好gdalcubes，可以通过定义时间立方体视图和自定义聚合表达式实现：

library(gdalcubes)

# 创建月度时间分辨率的立方体视图
cube_view <- cube_view(
  extent = list(t0 = "2023-01-01", t1 = "2023-01-31"),
  srs = crs(rast),  # 匹配原栅格的坐标系
  dx = res(rast)[1], dy = res(rast)[2],  # 匹配原栅格分辨率
  dt = "P1M",  # 按月份分组
  aggregation = "first"  # 占位，后续替换为自定义聚合
)

# 将terra栅格转为图像集合
img_col <- create_image_collection(rast, date_time = time(rast))

# 计算月度95分位数
monthly_95p <- aggregate_time(img_col, cube_view, "quantile(x, 0.95)")

内容的提问来源于stack exchange，提问作者matej