基于terra/stars空间包对NetCDF栅格逐格点逐月线性去趋势

需求说明 #

处理一个包含360个时间步（30年×12个月）的NetCDF栅格，需对每个网格单元、每个月份分别执行线性去趋势处理，最终得到与输入维度一致的残差栅格。

原有tidyverse实现（长格式数据框） #

针对[x,y,value,month,year]结构的长格式数据，可通过以下tidyverse代码实现：

library(tidyverse)

data %>%
  group_by(x, y, month) %>%
  nest() %>%
  mutate(model_trend = map(data,  ~ lm(value ~ year, data = .))) %>%
  mutate(augment_trend  = map(model_trend,  ~ broom::augment(.))) %>%
  unnest(cols = c(augment_trend)) %>%
  ungroup() %>%
  dplyr::select(x, y, year, month, .resid)

基于terra包的高效实现 #

terra针对栅格操作做了底层优化，适合处理大规模空间时序数据，内存效率与运行速度更优：

library(terra)

# 读取NetCDF栅格（替换为你的文件路径和变量名）
r <- rast("input_data.nc", var = "value")

# 提取时间维度的年份与月份信息
time_vec <- time(r)
year_vec <- as.integer(format(time_vec, "%Y"))
month_vec <- as.integer(format(time_vec, "%m"))

# 按月份分组，获取各月份对应的图层索引
month_layer_idx <- split(seq_along(r), month_vec)

# 初始化与输入维度一致的残差栅格
resid_rast <- rast(r)

# 遍历每个月份，执行分月去趋势
for (m in unique(month_vec)) {
  # 提取当前月份的所有图层
  m_layers <- r[[month_layer_idx[[as.character(m)]]]]
  # 匹配当前月份对应的年份
  m_years <- year_vec[month_vec == m]
  
  # 逐网格执行线性回归，计算残差
  resid_m <- app(m_layers, function(grid_vals) {
    if (all(is.na(grid_vals))) return(rep(NA, length(grid_vals)))
    trend_model <- lm(grid_vals ~ m_years)
    residuals(trend_model)
  })
  
  # 将当前月份的残差赋值回对应图层位置
  resid_rast[[month_layer_idx[[as.character(m)]]]] <- resid_m
}

# 保存残差栅格（可选）
writeRaster(resid_rast, "residuals_terra.nc", overwrite = TRUE)

基于stars包的实现 #

stars支持tidyverse风格的管道操作，代码可读性强，适合熟悉dplyr语法的用户：

library(stars)
library(dplyr)

# 读取NetCDF数据
s <- read_stars("input_data.nc")

# 转换格式并执行分月去趋势，最终转回stars对象
resid_stars <- s %>%
  as_tibble() %>%
  rename(value = everything()) %>%
  mutate(
    time = as.POSIXct(time),
    year = as.integer(format(time, "%Y")),
    month = as.integer(format(time, "%m"))
  ) %>%
  group_by(x, y, month) %>%
  mutate(resid = value - predict(lm(value ~ year, data = cur_data()))) %>%
  select(x, y, time, resid) %>%
  st_as_stars()

# 保存残差栅格（可选）
write_stars(resid_stars, "residuals_stars.nc", overwrite = TRUE)

方案对比 #

• terra：底层优化的栅格操作，处理大规模数据集（如千级以上网格）时速度更快，内存占用更合理，适合高性能需求场景。
• stars：语法贴近tidyverse，代码逻辑直观，无需手动管理图层索引，适合追求代码可读性与开发效率的场景。

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