按组查找多个多边形的交集答案

【问题标题】：Find intersection of multiple polygons by group按组查找多个多边形的交集
【发布时间】：2021-12-06 06:09:11
【问题描述】：

因此，使用 dplyr 中的 st_union 和 group_by() 找到两个以上多边形的并集是很简单的。找到与 st_intersection 的交集并不是那么简单，因为 st_intersection 仅适用于寻找两个对象的交集。按组查找多个多边形的交集的最佳方法是什么？这是一些示例数据

s1 <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1))
s2 <- s1 + 4
s3 <- s1 - 4

group <- c(rep(1, 3), rep(2, 3))

df <- data.frame("group" = c(rep(1, 3), rep(2, 3)))

df <- data.frame(cbind(df, rep((c(st_sfc(st_polygon(list(s1))), st_sfc(st_polygon(list(s2))), st_sfc(st_polygon(list(s3))))), 2)))

我想要的是将几何列替换为每个组中的交集，即由

创建的对象


p <- list(s1 = s1, s2 = s2, s3 = s3)

p  <-  lapply(p, function(x) st_sfc(st_polygon(list(x))))

intersection <- accumulate(p, st_intersection)$s3

intersection <- st_sfc(st_polygon(intersection))

所以最终的数据看起来像


df <- data.frame("group" = c(rep(1, 3), rep(2, 3)))

df <- data.frame(st_sf(cbind(df, rep(st_sfc(st_polygon(intersection)), 6))))

【问题讨论】：

标签： r dplyr intersection sf

【解决方案1】：

设法通过循环来完成，但更喜欢 dplyr 解决方案

agg <- as.data.frame(matrix(NA, 0, 2))
for (i in 1:nrow(df)) {
  group <- df$group[i]
  vec <- df$geometry[df$group == group]
  my.list <- as.list(vec)
  output <- Reduce(st_intersection, my.list) %>% st_geometry()
  my.df <- data.frame("group" = group)
  res <- as.data.frame(st_sf(cbind(my.df, output)))
  
  agg <- rbind(agg, res)
}

【讨论】：

【解决方案2】：

首先，st_intersection() 只能与单个 R 对象一起使用，例如，如果您为其提供 sf collection，它将返回它的自相交；查看 sf 页面 here 了解更多详情，以及另一个答案作为示例 here。生成的collection 包括features，其属性包括n.overlaps（层数）和origins（哪些多边形对每个重叠区域有贡献）。

因此，对于本示例，我们使用 st_intersection() 并仅保留 n.overlaps == 3 所在的几何图形（因为每组中有 3 个几何图形）。基本上步骤是：将df变成sf collection，然后我们使用split()将sf collection拆分为group，找到每个组的自交点，将两个组重新连接在一起，然后只保留有 3 个重叠的区域。

# example data from question ~~~~~

library(sf)
library(tidyverse)

s1 <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1))
s2 <- s1 + 4
s3 <- s1 - 4

group <- c(rep(1, 3), rep(2, 3))
df <- data.frame("group" = c(rep(1, 3), rep(2, 3)))
df <- data.frame(cbind(df, rep((c(st_sfc(st_polygon(list(s1))), st_sfc(st_polygon(list(s2))), st_sfc(st_polygon(list(s3))))), 2)))


# solution ~~~~~

df %>%
  st_as_sf %>% 
  split(., .$group) %>% 
  map(st_intersection) %>% 
  bind_rows %>% 
  filter(
    n.overlaps == 3
  )

# Simple feature collection with 2 features and 3 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 5 ymin: 5 xmax: 6 ymax: 6
# CRS:           NA
#     group n.overlaps origins                       geometry
# 1.3     1          3 1, 2, 3 POLYGON ((5 5, 5 6, 6 6, 6 ...
# 4.3     2          3 1, 2, 3 POLYGON ((5 5, 5 6, 6 6, 6 ...

不幸的是，这种方法并不完全是dplyr，因为它使用split() 和map() 来遍历每个group，但是它是可管道的。我认为group_by() 可能会代替split() %>% map() %>% bind_rows()，但不幸的是我无法让st_intersection() 尊重group_by()。

希望这仍然有用，也许有人可以找到将split() 和map() 替换为group_by() 的方法。

【讨论】：