如何用`terra`高效组装分层栅格数据立方体，实现灵活维度子集化分析？

1. 批量加载+自动关联维度元数据 #

不用手动解析文件夹存数据框，terra的rast()支持递归通配符路径加载，直接从文件路径/名称提取模拟、变量、时间跨度等维度信息：

假设你的文件夹结构为：

simulations/
├─ sim1/
│  ├─ temp/
│  │  ├─ temp_2000-2010.tif
│  │  └─ temp_2011-2020.tif
│  └─ precip/
│     ├─ precip_2000-2010.tif
│     └─ precip_2011-2020.tif
└─ sim2/
   ├─ temp/
│  │  └─ ...
│  └─ precip/
│     └─ ...

执行以下代码批量加载并绑定元数据：

# 递归加载所有tif文件（Windows系统将"/"替换为"\\"）
r <- rast("simulations/**/*.tif")

# 从文件路径提取维度信息
paths <- sources(r)
path_parts <- strsplit(paths, "/")

# 拆分出模拟、变量、时间跨度
sim <- sapply(path_parts, function(x) x[2])
var <- sapply(path_parts, function(x) x[3])
time_period <- gsub("\\.tif$", "", sapply(path_parts, function(x) x[4]))

# 将维度信息存入栅格元数据
terra::meta(r, "simulation") <- sim
terra::meta(r, "variable") <- var
terra::meta(r, "time_period") <- time_period

# 给图层命名（方便后续识别）
names(r) <- paste(sim, var, time_period, sep = "_")

2. 灵活子集化的高效方式 #

方式一：基于元数据直接筛选 #

# 筛选sim1的所有气温数据
subset_sim1_temp <- r[terra::meta(r, "simulation") == "sim1" & terra::meta(r, "variable") == "temp"]

# 筛选2000-2010时间跨度的所有变量
subset_2000_2010 <- r[terra::meta(r, "time_period") == "2000-2010"]

方式二：通过图层名正则匹配 #

利用之前命名的图层名规则，用正则快速定位：

# 提取sim2的所有降水数据
subset_sim2_precip <- r[[grep("sim2_precip", names(r), value = TRUE)]]

方式三：按维度分组操作 #

如果需要频繁按维度批量处理，可将栅格按维度分组：

# 按「模拟+变量」分组
groups <- paste(terra::meta(r, "simulation"), terra::meta(r, "variable"), sep = "_")
r_groups <- split(r, groups)

# 直接调用分组后的栅格，比如sim1的气温数据
r_groups[["sim1_temp"]]

3. 进阶：构建时间序列数据立方体 #

如果时间跨度是连续的（如逐年数据），可将时间维度转为Date格式，让terra识别为时间序列栅格，支持时间范围筛选：

# 从文件名提取年份并转为日期（示例文件名如temp_2000.tif）
time_dates <- as.Date(paste0(gsub(".*_(\\d{4})\\.tif$", "\\1", sources(r)), "-01-01"))
terra::time(r) <- time_dates

# 按时间范围子集化
subset_2005_2010 <- r[time(r) >= as.Date("2005-01-01") & time(r) <= as.Date("2010-12-31")]

4. 性能优化要点 #

• 利用延迟加载：terra默认延迟加载栅格，仅在计算时读取数据，直接用rast()加载路径不会占用大量内存，比手动存列表高效。
• 转为NetCDF格式：将构建好的数据立方体存为NetCDF，后续可直接按维度索引子集化，性能更优：

# 保存为NetCDF
writeCDF(r, "simulation_cube.nc", varname = "value", overwrite = TRUE)

# 后续加载并多维子集化
r_cube <- rast("simulation_cube.nc")
subset_cube <- r_cube[simulation == "sim1", variable == "temp", time_period == "2000-2010"]

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