如何用stars/terra构建含x、y、波段、时间维度的GeoTIFF时空立方体？

一、使用stars包构建x、y、band、time四维对象 #

问题核心是每个输入文件的band维度同时包含变量和时间信息，直接堆叠时因维度标签不匹配报错。解决思路是先拆分band维度为变量（band）和时间（time），再合并时序文件：

1. 逐个读取文件并提取波段名
2. 从波段名中拆分出变量名和时间标识
3. 重塑每个文件的维度，将band拆分为band和time
4. 合并所有文件到同一stars对象

示例代码：

library(stars)

# 定义时序文件路径
urls <- c("forecast_20240101.tif", "forecast_20240102.tif")

# 处理单个文件的函数
process_file <- function(file_path) {
  # 读取文件
  s <- read_stars(file_path)
  # 获取波段名（假设格式如"temp_3h", "precip_6h"）
  band_names <- dimnames(s)[[3]]
  
  # 拆分变量名和时间标识
  vars <- gsub("_\\d+h$", "", band_names)  # 提取变量名，去掉末尾的_3h/_6h
  times <- gsub("^.*_", "", band_names)    # 提取时间标识，如3h/6h
  
  # 重塑数组：将原band维度拆分为band（变量）和time（预报时次）
  # 按x,y分组，把每个像素的波段值转为[变量数×时间数]的矩阵
  s_reshaped <- st_apply(s, c("x", "y"), function(pixel_vals) {
    matrix(pixel_vals, nrow = length(unique(vars)), ncol = length(unique(times)),
           dimnames = list(band = unique(vars), time = unique(times)))
  })
  
  return(s_reshaped)
}

# 处理所有文件
processed_list <- lapply(urls, process_file)

# 合并所有处理后的对象（沿time维度，可根据实际调整合并方向）
combined_stars <- do.call(c, c(processed_list, list(along = "time")))

二、使用terra包重构为8个band + time维度，并修复波段名 #

terra读取后会将所有波段合并为一维，需要手动拆分维度，同时从原始文件元数据中恢复被修改的波段名：

1. 恢复原始波段名 #

terra默认会将波段名中的特殊字符（如空格、减号）替换为点，可通过stars读取原始名称来修复：

library(terra)

# 读取所有文件为SpatRaster
r <- rast(urls)

# 用stars读取原始波段名（避免terra的字符替换）
s <- read_stars(urls)
original_band_names <- unlist(dimnames(s)[[3]])

# 重新赋值给terra对象
names(r) <- original_band_names

2. 重构为band + time维度 #

通过拆分波段名分组，再重新组合为目标结构：

# 从波段名拆分变量和时间
band_names <- names(r)
vars <- gsub("_\\d+h$", "", band_names)
times <- unique(gsub("^.*_", "", band_names))

# 按变量分组，每个变量下包含不同时间的波段
var_groups <- split(r, vars)

# 每个变量组内的波段按时间排序，合并为多层Raster，再组合所有变量
stacked_rast <- rast(lapply(var_groups, function(var_rast) {
  # 按时间排序波段
  time_order <- match(times, gsub("^.*_", "", names(var_rast)))
  var_rast_sorted <- var_rast[order(time_order)]
  return(var_rast_sorted)
}))

# 若需要显式的四维结构，可转换为stars对象再调整维度名称
stars_from_terra <- st_as_stars(stacked_rast)
dimnames(stars_from_terra) <- list(
  x = dimnames(stars_from_terra)[[1]],
  y = dimnames(stars_from_terra)[[2]],
  band = unique(vars),
  time = times
)

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