基于R语言Terra计算多边形内外栅格相对相似度的技术咨询

本问题基于Stack Overflow上的SpatRaster唯一性计算方案，目标是计算研究区内每个栅格单元与一组多边形内栅格的相似度，同时考虑唯一性（定义为栅格单元与整个研究区的平均相似度）。以下是可正常运行的示例代码，现针对代码逻辑严谨性及相异度/适用性疑问进行解答：

一、示例数据准备 #

library(terra)

#### 构建示例数据 ####

# 创建多层示例栅格
set.seed(1)
r <- rast(ncol=100, nrow=100, xmin=-150, xmax=-80, ymin=20, ymax=60,
          nlyr=5, vals=runif(10000*5))
# 为部分栅格单元设置NA值
r[[1:3]][10:20] <- NA
r[100:120] <- NA

# 创建第一个多边形
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols1 <- vect(lonlat, type="polygons", crs="+proj=longlat +datum=WGS84")

# 创建第二个多边形
lon <- c(-140.8, -135.2, -134.9, -138.9, -137.2, -138.4, -140.0)
lat <- c(51.3, 52.9, 52.4, 49.8, 47.6, 48.3, 47.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols2 <- vect(lonlat, type="polygons", crs="+proj=longlat +datum=WGS84")

# 将多边形合并为单个SpatVector
pa <- vect(c(pols1, pols2))

二、分析代码 #

#### 分析流程 ####

# 用任意截距值初始化GDM模型结构（实际应使用训练好的模型参数）
gdmRastMod <- list(intercept = 2)

#### 计算研究区全局相异度 ####

# 计算栅格各层的均值（忽略NA值）
m <- global(r, "mean", na.rm=T)

# 计算栅格单元与全局均值的差值之和
edist <- sum(r - unlist(m))

# 生成全局相异度栅格（基于GDM形式的转换）
dissim <- 1 - exp(-1 * (gdmRastMod$intercept + edist))

# 替代方案：生成全局相似度栅格（结果无实际意义）
# sim <- exp(-1 * (gdmRastMod$intercept + edist))

#### 计算多边形区域的相异度 ####

# 裁剪并掩膜栅格至多边形范围
pa_rast <- crop(r, pa, mask = TRUE)

# 计算多边形范围内栅格各层的均值（忽略NA值）
m_pa <- global(pa_rast, "mean", na.rm=T)

# 计算栅格单元与多边形均值的差值之和
edist_pa <- sum(r - unlist(m_pa))

# 生成与多边形的相异度栅格
dissim_pa <- 1 - exp(-1 * (gdmRastMod$intercept + edist_pa))

#### 计算相对相似度（基于相异度比值） ####
rel_sim <- dissim_pa / dissim

三、代码逻辑严谨性反馈 #

1. 距离计算的合理性问题
   代码中用sum(r - unlist(m))直接对各层差值求和，未考虑栅格层的量纲差异。若不同层的数值范围差异较大，数值范围大的层会主导最终结果，导致距离计算的权重失衡。建议采用标准化后的距离（如欧氏距离、曼哈顿距离），或对各层进行归一化处理后再求和。

2. NA值处理的局限性
   当前逻辑中，只要栅格单元某一层存在NA，该单元的差值求和结果就会为NA。若业务允许在部分层缺失值的情况下计算相似度，需调整NA处理逻辑（如仅对有效层求和，或用均值填充NA后计算）。

3. GDM模型参数的随意性
   代码中使用了任意指定的intercept=2，但实际场景中应使用训练完成的GDM模型参数，否则相异度/相似度的计算结果仅为形式上的数值，不具备实际的模型解释性，逻辑支撑不足。

4. 相对相似度的定义匹配
   rel_sim <- dissim_pa / dissim的逻辑是“单元与多边形的相异度占其与全局均值相异度的比例”，比值越小代表单元与多边形的相似度相对于全局越高。需确保该逻辑与业务中“相对相似度”的定义完全匹配，避免概念偏差。

四、为何流程仅适用于相异度而非相似度？ #

核心原因在于相对值计算的逻辑合理性与数值特性：

1. 数值范围与解释性
   相异度dissim的取值范围是(0,1]（由1 - exp(-x)转换而来，x≥0），其比值dissim_pa / dissim的意义清晰：比值<1时，单元比全局平均更接近多边形；比值>1时则更远离。而相似度sim的取值范围同样是(0,1]，但当sim（单元与全局均值的相似度）趋近于0时，sim_pa / sim的比值会趋近于无穷大，出现极端无意义的数值，无法进行有效解释。

2. GDM模型的原生逻辑
   GDM（广义相异度模型）的核心输出是相异度，通过1 - exp(...)的转换符合模型的设计逻辑。直接用exp(...)得到的相似度，其数值的相对关系在做比值运算时会出现非线性失真，导致结果失去业务价值。

3. 代码笔误的潜在影响
   你提供的相似度代码注释中存在笔误（mod应为gdmRastMod$intercept），但即使修正笔误，上述数值特性与逻辑问题依然存在，这是导致相似度结果无意义的根本原因。

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