如何在R中遍历数组并计算单像素与其余像素的时间序列相关性

准备工作：加载数据 #

首先加载ncdf4包并读取NetCDF文件：

library(ncdf4)
# 替换为你的NetCDF文件本地路径
nc_data <- nc_open("ncdd-195101-grd-scaled.nc")

读取tmax变量并查看数组结构：

tas <- ncvar_get(nc_data, "tmax")
str(tas)
# 输出示例：num [1:1385, 1:596, 1:31] NaN NaN NaN ...

tas是三维数组，维度依次为经度、纬度、时间，包含1385×596个地理像素，每个像素对应31天的时间序列。

提取目标像素的时间序列（示例为第100个经度、第245个纬度的像素）：

gg <- tas[100, 245, ]

计算所有像素与目标序列的相关性 #

方法1：向量化处理（高效推荐） #

用apply函数批量遍历前两个维度，自动提取每个像素的时间序列并计算相关性，最终输出与tas[,,1]维度一致的结果矩阵：

# 遍历每个地理像素，计算与gg的Pearson相关系数
cor_matrix <- apply(tas, MARGIN = c(1, 2), FUN = function(x) {
  # 跳过全为缺失值的像素，避免报错
  if (all(is.na(x))) return(NA)
  cor(gg, x, use = "pairwise.complete.obs")
})

# 验证结果维度，与tas[,,1]完全匹配
str(cor_matrix)

MARGIN = c(1,2)指定对数组的经度、纬度维度遍历，use = "pairwise.complete.obs"用于自动忽略两个序列中同时为缺失值的位置。

方法2：双重循环（直观但效率低） #

如果需要显式遍历逻辑，可使用嵌套循环逐个处理像素，适合小数据集测试：

# 初始化结果矩阵，维度与tas[,,1]一致
cor_matrix <- matrix(NA, nrow = nrow(tas), ncol = ncol(tas))

# 遍历所有经度
for (i in 1:nrow(tas)) {
  # 遍历所有纬度
  for (j in 1:ncol(tas)) {
    current_pixel <- tas[i, j, ]
    if (all(is.na(current_pixel))) next
    cor_matrix[i, j] <- cor(gg, current_pixel, use = "pairwise.complete.obs")
  }
}

结果说明 #

最终的cor_matrix是二维矩阵，每个位置的值对应该地理像素时间序列与目标像素gg的相关性系数，目标像素自身的相关性为1，完全匹配tas[,,1]的格式要求。

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