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

Posted 2023-02-14

【中文标题】高效提取 MultiPolygon 中自相交特征生成的所有子多边形

【英文标题】：Efficient extraction of all sub-polygons generated by self-intersecting features in a MultiPolygon

【发布时间】：2017-11-21 16:45:32

【问题描述】：

从包含相当多（大约 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 或 sp 多多边形，其中包含所有且仅包含交叉点生成的多边形，例如：

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

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

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

更新：

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

更新 2：

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

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

更新 3：

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

【问题讨论】：

您是对所有的子多边形感兴趣还是只对交叉点产生的子多边形感兴趣？ （在示例中，底部中心“组”，您想要绿色、蓝色和橄榄色，还是只想要蓝色？）

我对所有的子多边形都感兴趣。但是，不与任何东西相交的“区域”很容易在每个多边形及其联合之间使用 sym_difference 运算符提取。 I0m 苦苦挣扎的是交叉口，尤其是 multiple 交叉口。

@LoBu ，看看this ，我也在尝试解决它

【参考方案1】：

不确定它是否对您有帮助，因为它不在 R 中，但我认为有一种使用 Python 解决此问题的好方法。有一个名为 GeoPandas (http://geopandas.org/index.html) 的库可以让您轻松地进行地理操作。在步骤中，您需要执行以下操作：

将所有多边形加载到 geopandas GeoDataFrames 循环所有运行联合叠加操作的 GeoDataFrame (http://geopandas.org/set_operations.html)

文档中显示了确切的示例。

操作前 - 2 个多边形

操作后 - 9 个多边形

如果有任何不清楚的地方，请随时告诉我！希望对您有所帮助！

【讨论】：

您好，感谢您的回答。如果我正确解释文档，这将计算两个不同“多边形”之间的 set_operations，这也可以在“R”中使用sf 或rgeos 完成。相反，我的问题是我有大量可以相交的“单个多边形”，并且相交区域可以随机“重叠”。这会生成额外的“子多边形”，这是我想“提取”的。例如，看看示例图像中间的三个“大”圆圈会发生什么。

@LoBu 你是对的！这并不像我想的那么简单。您可以采用的另一种方法是为每个多边形分配一个数字 1 值，然后对所有多边形进行聚合。这样，您将在每个子多边形中生成第三个维度。完成此操作后，您只需要过滤掉所有仍保持分配值等于 1 的子多边形。每个重叠的子多边形的编号为 2 或更高。我没有时间玩很多地理查询来确保这种方法有效，或者它是否足够快来运行你的数据集。你怎么看？

嗨。抱歉，“聚合”是什么意思？

我的意思是将分配给每个多边形的数字 1 相加。所以以你图为例。中间三个大圆圈之间的重叠值等于 3； 2个圆圈之间的重叠将具有等于2的值，依此类推。实际上，这将是每个重叠区域的计数。

这很有希望。现在的问题是通过“R”中的脚本找到一种方法来做到这一点。有什么提示吗？

【参考方案2】：

输入

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

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，则流程结束。

总结

我们把所有的差结果放入close list，把所有的交集结果放入open list。然后我们有：

当打开列表为空时，我们停止进程 结果是关闭列表

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

# 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])

【讨论】：

非常有趣，谢谢！我会尽快测试它。顺便说一句，在我开始研究之前：在我们有“四重”甚至“五重”交叉点的情况下，这是否也能很好地处理，还是“限制”最多三个重叠的多边形？

我认为它适用于一般情况，因为重复的过程与“最大三个重叠”无关。但请记住使用 distinct 交叉点作为下一级的输入。我花了一些时间才弄明白。

好的，谢谢！我会在不同的条件下进行测试并回复您！

嗨。我只能在今天看到这个。非常适合手头的案例，谢谢！然而，这将需要推广到“三重”交叉点以上的情况。将在接下来的几天内调查它！

【参考方案3】：

这现在已在 R 包 sf 中实现，作为使用单个参数（sf 或 sfc）调用 st_intersection 时的默认结果，有关示例，请参见 https://r-spatial.github.io/sf/reference/geos_binary_ops.html。 （我不确定origins 字段是否包含有用的索引；理想情况下，它们应该只指向x 中的索引，现在它们有点自我引用）。

【讨论】：

太棒了。我刚刚在一个 1000 个多边形的测试数据集上进行了尝试，它完美地解决了这个问题！谢谢！

@Edzer Pebesma，由于自相交，我得到了评估错误。 sf包里有没有简单的方法解决这个问题？

在这里找到答案：github.com/r-spatial/sf/issues/347：使用st_make_valid 或st_buffer(x, 0)