【问题标题】:R sf: st_intersection on a list of polygonsR sf:多边形列表上的 st_intersection
【发布时间】:2020-12-29 20:30:57
【问题描述】:

我有一个名为 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中所有多边形的公共交集之外,还衍生出很多其他的多边形。

这三种方法都不能提供令人满意的解决方案。在并行化框架中循环通过成对比较的层次结构可能会更快。但是,我认为还有比这更简单、更有效的解决方案。

欢迎任何cmets和建议。

给昨天关闭此问题的人的说明:请勿关闭此问题。如果您个人对此有任何问题,请发表评论或给我发私信。但不要关闭它。

【问题讨论】:

    标签: r gis geospatial sf


    【解决方案1】:

    我不认为遍历列表的开销在这里是一个问题:找到多个多边形的交集只是计算成本很高。但是,使用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 确保所有多边形实际相交。这是一个更快的计算,前提是没有不规则的交叉点。

    【讨论】:

    • 谢谢。根据您的评论,我决定使用类似于purrr::accumulate()purrr::reduce()。但它不会存储所有多边形,而是只生成感兴趣的多边形 - 您示例中的红色正方形。
    猜你喜欢
    • 2022-01-03
    • 1970-01-01
    • 2019-09-30
    • 2019-03-02
    • 2019-11-22
    • 1970-01-01
    • 1970-01-01
    • 2020-11-14
    • 1970-01-01
    相关资源
    最近更新 更多