logo标签列表

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

rpolygonintersectionr-spr-sf
27
从包含相当数量(约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元组”相交的情况。如果我成功，我将尝试在未来几天内解决这个问题，并发布更通用的解决方案！

你对所有子多边形都感兴趣，还是只对由交集产生的一个感兴趣？（在示例中，底部中心的“group”，你想要绿色、蓝色和橄榄色，还是只想要蓝色？） - GGamba
我对所有的子多边形都很感兴趣。然而，那些与其他任何东西都没有交集的“区域”可以使用每个多边形和它们的并集之间的sym_difference运算符简单地提取出来。我正在努力解决的是交叉点，特别是多个交叉点。 - lbusett
@LoBu，请看一下这个链接，我也在尝试它。 - parth
3个回答
15

输入

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

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，然后该过程应该结束。

总结

我们将所有的差异结果放入关闭列表中，将所有的交集结果放入打开列表中。然后我们有：

  • 当打开列表为空时，我们停止该过程
  • 结果是闭合列表

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

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

final result

enter image description here

非常有趣，谢谢！我会尽快测试它。顺便问一下，在我开始之前：如果我们有“四重”甚至“五重”交叉点，这个程序也能处理得好吗？还是它只局限于最多三个重叠的多边形？ - lbusett
我认为这对一般情况应该是可行的，因为重复的过程与“最大三重叠”没有关系。但是记得使用“独特”的交集作为下一级的输入。我花了一些时间才想明白。 - Shuiping Chen
好的，谢谢！我会在不同的条件下进行测试，并回复您！ - lbusett
嗨。我今天才得以看一下这个。对于手头的问题，它运行得非常好，谢谢！但是，这需要为具有超过“三重”交叉点的情况进行概括。我会在接下来的几天里研究一下它！ - lbusett
9

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

太好了。我刚刚在一个包含1000个多边形的测试数据集上尝试了一下，它完美地解决了这个问题！谢谢！ - lbusett
@Edzer Pebesma，我因为自相交而得到了评估错误。在sf包中有没有简单的方法来解决这个问题？ - JdP
在这里找到答案：https://github.com/r-spatial/sf/issues/347：使用 st_make_valid 或 st_buffer(x, 0)。 - JdP
2
不确定这是否对您有帮助，因为它不是R语言，但我认为使用Python有一种很好的方法来解决这个问题。有一个叫做GeoPandas（http://geopandas.org/index.html）的库，它允许您轻松地进行地理操作。您需要执行以下步骤：

  1. 将所有多边形加载到geopandas GeoDataFrames中
  2. 循环所有GeoDataFrames，运行联合叠加操作（http://geopandas.org/set_operations.html）

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

操作之前 - 2个多边形

2 polygons

手术后 - 9个多边形

9 polygons

如果有任何不清楚的地方，请随时让我知道！希望能帮到你！

@LoBu 你是对的！这并不像我想象的那么简单。另一种方法是为每个多边形分配一个数字1值，然后在所有子多边形上进行聚合。通过这样做，您将为每个子多边形生成第三个维度。完成此操作后，您只需要过滤掉所有仍然具有分配值等于1的子多边形即可。每个重叠的子多边形将具有编号2或更高。我还没有花太多时间玩地理查询，以确保这种方法是否可行或者运行数据集的速度是否足够快。你认为呢？ - fernandosjp
嗨。抱歉，您所说的“聚合”是什么意思？ - lbusett
我的意思是将分配给每个多边形的数字1相加。以您的图形为例，中间三个大圆之间的重叠区域将具有值等于3；2个圆之间的重叠区域将具有值等于2，依此类推。实际上，这将是每个重叠区域的计数。 - fernandosjp
没问题！你可能应该看一下如何对多边形进行空间操作。在这个链接中有一个非常相似的案例讨论：https://gis.stackexchange.com/questions/190648/summing-attribute-values-for-areas-where-multiple-polgons-overlap-using-arcgis-d - fernandosjp
1
很有前途。现在的问题是如何通过“R”脚本来实现这一点。有什么提示吗？ - lbusett
显示剩余3条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接