基于空间链接距离计算多边形质心距离矩阵的方法优化咨询

需求背景 #

计算区域内多边形质心的空间链接距离矩阵：距离定义为沿空间邻域链接的欧氏距离之和（而非直接欧氏距离），即通过多边形邻接拓扑，计算任意两个质心之间的最短路径总长度。

原方案合理性验证 #

你提出的基于sfdep+cppRouting的实现逻辑完全合理，核心流程符合空间链接距离的计算要求：

1. 提取多边形邻接关系，计算相邻质心的欧氏距离
2. 将邻接关系转换为网络拓扑的OD（起点-终点）格式
3. 基于拓扑网络计算所有节点对的最短路径，得到空间链接距离矩阵
4. 选择cppRouting而非igraph是正确的，其底层C++实现在大规模数据场景下效率远高于igraph

优化方法与高效实现 #

针对原方案的可优化点，提供以下改进方向：

1. 简化OD格式转换，降低内存消耗 #

原方案通过生成全邻接矩阵再筛选非零值，会占用O(n²)的内存（n为多边形数量），当n较大时内存开销极高。直接从邻接列表生成OD对，仅占用O(m)内存（m为邻接边数，远小于n²）：

nb_list_to_df <- function(nb, dists) {
  from <- rep(seq_along(nb), lengths(nb))
  to <- unlist(nb)
  dist <- unlist(dists)
  data.table(from = from, to = to, dist = dist)
}

2. 自定义邻接规则与质心计算 #

• 邻接关系：st_contiguity默认使用queen邻接（共享边或顶点即视为邻接），若业务需要仅共享边的rook邻接，可修改参数：

  nb <- sfdep::st_contiguity(geo, queen = FALSE)
  

• 质心计算：默认质心可能不适用于带洞多边形，可手动计算最大面的质心：

  centroids <- sf::st_centroid(geo, of_largest_polygon = TRUE)
  dists <- sfdep::st_nb_dists(centroids, nb)
  

3. 完整优化后代码示例 #

library(sf)
library(sfdep)
library(data.table)
library(cppRouting)

# 加载示例数据
data(guerry)
geo <- sf::st_geometry(guerry)

# 自定义质心（针对带洞多边形）
centroids <- sf::st_centroid(geo, of_largest_polygon = TRUE)

# 获取邻接关系（可选择queen/rook邻接）
nb <- sfdep::st_contiguity(geo, queen = TRUE)

# 计算相邻质心的欧氏距离
dists <- sfdep::st_nb_dists(centroids, nb)

# 高效生成OD格式数据框
nb_list_to_df <- function(nb, dists) {
  from <- rep(seq_along(nb), lengths(nb))
  to <- unlist(nb)
  dist <- unlist(dists)
  data.table(from = from, to = to, dist = dist)
}
od_df <- nb_list_to_df(nb, dists)

# 构建无向拓扑网络
graph <- makegraph(od_df, directed = FALSE)

# 获取所有多边形ID，计算全节点对的空间链接距离矩阵
all_nodes <- unique(c(od_df$from, od_df$to))
dist_link <- get_distance_matrix(Graph = graph, from = all_nodes, to = all_nodes)

4. 超大规模数据的并行优化 #

当多边形数量超过1000时，可结合parallel包分块计算距离矩阵，进一步提升效率：

library(parallel)
# 按节点分块
node_chunks <- split(all_nodes, cut(seq_along(all_nodes), 4)) # 分成4块
# 并行计算每块的距离
cl <- makeCluster(4)
clusterExport(cl, c("graph", "get_distance_matrix"))
dist_list <- parLapply(cl, node_chunks, function(chunk) {
  get_distance_matrix(Graph = graph, from = chunk, to = all_nodes)
})
stopCluster(cl)
# 合并结果
dist_link <- do.call(rbind, dist_list)

总结 #

原方案逻辑严谨，完全满足需求；优化后的方案在内存效率、自定义灵活性上有显著提升，更适合处理大规模空间数据集。

内容的提问来源于stack exchange，提问作者rafa.pereira