如何使用sf包计算shapefile中各多边形的最大/最小距离

实现思路 #

• 先生成所有不重复的多边形两两配对，避免重复计算、排除自身配对
• 最小距离直接使用st_distance默认配置，返回两个几何的最近点距离
• 最大距离（最远点对距离）通过st_distance的which参数指定为"Greatest"实现，若需要Hausdorff距离可替换为which = "Hausdorff"
• 最终将配对信息、两类距离整合为data.frame格式输出

代码实现 #

基础版（遍历配对，内存友好） #

适合多边形数量较多的场景，内存占用更低：

library(sf)
# 加载示例shapefile
counties <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
# 演示用取前3个多边形，全量计算直接删除该行即可
counties <- counties[1:3, ]

# 生成所有不重复的两两组合索引
pair_index <- combn(nrow(counties), 2, simplify = FALSE)

# 遍历计算距离并整合为数据框
dist_result <- do.call(rbind, lapply(pair_index, function(p) {
  p1 <- counties[p[1], ]
  p2 <- counties[p[2], ]
  # 计算最小距离
  min_d <- st_distance(p1, p2, by_element = TRUE)
  # 计算最大距离
  max_d <- st_distance(p1, p2, by_element = TRUE, which = "Greatest")
  # 可根据需求自定义输出字段，比如添加多边形的名称、属性等
  data.frame(
    poly1_id = p[1],
    poly1_name = p1$NAME,
    poly2_id = p[2],
    poly2_name = p2$NAME,
    min_distance_m = as.numeric(min_d),
    max_distance_m = as.numeric(max_d)
  )
}))

# 查看结果
print(dist_result)

简洁版（全矩阵提取，代码更短） #

适合多边形数量不多的场景，代码更易读：

library(sf)
counties <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
counties <- counties[1:3, ]

# 生成最小、最大距离全矩阵
min_dist_mat <- st_distance(counties)
max_dist_mat <- st_distance(counties, which = "Greatest")

# 提取矩阵上三角（即所有不重复两两配对的距离）
upper_pos <- upper.tri(min_dist_mat)

# 整理为数据框
dist_result <- data.frame(
  poly1_id = rep(1:nrow(counties), nrow(counties))[upper_pos],
  poly2_id = rep(1:nrow(counties), each = nrow(counties))[upper_pos],
  min_distance_m = as.numeric(min_dist_mat[upper_pos]),
  max_distance_m = as.numeric(max_dist_mat[upper_pos])
)

注意事项 #

• which参数要求sf版本≥1.0.0，版本过低可运行install.packages("sf")升级后使用
• 距离单位和输入shapefile的空间参考单位一致，示例nc.shp采用NAD83投影，单位为米
• 若多边形数量超过10000，全矩阵法会占用较多内存，推荐使用基础版遍历实现

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