为啥使用 st_intersection 而不是 st_intersects？

Posted 2023-02-16

【中文标题】为啥使用 st_intersection 而不是 st_intersects？

【英文标题】：Why use st_intersection rather than st_intersects?为什么使用 st_intersection 而不是 st_intersects？

【发布时间】：2020-10-08 01:09:33

【问题描述】：

st_intersection 与st_intersects 相比非常慢。那么为什么不使用后者而不是前者呢？这是一个带有小型玩具数据集的示例，但对于我的实际集合只有 62,020 个与实际地理区域边界相交的点，执行时间的差异是巨大的。我有 24Gb 的 RAM，st_intersects 代码需要几秒钟，而st_intersection 代码需要超过 15 分钟（可能更多，我没有耐心等待......）。 st_intersection 是否会做我无法使用 st_intersects 做的事情？

以下代码处理sfc 对象，但我相信sf 对象同样适用。

library(sf)
library(dplyr)

# create square
s <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1)) %>% list %>% st_polygon %>% st_sfc
# create random points
p <- runif(50, 0, 11) %>% cbind(runif(50, 0, 11)) %>% st_multipoint %>% st_sfc %>% st_cast("POINT")

# intersect points and square with st_intersection
st_intersection(p, s)

# intersect points and square with st_intersects (courtesy of https://***.com/a/49304723/7114709)
p[st_intersects(p, s) %>% lengths > 0,]

【参考方案1】：

答案是，通常这两种方法做不同的事情，但在您的特定情况下（查找点集合和多边形的交集），st_intersects 可用于有效地完成相同的工作。

我们可以通过您自己修改的简单示例来展示差异。我们从一个正方形开始：

library(sf)
library(dplyr)

# create square
s <- rbind(c(1, 1), c(10, 1), c(10, 10), c(1, 10), c(1, 1)) %>% 
  list %>% 
  st_polygon %>% 
  st_sfc

plot(s)

现在我们将创建一个矩形并用虚线轮廓将其绘制在同一个图上：

# create rectangle
r <- rbind(c(-1, 2), c(11, 2), c(11, 4), c(-1, 4), c(-1, 2)) %>% 
  list %>% 
  st_polygon %>% 
  st_sfc

plot(r, add= TRUE, lty = 2)

现在我们找到两个多边形的交点并将其绘制为红色：

# intersect points and square with st_intersection
i <- st_intersection(s, r)

plot(i, add = TRUE, lty = 2, col = "red")

当我们检查对象i时，我们会看到它是一个新的多边形：

i
#> Geometry set for 1 feature 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 1 ymin: 2 xmax: 10 ymax: 4
#> epsg (SRID):    NA
#> proj4string:    NA
#> POLYGON ((10 4, 10 2, 1 2, 1 4, 10 4))

然而，如果我们使用st_intersects，我们只会得到一个逻辑结果，告诉我们r 和s 之间是否确实存在交集。如果我们尝试使用它对r 进行子集化以找到交点，我们不会得到相交的形状，我们只会得到原来的矩形：

r[which(unlist(st_intersects(s, r)) == 1)]
#> Geometry set for 1 feature 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -1 ymin: 2 xmax: 11 ymax: 4
#> epsg (SRID):    NA
#> proj4string:    NA
#> POLYGON ((-1 2, 11 2, 11 4, -1 4, -1 2))

您的情况有所不同，因为您试图找到与多边形相交的点的子集。是不是这种情况，一组点与一个多边形的交集与满足条件st_intersects的子集相同。

所以很高兴您找到了一种获得更快交叉路口的有效方法。请注意，这只适用于与多边形相交的点集合。

【讨论】：

非常清晰，追溯起来很明显！我应该认为我的用例在地理学家中很常见，所以我还会用“地理空间”标记我的问题。