使用 TERRA 包根据插值栅格计算长期平均月降雨量？

Question

我在 R 中有一个名为 IDW3 的 SpatRaster 对象，使用 IDW 插值方法进行估计。我有 nlyr = 240，包含 12 个月 x 20 年。我需要计算各图层的长期平均月降雨量，以便最终得到 nlyr = 12，其中每一图层代表一个日历月（一月 - 十二月）。

我尝试使用下面的代码，按照这个线程 计算长期每日意味着来自 R 中的光栅，但我想验证我使用的代码。

有什么想法和意见吗？

idw3
#> class       : SpatRaster 
#> dimensions  : 723, 449, 240  (nrow, ncol, nlyr)
#> resolution  : 100, 100  (x, y)
#> extent      : 624698.7, 669598.7, 640507.8, 712807.8  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> sources     : May 1998_masked_idw3.asc  
#>               May 1999_masked_idw3.asc  
#>               May 2000_masked_idw3.asc  
#>               ... and 237 more source(s)
#> names       :     Jan 1998,     Jan 1999,     Jan 2000,     Jan 2001,     #> Jan 2002,     Jan 2003, ... 
#> min values  :           ? ,           ? ,           ? ,           ? ,           ? ,           ? , ... 
#> max values  :           ? ,           ? ,           ? ,           ? ,           ? ,           ? , ... 

## CALCULATE THE LONGTERM MONTHLY MEANS
# get the months substring
month.ltm <- substr(my, 1,3)

# calculate the ltm using tapp funtion in terra
idw3.ltm <- tapp(idw3, month.ltm, mean)
names(idw3.ltm)
#> [1] "May" "Apr" "Aug" "Jan" "Sep" "Jul" "Jun" "Feb" "Dec"
#> [10] "Nov" "Oct" "Mar"

Answer 1

您可以使用 tapp 来实现。以下是一些示例数据

library(terra)
r <- rast(ncols=10, nrows=10, nlyr=24)
values(r) <- 1:size(r)

如果数据按年排序，然后按月排序，您现在可以这样做。

x <- tapp(r, 1:12, mean)

在其他情况下，您可能必须创建另一个索引来匹配需要组合的图层。如果您的数据有时间戳，则有一些快捷方式。在这种情况下，您可以使用 index="months"

time(r) <- as.Date("2000-01-15") + 0:23 * 30.5
y <- tapp(r, "months", mean)

Answer 2

自从您发布此问题以来已经有一段时间了，但我会尝试在下一个问题中解决此问题。由于您的示例似乎无法完全重现，因此让我使用自己的数据来实现此目的。

我使用了 GPCC v2022 数据 - 过去二十年的月度数据准确地说，分辨率为 0.5° 的数据 - 来自德国气象局。

library(terra)
#> terra 1.5.21

# define filenames
files <- c("full_data_monthly_v2022_2001_2010_05.nc",
           "full_data_monthly_v2022_2011_2020_05.nc")

# create SpatRaster object
nc_data <- rast(files)

# get variable names
varnames(nc_data)
#>  [1] "precip"                       "numgauge"                    
#>  [3] "infilled_numgauges"           "interpolation_error"         
#>  [5] "interpolation_error_infilled" "diff_new_old_method"         
#>  [7] "precip"                       "numgauge"                    
#>  [9] "infilled_numgauges"           "interpolation_error"         
#> [11] "interpolation_error_infilled" "diff_new_old_method"

# subset dataset to precipitation only
nc_precip <- nc_data["precip"]

# sneak peek
nc_precip
#> class       : SpatRaster 
#> dimensions  : 360, 720, 240  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> sources     : full_data_monthly_v2022_2001_2010_05.nc:precip  (120 layers) 
#>               full_data_monthly_v2022_2011_2020_05.nc:precip  (120 layers) 
#> varnames    : precip (gpcc full data monthly product version 2022, precipitation per grid) 
#>               precip (gpcc full data monthly product version 2022, precipitation per grid) 
#> names       : precip_1, precip_2, precip_3, precip_4, precip_5, precip_6, ... 
#> unit        : mm/month, mm/month, mm/month, mm/month, mm/month, mm/month, ... 
#> time        : 2001-01-01 to 2020-12-01

正如您所看到的，这个数据集至少在信息方面与您的数据集非常相似（除了 crs、范围和分辨率）。 nlyr = 240 包含每月降水数据的 SpatRaster 对象堆栈。最显着的不同是 time 属性范围从 2001-01-01 到 2020-12-01。

但是，基本上我解决了您的问题，使用 fun = Mean 构建适当的基于时间的索引向量作为 tapp 的输入：

# get timestamps from your SpatRaster object
tst <- terra::time(nc_precip) 

# calculate monthly means
lta <- tapp(nc_precip, index = 1:12, fun = mean)

# tidy your names a little bit
names(lta) <- format(tst, "%B") |> unique()

# inspect result
lta
#> class       : SpatRaster 
#> dimensions  : 360, 720, 12  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> sources     : memory  
#>               memory  
#>               memory  
#>               ... and 9 more source(s)
#> names       :   January,  February,     March,     April,       May,      June, ... 
#> min values  :         0,         0,         0,         0,         0,         0,     ... 
#> max values  :  979.1880,  852.0020,  720.6245,  739.8225,  884.2455, 1590.6805,     ...

从我的角度来看，结果似乎是合理的，但因为这是我第一次使用 tapp，我想通过手动重新计算来确保该函数的行为符合预期：

# init an empty list for temporary storage purposes
lta <- list()

# loop monthly and calculate the long-term mean
for (i in 1:12) {
  
  idx <- seq(from = i, by = 12, length.out = n_years)
  
  lta[[i]] <- nc_precip[[idx]] |> terra::mean()
}

# create a SpatRast object with nlyr = 12
lta <- terra::rast(lta)

lta
#> class       : SpatRaster 
#> dimensions  : 360, 720, 12  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> sources     : memory  
#>               memory  
#>               memory  
#>               ... and 9 more source(s)
#> names       :      mean,      mean,      mean,      mean,      mean,      mean, ... 
#> min values  :         0,         0,         0,         0,         0,         0, ... 
#> max values  :  979.1880,  852.0020,  720.6245,  739.8225,  884.2455, 1590.6805, ...

相同的结果，唷。

编辑：

昨天发生了一些今天无法重现的奇怪行为之后，我可以确认使用 index = Months.abb 会得到与使用 index = "months 相同的结果“（正如罗伯特在下面的评论中建议的那样）：

tapp(nc_precip, index = month.abb, fun = mean)
#> class       : SpatRaster 
#> dimensions  : 360, 720, 12  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : memory 
#> names       :       Jan,       Feb,       Mar,       Apr,       May,       Jun, ... 
#> min values  :         0,         0,         0,         0,         0,         0, ... 
#> max values  :  979.1880,  852.0020,  720.6245,  739.8225,  884.2455, 1590.6805, ...

tapp(nc_precip, index = "months", fun = mean)
#> class       : SpatRaster 
#> dimensions  : 360, 720, 12  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : memory 
#> names       :        X1,        X2,        X3,        X4,        X5,        X6, ... 
#> min values  :         0,         0,         0,         0,         0,         0, ... 
#> max values  :  979.1880,  852.0020,  720.6245,  739.8225,  884.2455, 1590.6805, ...