在 R 中按周聚合每日栅格数据

Question

我有来自 NOAA 的 2022 年每日冰雪范围栅格文件 (https://noaadata.apps.nsidc.org/NOAA/G02156/GIS/4km/2022/)。我需要获得一年中每周的最大冰面积。

我的任务与此处描述的类似，但我需要保留日期以供将来使用。

我尝试使用tapply，但我需要保留一周的开始和结束日期来重命名新的每周栅格，最好采用“YYYYDDDYYYYDDD_ice”格式。我需要使用这些栅格作为其他每周栅格数据的掩码，因此我需要能够匹配日期。

这是我正在处理的内容：

library(terra)
library(dplyr)

ICE_dat <- list.files(ICE_path,
                      full.names = TRUE,
                      pattern = ".tif$")

ICE_stack <- rast(ICE_dat)
print(ICE_stack)
#> class       : SpatRaster 
#> dimensions  : 6144, 6144, 365  (nrow, ncol, nlyr)
#> resolution  : 4000, 4000  (x, y)
#> extent      : -12288000, 12288000, -12288000, 12288000  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +x_0=0 +y_0=0 +a=6378137 +rf=291.505347349177 +units=m +no_defs 
#> sources     : ims2022001_4km_GIS_v1.3.tif  
#>               ims2022002_4km_GIS_v1.3.tif  
#>               ims2022003_4km_GIS_v1.3.tif  
#>               ... and 362 more source(s)
#> names       : ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ... 
#> min values  :          ? ,           0,           0,           0,           0,           0, ... 
#> max values  :          ? ,           4,           4,           4,           4,           4, ...

# Extract dates
n <- names(ICE_stack)
year <- as.numeric(substr(n, 4, 7))
doy  <- as.numeric(substr(n, 8, 10))
date <- as.Date(doy, origin=paste(year-1, "-12-31", sep=""))
week <- floor(doy / 7) + 1 #set my own weeks
newweek <- c(0, 1, week[-c(364,365)]) # turn Jan 1 to zero then add another day to week 1 and get rid of 364 and 365 values, which are week 53
myweek <- formatC(newweek, width=2, flag=0)
myweek
#>   [1] "00" "01" "01" "01" "01" "01" "01" "01" "02" "02" "02" "02" "02" "02" "02"

# Assign date to rasters
terra::time(ICE_stack) <- date

# Summarize by week
wk <- tapp(ICE_stack, myweek, max, na.rm=TRUE)

print(wk)
#> class       : SpatRaster 
#> dimensions  : 6144, 6144, 53  (nrow, ncol, nlyr)
#> resolution  : 4000, 4000  (x, y)
#> extent      : -12288000, 12288000, -12288000, 12288000  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +x_0=0 +y_0=0 +a=6378137 +rf=291.505347349177 +units=m +no_defs 
#> source      : spat_198457836cc_408.tif 
#> names       : X00, X01, X02, X03, X04, X05, ... 
#> min values  :   0,   0,   0,   0,   0,   0, ... 
#> max values  :   4,   4,   4,   4,   4,   4, ...
names(wk)
#>  [1] "X00" "X01" "X02" "X03" "X04" "X05" "X06" "X07" "X08" "X09" "X10" "X11"
#> [13] "X12" "X13" "X14" "X15" "X16" "X17" "X18" "X19" "X20" "X21" "X22" "X23"
#> [25] "X24" "X25" "X26" "X27" "X28" "X29" "X30" "X31" "X32" "X33" "X34" "X35"
#> [37] "X36" "X37" "X38" "X39" "X40" "X41" "X42" "X43" "X44" "X45" "X46" "X47"
#> [49] "X48" "X49" "X50" "X51" "X52"

我希望名称（周）包含日期。如果可能的话，开始日期和结束日期，或者至少是一周的结束日期。 新的“rts”或“tidyterra”包可以提供一些解决方案吗？

Answer 1

据我判断，

names(wk)

仅包含传递给

index

参数的（聚合）级别

terra::tapp()

，因此您无法在此处本地自定义间隔定义，但当然您可以自定义名称之后手动生成结果层。

顺便说一句，在当前状态下使用

myweek

的方法会在处理多年时汇总跨年的数据，因此除非这是您的意图，否则您可能应该只附加年份。

仅使用 2022 年 1 月的数据，您将得到以下结果：

library(terra)
#> terra 1.7.78
library(dplyr)

fnames <- list.files(pattern = "tif$")

r <- rast(fnames)

# set time attribute
year <- substr(fnames, 4, 7)
doy <- substr(fnames, 8, 10)

time(r) <- paste0(year, "-", doy) |> strptime("%Y-%j") |> as.Date()
r
#> class       : SpatRaster 
#> dimensions  : 6144, 6144, 31  (nrow, ncol, nlyr)
#> resolution  : 4000, 4000  (x, y)
#> extent      : -12288000, 12288000, -12288000, 12288000  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +x_0=0 +y_0=0 +a=6378137 +rf=291.505347349177 +units=m +no_defs 
#> sources     : ims2022001_4km_GIS_v1.3.tif  
#>               ims2022002_4km_GIS_v1.3.tif  
#>               ims2022003_4km_GIS_v1.3.tif  
#>               ... and 28 more source(s)
#> names       : ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ims20~_v1.3, ... 
#> time (days) : 2022-01-01 to 2022-01-31

# built index
week_levels <- time(r) |> format("%Y-%W")
week_levels
#>  [1] "2022-00" "2022-00" "2022-01" "2022-01" "2022-01" "2022-01" "2022-01"
#>  [8] "2022-01" "2022-01" "2022-02" "2022-02" "2022-02" "2022-02" "2022-02"
#> [15] "2022-02" "2022-02" "2022-03" "2022-03" "2022-03" "2022-03" "2022-03"
#> [22] "2022-03" "2022-03" "2022-04" "2022-04" "2022-04" "2022-04" "2022-04"
#> [29] "2022-04" "2022-04" "2022-05"

# aggregate by time
wk <- tapp(r, week_levels, max, na.rm = TRUE)

# not bad, but let's tidy this a bit
names(wk)
#> [1] "X2022.00" "X2022.01" "X2022.02" "X2022.03" "X2022.04" "X2022.05"

我很确定有一种更优雅的方法可以根据您的需要“修复”名称，但看起来天真的 for 循环也可以完成工作：

# get individual indices used
wl <- week_levels |> unique()

for (i in 1:length(wl)) {
  
  # extract relevant year and doy
  ind <- week_levels == wl[i]
  
  y <- year[ind] |> max()
  
  d <- doy[ind] |> range() |> stringr::str_pad(width = 3, side = "left", pad = "0")
  
  # construct your interval definition to be used as layer name 
  rname <- paste0(y, d) |> paste0(collapse = "_")
  
  # and assign this to your aggregated layer
  names(wk[[i]]) <- rname
}

names(wk)
#> [1] "2022001_2022002" "2022003_2022009" "2022010_2022016" "2022017_2022023"
#> [5] "2022024_2022030" "2022031_2022031"