在R中按距离融合同属性多边形的技术求助

我是R语言空间分析领域的新手，现需解决按位置融合多边形的问题。我有一个包含约1500个多边形的shapefile，每个多边形带有属性码，属性表中共有4种属性。需要将属性码相同且相互距离不超过20米的所有多边形进行融合，数据子集的dput输出如下：

> dput(SpatialPolygons(list(Polygons(list(Polygon(shp_sub)), 1))), control="all")
new("SpatialPolygons", polygons = list(new("Polygons", Polygons = list(
    new("Polygon", labpt = c(662107.09446353489, 670702.83409774397
    ), area = 15642909.541120354, hole = FALSE, ringDir = 1L, 
        coords = structure(c(662914.16988108843, 660819.24410834303, 
        661686.14174023282, 655590.40998587455, 662663.06165560661, 
        654986.60489466891, 656554.61134749465, 662744.22298060439, 
        653763.25605930295, 654452.19243265921, 662809.59039208689, 
        660147.03013036621, 662176.09769302572, 660630.43945095874, 
        664144.72288613662, 662914.16988108843, 663103.61422573915, 
        666180.53174116358, 665542.97683514666, 665529.15833193157, 
        664023.14232105017, 663416.57999173785, 664675.82179644238, 
        663715.85020686744, 663411.2359451541, 663617.32660092064, 
        663030.70651094906, 676394.80790592649, 675262.54310305626, 
        672662.00823758543, 671630.53395956859, 663103.61422573915
        ), .Dim = c(16L, 2L)))), plotOrder = 1L, labpt = c(662107.09446353489, 
670702.83409774397), ID = "1", area = 15642909.541120354)), plotOrder = 1L, 
    bbox = structure(c(653763.25605930295, 663030.70651094906, 
    664144.72288613662, 676394.80790592649), .Dim = c(2L, 2L), .Dimnames = list(
        c("x", "y"), c("min", "max"))), proj4string = new("CRS", 
        projargs = NA_character_))

实现思路与代码 #

核心逻辑 #

• 按属性码分组，仅处理同属性的多边形
• 对每组内的多边形，通过距离计算识别出间距≤20米的连通集群
• 对每个集群内的多边形执行融合操作
• 合并所有组的处理结果

具体代码（基于主流空间分析包sf） #

步骤1：加载依赖与数据

# 安装并加载sf包
if (!require(sf)) install.packages("sf")
library(sf)

# 读取shapefile，替换为你的文件路径
shp_data <- st_read("your_shapefile.shp")

# 关键检查：距离计算要求投影单位为米，若当前是经纬度需转换（示例为UTM投影，需替换为对应EPSG码）
# shp_data <- st_transform(shp_data, crs = 32633) # 示例EPSG，根据你的数据区域调整

步骤2：分组融合处理

# 替换'attr_code'为你的属性码列名
result_list <- lapply(split(shp_data, shp_data$attr_code), function(group) {
  # 计算组内多边形的距离矩阵
  dist_matrix <- st_distance(group)
  
  # 用层次聚类（最近邻法）识别≤20米的连通集群
  hclust_result <- hclust(as.dist(dist_matrix), method = "single")
  cluster_ids <- cutree(hclust_result, h = 20)
  
  # 给组内多边形添加集群ID
  group$cluster_id <- cluster_ids
  
  # 按集群ID融合多边形
  fused_geometry <- st_cast(st_union(group, by = "cluster_id"), "POLYGON")
  
  # 保留原属性（取每个集群第一个多边形的属性，如需聚合可修改此部分）
  attr_info <- group %>% 
    group_by(cluster_id) %>% 
    slice(1) %>% 
    st_drop_geometry() %>% 
    select(-cluster_id)
  
  # 合并属性与融合后的几何
  st_sf(attr_info, geometry = fused_geometry)
})

# 合并所有组的结果
final_fused_data <- do.call(rbind, result_list)

# 保存结果到新shapefile
st_write(final_fused_data, "fused_shapefile.shp")

重要说明 #

• 投影验证：必须确保数据投影单位为米，否则st_distance计算的距离会出错（比如经纬度的单位是度，无法直接用20米阈值）
• 聚类方法：选用层次聚类的"single"方法，能精准识别所有间距≤20米的连通多边形集群
• 属性处理：代码默认保留集群第一个多边形的属性，若需聚合属性（如求和、均值），可修改attr_info的生成逻辑

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