基于sf包的多边形邻近分析、距离计算与合并实现请求

需求说明 #

• 识别每个多边形的最近多边形及两者间距离
• 获取两个多边形最近部分的坐标（用于绘制验证线段）
• 将距离≤800米的多边形合并为多部件多边形（multipart polygons）

现有代码基础 #

用户已通过st_nearest_feature完成每个多边形最近UID的匹配，在此基础上实现剩余需求：

library(sf)
library(dplyr)
library(igraph) # 用于处理连通分量合并

# 加载数据集
download.file("https://drive.google.com/uc?export=download&id=1-I4F2NYvFWkNqy7ASFNxnyrwr_wT0lGF" , destfile="ProximityAreas.zip")
unzip("ProximityAreas.zip")
Proximity_Areas <- st_read("Proximity_Areas.gpkg") 

# 匹配每个多边形的最近UID
Nearest_UID <- st_nearest_feature(Proximity_Areas)
Proximity_Areas <- Proximity_Areas %>% 
  select(UID) %>% 
  mutate(NearUID = UID[Nearest_UID])

实现两个目标输出 #

输出1：附加距离与最近点坐标的数据集 #

利用st_nearest_points获取每对多边形的最近点对，计算距离并提取坐标：

# 获取每对多边形的最近点线段
nearest_points <- st_nearest_points(Proximity_Areas, Proximity_Areas[Nearest_UID, ], pairwise = TRUE)

# 计算距离并提取最近点坐标
Proximity_Areas <- Proximity_Areas %>%
  mutate(
    distance_m = st_length(nearest_points),
    # 当前多边形的最近点坐标
    start_x = st_coordinates(st_cast(nearest_points, "POINT"))[, 1],
    start_y = st_coordinates(st_cast(nearest_points, "POINT"))[, 2],
    # 最近多边形的最近点坐标
    end_x = st_coordinates(st_cast(nearest_points, "POINT"))[, 3],
    end_y = st_coordinates(st_cast(nearest_points, "POINT"))[, 4]
  )

# 查看结果示例
head(Proximity_Areas)

注意：需确保数据集使用投影坐标系（如UTM），否则st_length返回的距离单位为度，需转换为米。

输出2：合并距离≤800米的多边形 #

需处理双向邻近关系，将所有互相邻近（距离≤800米）的多边形合并为同一多部件多边形：

# 1. 筛选有效配对并去重（避免A-B、B-A重复）
valid_pairs <- Proximity_Areas %>%
  st_drop_geometry() %>%
  filter(distance_m <= 800) %>%
  mutate(
    min_uid = pmin(UID, NearUID),
    max_uid = pmax(UID, NearUID)
  ) %>%
  distinct(min_uid, max_uid, .keep_all = TRUE) %>%
  select(UID, NearUID)

# 2. 构建连通图，识别所有连通分量
connect_graph <- graph_from_data_frame(valid_pairs, directed = FALSE)
component_mapping <- components(connect_graph)
Proximity_Areas <- Proximity_Areas %>%
  mutate(component_id = component_mapping$membership[as.character(UID)])

# 3. 按连通分量合并多边形
merged_polygons <- Proximity_Areas %>%
  group_by(component_id) %>%
  summarise(geometry = st_union(geometry), .groups = "drop") %>%
  # 保持与原数据结构一致，生成合并后的UID
  mutate(UID = paste0("merged_", component_id)) %>%
  select(UID, geometry)

# 可视化合并结果
plot(merged_polygons["geometry"])

说明：st_union会将同一分量内的多边形合并为单个multipart polygon，可根据原数据结构调整字段保留逻辑。

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