R中如何合并地理一致但ID不匹配的两个多面图层空间数据

场景说明 #

两个存储美国人口普查区的多面shapefile图层ID编码体系不互通，无法通过属性ID做匹配合并，需要完成空间属性连接。

该操作等价于QGIS中「Data Management Tools > Join Attributes by Location」功能，最终输出的空间对象需完整保留两个图层的所有属性，同时包含people和homes字段。

以下是示例数据构建代码：

library(sp)

# 构建第一个多边形图层
square <- rbind(c(255842.4, 4111578, 255862.4, 4111578, 255862.4, 4111558, 
                  255842.4, 4111558, 255842.4, 4111578, 255842.4, 4111578),
                c(257397.0, 4111309, 257417.0, 4111309, 257417.0, 4111289, 
                  257397.0, 4111289, 257397.0, 4111309, 257397.0, 4111309))

ID <- c("C1", "C2")
people <- c(5000, 3000)

polys <- SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID[2])
))

p.df <- data.frame( ID=1:length(polys)) 
pid <- sapply(slot(polys, "polygons"), function(x) slot(x, "ID"))
p.df <- data.frame( ID=1:length(polys), row.names = pid)    
p <- SpatialPolygonsDataFrame(polys, p.df)
p@data$people <- c(7500,3000) 


# 构建第二个多边形图层
square2 <- rbind(c(255842.4, 4111578, 255862.4, 4111578, 255862.4, 4111558,
                   255842.4, 4111558, 255842.4, 4111578, 255842.4, 4111578),
                 c(257397.0, 4111309, 257417.0, 4111309, 257417.0, 4111289,
                   257397.0, 4111289, 257397.0, 4111309, 257397.0, 4111309))

ID2 <- c("30067893", "30794385")

polys2 <- SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID2[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID2[2])
))

p2.df <- data.frame( ID=1:length(polys2)) 
pid2 <- sapply(slot(polys2, "polygons"), function(x) slot(x, "ID"))
p2.df <- data.frame( ID=1:length(polys2), row.names = pid2)    
p2 <- SpatialPolygonsDataFrame(polys2, p2.df)
p2@data$homes <- c(500,250) 

实现逻辑 #

放弃基于属性ID的匹配思路，直接以多边形的空间位置、几何重叠关系为匹配依据，和QGIS按位置连接的底层逻辑完全一致，不需要两个图层存在共通的属性字段。 #

方案1：兼容sp对象的原生实现 #

如果需要继续沿用示例中SpatialPolygonsDataFrame格式的对象，用rgeos包做空间计算即可完成匹配：

library(rgeos)

# 逐要素匹配空间位置完全重合的多边形
match_index <- sapply(seq_along(p), function(i){
  # 计算当前多边形与另一图层所有多边形的重叠面积，取完全重叠的匹配项
  intersect_res <- gIntersection(p[i,], p2, byid = TRUE)
  if(is.null(intersect_res)) return(NA)
  overlap <- gArea(intersect_res, byid = TRUE)
  which.max(overlap)
})

# 绑定匹配到的属性字段
p@data$homes <- p2@data$homes[match_index]

执行完成后p就是合并完成的对象，属性表同时包含ID、people、homes三个字段，几何保持不变。 #

方案2：基于sf包的简洁实现（推荐） #

sf是R当前主流的空间数据处理工具，不需要操作S4对象的底层slot，语法更简洁，处理大文件速度更快：

library(sf)

# 格式转换（如果直接用st_read读shapefile可跳过这步，直接得到sf对象）
p_sf <- st_as_sf(p)
p2_sf <- st_as_sf(p2)

# 按几何完全相等规则做空间连接
merged <- st_join(
  x = p_sf,
  y = p2_sf[, "homes"],
  join = st_equals
)

# 需要转回sp格式的话执行以下代码即可
merged_sp <- as(merged, "Spatial")

如果实际数据存在边界数字化带来的微小偏差，可将st_equals替换为st_is_within_distance，设置合理的距离阈值即可完成容错匹配。 #

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