高效提取多多边形中自相交特征生成的所有子多边形

从包含大量（大约 20000 个）可能部分重叠的多边形的 shapefile 开始，我需要提取通过相交不同“边界”而产生的所有子多边形。

在实践中，从一些模拟数据开始：

library(tibble)
library(dplyr)
library(sf)

ncircles <- 9
rmax     <- 120
x_limits <- c(-70,70)
y_limits <- c(-30,30)
set.seed(100) 
xy <- data.frame(
  id = paste0("id_", 1:ncircles), 
  x = runif(ncircles, min(x_limits), max(x_limits)),
  y = runif(ncircles, min(y_limits), max(y_limits))) %>% 
  as_tibble()

polys <- st_as_sf(xy, coords = c(2,3)) %>% 
  st_buffer(runif(ncircles, min = 1, max = 20)) 

plot(polys[1])  

我需要导出一个sf or sp包含所有且仅由相交生成的多边形的多重多边形，例如：

（请注意，颜色只是为了举例说明预期结果，其中每个“不同颜色”区域都是一个单独的多边形，不会覆盖任何其他多边形）

我知道我可以通过一次分析一个多边形，识别并保存其所有交点，然后“擦除”这些区域形成完整的多边形并循环进行来解决问题，但这非常慢。

我觉得应该有一个更有效的解决方案，但我无法弄清楚，所以任何帮助将不胜感激！ （两个都sf and sp欢迎基于解决方案）

UPDATE:

最后，我发现即使是“一次一个多边形”，任务也远非简单！我真的在这个看似“简单”的问题上苦苦挣扎！有什么提示吗？即使是缓慢的解决方案或开始正确路径的提示也将受到赞赏！

UPDATE 2:

也许这会澄清一些事情：所需的功能将类似于此处描述的功能：

https://it.mathworks.com/matlabcentral/fileexchange/18173-polygon-intersection?requestedDomain=www.mathworks.com https://it.mathworks.com/matlabcentral/fileexchange/18173-polygon-intersection?requestedDomain=www.mathworks.com

UPDATE 3:

我将赏金授予了@shuiping-chen（谢谢！），他的答案正确地解决了所提供的示例数据集上的问题。然而，该“方法”必须推广到可能存在“四重”或“n重”交叉的情况。我将在接下来的几天里尝试解决这个问题，如果我能做到的话，我会发布一个更通用的解决方案！

Input

我稍微修改了模型数据，以说明处理多个属性的能力。

library(tibble)
library(dplyr)
library(sf)

ncircles <- 9
rmax     <- 120
x_limits <- c(-70,70)
y_limits <- c(-30,30)
set.seed(100) 
xy <- data.frame(
  id = paste0("id_", 1:ncircles), 
  val = paste0("val_", 1:ncircles),
  x = runif(ncircles, min(x_limits), max(x_limits)),
  y = runif(ncircles, min(y_limits), max(y_limits)),
  stringsAsFactors = FALSE) %>% 
  as_tibble()

polys <- st_as_sf(xy, coords = c(3,4)) %>% 
  st_buffer(runif(ncircles, min = 1, max = 20)) 
plot(polys[1])

基本操作

然后定义以下两个函数。

• cur：基础多边形的当前索引
• x：与 相交的多边形的索引cur
• input_polys：多边形的简单特征
• keep_columns：几何计算后需要保留的属性名称向量

get_difference_region()获取基础多边形与其他相交多边形之间的差异；get_intersection_region()获取相交多边形之间的交点。

library(stringr)
get_difference_region <- function(cur, x, input_polys, keep_columns=c("id")){
  x <- x[!x==cur] # remove self 
  len <- length(x)
  input_poly_sfc <- st_geometry(input_polys)
  input_poly_attr <- as.data.frame(as.data.frame(input_polys)[, keep_columns])

  # base poly
  res_poly <- input_poly_sfc[[cur]]
  res_attr <- input_poly_attr[cur, ]

  # substract the intersection parts from base poly
  if(len > 0){
    for(i in 1:len){
      res_poly <- st_difference(res_poly, input_poly_sfc[[x[i]]])
    }
  }
  return(cbind(res_attr, data.frame(geom=st_as_text(res_poly))))
}


get_intersection_region <- function(cur, x, input_polys, keep_columns=c("id"), sep="&"){
  x <- x[!x<=cur] # remove self and remove duplicated obj 
  len <- length(x)
  input_poly_sfc <- st_geometry(input_polys)
  input_poly_attr <- as.data.frame(as.data.frame(input_polys)[, keep_columns])

  res_df <- data.frame()
  if(len > 0){
    for(i in 1:len){
      res_poly <- st_intersection(input_poly_sfc[[cur]], input_poly_sfc[[x[i]]])
      res_attr <- list()
      for(j in 1:length(keep_columns)){
        pred_attr <- str_split(input_poly_attr[cur, j], sep, simplify = TRUE)
        next_attr <- str_split(input_poly_attr[x[i], j], sep, simplify = TRUE)
        res_attr[[j]] <- paste(sort(unique(c(pred_attr, next_attr))), collapse=sep)
      }
      res_attr <- as.data.frame(res_attr)
      colnames(res_attr) <- keep_columns
      res_df <- rbind(res_df, cbind(res_attr, data.frame(geom=st_as_text(res_poly))))
    }
  }
  return(res_df)
}

第一级

不同之处

让我们看看差异函数对模型数据的影响。

flag <- st_intersects(polys, polys)

first_diff <- data.frame()
for(i in 1:length(flag)) {
  cur_df <- get_difference_region(i, flag[[i]], polys, keep_column = c("id", "val"))
  first_diff <- rbind(first_diff, cur_df)
}
first_diff_sf <- st_as_sf(first_diff, wkt="geom")
first_diff_sf
plot(first_diff_sf[1])

路口

first_inter <- data.frame()
for(i in 1:length(flag)) {
  cur_df <- get_intersection_region(i, flag[[i]], polys, keep_column=c("id", "val"))
  first_inter <- rbind(first_inter, cur_df)
}
first_inter <- first_inter[row.names(first_inter %>% select(-geom) %>% distinct()),]
first_inter_sf <- st_as_sf(first_inter, wkt="geom")
first_inter_sf
plot(first_inter_sf[1])

第二级

使用第一级的交集作为输入，并重复相同的过程。

不同之处

flag <- st_intersects(first_inter_sf, first_inter_sf)
# Second level difference region
second_diff <- data.frame()
for(i in 1:length(flag)) {
  cur_df <- get_difference_region(i, flag[[i]], first_inter_sf, keep_column = c("id", "val"))
  second_diff <- rbind(second_diff, cur_df)
}
second_diff_sf <- st_as_sf(second_diff, wkt="geom")
second_diff_sf
plot(second_diff_sf[1])

路口

second_inter <- data.frame()
for(i in 1:length(flag)) {
  cur_df <- get_intersection_region(i, flag[[i]], first_inter_sf, keep_column=c("id", "val"))
  second_inter <- rbind(second_inter, cur_df)
}
second_inter <- second_inter[row.names(second_inter %>% select(-geom) %>% distinct()),]  # remove duplicated shape
second_inter_sf <- st_as_sf(second_inter, wkt="geom")
second_inter_sf
plot(second_inter_sf[1])

获取第二层的不同交集，并将它们用作第三层的输入。可以得到第三层的交集结果为NULL，那么该过程应该结束。

Summary

我们将所有差异结果放入 close 列表中，并将所有交集结果放入 open 列表中。然后我们有：

• 当打开列表为空时，我们停止进程
• 结果已接近列表

因此，我们在这里得到最终的代码（应该声明基本的两个函数）：

# init
close_df <- data.frame()
open_sf <- polys

# main loop
while(!is.null(open_sf)) {
  flag <- st_intersects(open_sf, open_sf)
  for(i in 1:length(flag)) {
    cur_df <- get_difference_region(i, flag[[i]], open_sf, keep_column = c("id", "val"))
    close_df <- rbind(close_df, cur_df)
  }
  cur_open <- data.frame()
  for(i in 1:length(flag)) {
    cur_df <- get_intersection_region(i, flag[[i]], open_sf, keep_column = c("id", "val"))
    cur_open <- rbind(cur_open, cur_df)
  }
  if(nrow(cur_open) != 0) {
    cur_open <- cur_open[row.names(cur_open %>% select(-geom) %>% distinct()),]
    open_sf <- st_as_sf(cur_open, wkt="geom")
  }
  else{
    open_sf <- NULL
  }
}

close_sf <- st_as_sf(close_df, wkt="geom")
close_sf
plot(close_sf[1])