基于R语言terra包计算大GDM栅格的(不)相似性与独特性

背景与需求 #

需将Mokany等人（2022, Global Ecology and Biogeography）提出的基于raster包的GDM空间分析工作流，迁移至terra包以处理超大规模栅格（测算需102GB内存）。核心目标是基于GDM输出的转换预测因子多层栅格，完成：

• 计算网格单元间的预测不相似性，输出单图层栅格
• 计算网格单元的独特性（与研究区采样点的平均相似性），输出单图层栅格

示例数据与初始尝试 #

构建示例多层栅格 #

library(terra)
# 创建带NA值的5层示例栅格
r <- rast(ncol=100, nrow=100, xmin=-150, xmax=-80, ymin=20, ymax=60,
          nlyr=5, vals=runif(10000*5))
r[[1:3]][10:20] <- NA
r[100:120] <- NA

初始尝试代码（存在内存与效率问题） #

# 采样50%的栅格单元
n.sub <- ncell(r) * 0.5
x <- na.omit(spatSample(r, n.sub, "regular", as.df=FALSE)) 

# 初始化不相似性矩阵（大规模数据下此步骤会直接内存溢出）
sub.dissimilarity <- matrix(0, n.sub, n.sub)

# 模拟GDM模型参数
gdmRastMod <- list(intercept = 2)

# 双重循环计算不相似性（效率极低，大规模数据不可行）
for(i in 1:(n.sub-1)) {
    for(j in (i+1):n.sub) {
        ecol.dist <- sum(abs(x[i, ] - x[j, ]))
        sub.dissimilarity[j,i] <- 1 - exp(-1 * (gdmRastMod$intercept + ecol.dist))
    }
}

当前问题（更新） #

已参考方案成功生成可扩展的独特性栅格，但仍存在两个问题：

1. 计算单元间不相似性时输出为多层栅格，需生成与独特性栅格一致的单图层不相似性栅格
2. 无法将不相似性计算逻辑扩展至所有网格单元（替换采样数为总单元格数时触发内存报错）

解决方案 #

1. 独特性栅格计算（已验证可扩展） #

独特性定义为“单元与研究区采样点的平均相似性”，可通过向量化操作替代循环，结合采样近似实现高效计算：

# 采样研究区代表性点集（根据内存调整采样量，示例用1000个点）
sample_set <- na.omit(spatSample(r, 1000, "random", as.df = FALSE))

# 定义GDM相似性计算函数（相似性=1-不相似性）
gdm_similarity <- function(x) {
  if (any(is.na(x))) return(NA)
  ecol_dists <- rowSums(abs(t(sample_set) - x))
  diss_vals <- 1 - exp(-1 * (gdmRastMod$intercept + ecol_dists))
  mean(1 - diss_vals)
}

# 应用到整个栅格，生成独特性栅格
uniqueness_rast <- app(r, fun = gdm_similarity)

2. 单图层不相似性栅格计算 #

场景1：单个参考单元与所有单元的不相似性

若需计算每个单元与指定参考单元的不相似性，生成单图层栅格：

# 选取一个非NA的参考单元（示例取第一个非NA单元）
ref_cell <- na.omit(values(r))[1, ]

# 定义单对多不相似性函数
single_ref_diss <- function(x) {
  if (any(is.na(x))) return(NA)
  ecol_dist <- sum(abs(x - ref_cell))
  1 - exp(-1 * (gdmRastMod$intercept + ecol_dist))
}

# 生成单图层不相似性栅格
ref_diss_rast <- app(r, fun = single_ref_diss)

场景2：全局平均不相似性（每个单元与区域内采样点的平均不相似性）

直接计算所有单元对不可行，采用采样近似法，生成全局平均不相似性的单图层栅格：

# 复用之前的sample_set采样点集
avg_diss <- function(x) {
  if (any(is.na(x))) return(NA)
  ecol_dists <- rowSums(abs(t(sample_set) - x))
  diss_vals <- 1 - exp(-1 * (gdmRastMod$intercept + ecol_dists))
  mean(diss_vals)
}

# 生成单图层平均不相似性栅格
avg_diss_rast <- app(r, fun = avg_diss)

3. 超大规模栅格的分块处理 #

针对102GB级别的栅格，使用terra的分块读写功能避免内存溢出：

# 初始化输出栅格
avg_diss_rast <- rast(r, nlyr = 1)
names(avg_diss_rast) <- "average_dissimilarity"

# 获取分块参数
bs <- blockSize(r)

# 分块计算并写入
for (i in 1:bs$n) {
  cat("Processing block", i, "/", bs$n, "\n")
  # 读取当前块的数值
  block_vals <- values(r, row = bs$row[i], nrow = bs$nrow[i])
  # 计算当前块的平均不相似性
  block_results <- apply(block_vals, 1, avg_diss)
  # 写入结果到输出栅格
  writeValues(avg_diss_rast, block_results, bs$row[i])
}

# 保存最终结果
writeRaster(avg_diss_rast, "global_average_dissimilarity.tif", overwrite = TRUE)

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