我有一个名为 p_1、p_2、...、p_n 的多边形(和多边形)列表。我想获得它们相交的区域。由于 st_intersection()
不接受列表作为参数,我尝试了以下三种方法。它们都没有提供令人满意的解决方案,这就是为什么我正在寻找替代的、更有效的技术。
(i) 我可以遍历列表
for(i in P) p_1 <- st_intersection(p_1, i)
其中 P 是包含多边形 p_2 到 p_n 的列表。但这相当慢。
(ii) do.call()
方法,即
p <- do.call(st_intersection, P)
其中 P 是多边形 p_1 到 p_n 的列表,仅计算列表中前两个多边形之间的交集。
(iii) 我可以将多边形组合成一个 sf 对象,然后运行 st_intersection()
:
p <- do.call(c, P) %>%
st_sf() %>%
st_intersection()
它有效但速度很慢。大概是因为除了P中所有多边形的公共(public)交集之外,它还派生出很多其他多边形。
这三种方法都没有提供令人满意的解决方案。在并行化框架中循环遍历成对比较的层次结构可能会更快。但是,我认为有比这更简单、更有效的解决方案。
欢迎大家提出意见和建议。
给昨天关闭这个问题的人的注意事项:不要关闭这个问题。如果您个人对此有疑问,请发表评论或给我发私信。但不要关闭它。
最佳答案
我认为遍历列表的开销在这里不是问题:找到多个多边形的交集只是计算量大。但是,使用 purrr::accumulate
可以轻松管理按顺序将函数应用于列表成员的方法(实际上是您尝试使用 do.call
执行的操作):
您没有可重现的示例供这里的人测试可能的解决方案,并且从头开始创建 sf 多边形涉及一些工作,因此可能是您之前的问题被关闭的原因 - 我不知道不知道。
无论如何,让我们在列表中创建三个重叠的正方形并绘制它们:
library(sf)
library(purrr)
# create square
s1 <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1))
p <- list(s1 = s1, s2 = s1 + 4, s3 = s1 - 4)
p <- lapply(p, function(x) st_sfc(st_polygon(list(x))) )
plot(p[[1]], xlim = c(-5, 15), ylim = c(-5, 15))
plot(p[[2]], add = TRUE)
plot(p[[3]], add = TRUE)
我们的目标是找到所有三个正方形的交点,当然是中间的小正方形。使用 purrr
,这很简单:
intersection <- accumulate(p, st_intersection)$s3
所以当我们添加我们的结果时,颜色为红色,我们得到:
plot(intersection, col = "red", add = TRUE)
就性能而言,accumulate
仅比原始循环快 10% 左右,因此如果性能是个大问题,您可能需要将其并行化。此外,如果所有多边形之间可能没有交点,您可以找到最小的多边形并使用 st_intersects
来确保所有多边形实际上与它相交。这是一个更快的计算,前提是很有可能没有唯一的交叉点。
关于R SF : st_intersection on a list of polygons,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/63844990/