如何为全球栅格时间序列与单一变量构建线性回归并生成R²栅格图

我拥有一套全球栅格数据集的时间序列，以及对应年份的、所有栅格单元通用的关联变量。

修正后的示例数据生成代码 #

library(terra)
library(data.table)

# 设置栅格维度
nrow <- 360
ncol <- 720
# 时间序列长度
n <- 4

# 初始化栅格栈，直接指定范围和分辨率
raster_final <- rast(nrows = nrow, ncols = ncol, nlyr = n, 
                     ext = c(-180, 180, -90, 90), res = 0.5)
# 为每个图层赋值（模拟时间序列数据）
for(i in 1:n) {
  values(raster_final[[i]]) <- i
}

# 对应年份的关联变量
variable <- data.table(Year = 1:4, value = 0:3)

实现步骤 #

1. 定义逐像元计算R²的函数 #

针对单个栅格单元的时间序列值，拟合线性回归并提取R²：

calc_r_squared <- function(cell_values) {
  # 处理缺失值
  if (any(is.na(cell_values))) {
    return(NA)
  }
  # 拟合线性回归模型
  model <- lm(cell_values ~ variable$value)
  # 返回R²值
  return(summary(model)$r.squared)
}

2. 批量计算所有栅格单元的R² #

使用terra::app()函数对栅格栈逐像元应用上述函数：

# 生成R²栅格
r_squared_raster <- app(raster_final, calc_r_squared)

3. 可视化全球R²分布 #

直接绘制结果栅格：

plot(r_squared_raster, main = "全球栅格单元线性回归R²分布图", 
     col = hcl.colors(10, "Blues"))

注意事项 #

• 若栅格存在大量缺失值，函数会自动返回NA，不影响整体计算
• 针对大尺寸栅格或长时序数据，可通过app()的cores参数开启多线程加速（如cores = 4）
• 示例中栅格值与variable$value完全线性相关，因此所有像元R²均为1，实际数据会呈现差异化结果

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