基于terra工具合并栅格中对角相邻的多边形

问题背景 #

现有二进制栅格包含多形状斑块，通过terra包的as.polygons()转为矢量后，使用disagg()拆分得到的仅通过顶点对角相邻的多边形会被识别为独立几何对象。需要将这类对角相邻的斑块合并为同一几何对象。当前已实现循环遍历的解决方案，但因需处理3000+个各含数千多边形的矢量文件，循环效率极低，亟需更高效的实现方法。

复现代码 #

library(terra)

# 创建测试栅格
m <- matrix(
  c(1,1,0,1,1,0,1,
    1,0,0,1,1,0,0,
    0,0,0,0,0,1,0,
    0,1,0,0,0,0,0,
    0,1,1,0,1,1,1,
    0,0,0,1,0,0,0,
    0,0,0,0,0,1,1),
  nrow=7,
  ncol=7,
  byrow=TRUE
)

r <- rast(m)
plot(r)

# 转矢量并筛选目标斑块
v <- as.polygons(r)
v <- v[which(v$lyr.1 == 1),]
plot(v, "lyr.1")

# 拆分多边形并添加ID
v <- disagg(v)
v$id <- 1:nrow(v)
plot(v)
text(v, labels="id")

现有循环解决方案（效率低） #

nv <- nrow(v)

mv_df <- data.frame(id=1:nv, mval=1:nv)

for (i in 1:nv) {
  i_poly <- v[i,]
  i_ext <- extend(ext(i_poly), 1)
  i_crop <- crop(v, i_ext)
  i_touch <- is.related(i_crop, i_poly, relation="touches")
  i_update <- append(i_crop$id[i_touch], i_poly$id)
  mv_df$mval[i_update] <- min(mv_df$mval[i_update])
}

v <- merge(v, mv_df, by.x="id", by.y="id")
v <- aggregate(v, "mval")
plot(v, "mval")
text(v, labels="mval")

高效优化方案 #

利用空间邻接关系+图论连通分量的思路替代循环，可大幅提升效率：

步骤说明 #

1. 提取所有多边形的邻接关系（包含顶点接触的对角相邻）
2. 基于邻接关系构建图结构，识别连通分量（即需要合并的斑块组）
3. 按连通分量ID聚合多边形

实现代码 #

library(terra)
library(igraph)

# 复用前面创建的v对象（已拆分并添加id）
# 获取所有多边形的邻接关系矩阵
adj_matrix <- is.related(v, v, relation = "touches")

# 构建图对象
g <- graph_from_adjacency_matrix(adj_matrix, mode = "undirected")

# 识别连通分量，得到每个多边形对应的组ID
v$group_id <- components(g)$membership

# 按组ID聚合多边形
v_merged <- aggregate(v, by = "group_id")

# 可视化结果
plot(v_merged, "group_id")
text(v_merged, labels="group_id")

效率优势 #

• 避免了循环遍历每个多边形的O(n²)时间复杂度，改用图论连通分量算法（近似线性时间）
• 无需反复裁剪矢量图层，减少IO和空间计算开销
• 适合处理大规模矢量数据（上千/上万多边形）

批量处理封装 #

若需要处理批量文件，可将上述逻辑封装为函数，结合批量处理工具快速完成：

merge_diagonal_polygons <- function(vect_obj) {
  v <- disagg(vect_obj)
  adj_matrix <- is.related(v, v, relation = "touches")
  g <- graph_from_adjacency_matrix(adj_matrix, mode = "undirected")
  v$group_id <- components(g)$membership
  aggregate(v, by = "group_id")
}

# 批量处理示例（假设files是矢量文件路径列表）
# library(purrr)
# merged_list <- map(files, ~merge_diagonal_polygons(vect(.)))

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