检查 R 中的多边形交集：sf::intersects 错误地返回 TRUE 结果，而 rgeos::gIntersects 正确返回 FALSE

Question

我正在尝试检查两个多边形是否在 R 中相交。绘图时，它们显然没有相交。检查交集时，rgeos::gIntersects()当前返回FALSE，而sf::intersects()返回TRUE。我想这是因为多边形（1）大并且（2）靠得很近，所以当在平坦的表面上时它们似乎不相交，但在球体上它们似乎相交？

理想情况下，我可以将我的工作流程全部保留下来sf——但我想知道是否有一种方法可以使用sf::intersects()（或另一个sf函数）返回FALSE这里？

这是一个例子：

library(sf)
library(rgeos)
library(leaflet)
library(leaflet.extras)

#### Make Polygons
poly_1 <- c(xmin = -124.75961, ymin = 49.53330, xmax = -113.77328, ymax = 56.15249) %>% 
  st_bbox() %>%
  st_as_sfc()
st_crs(poly_1) <- 4326

poly_2 <- c(xmin = -124.73214, ymin = 25.11625, xmax = -66.94889, ymax = 49.38330) %>% 
  st_bbox() %>%
  st_as_sfc()
st_crs(poly_2) <- 4326

#### Plotting
# Visually, the polygons do not intersect
leaflet() %>%
  addTiles() %>%
  addPolygons(data = poly_1) %>%
  addPolygons(data = poly_2)

#### Check Intersection
# returns FALSE
gIntersects(poly_1 %>% as("Spatial"), 
            poly_2 %>% as("Spatial"))

# returns TRUE
st_intersects(poly_1,
              poly_2,
              sparse = F)

这是多边形，它们在视觉上不相交。

Answer 1

这是一个有趣的问题，其根本原因是平面（在平坦的表面上）和球形（在地球仪上）几何形状之间的差异。

\n

在平面上 - 这是 GEOS 所采用的简化方法 - 多边形的四个角由四条直线连接，角度之和为 360\xc2\xb0 度等。几何学的原理就像欧几里得很久以前所教的那样。

\n

但是，至关重要的是，世界并不是这样运作的。在地球仪上，多边形的四个连接点不是直线，而是大圆。或者更确切地说，当它们在地球上绘制时是直的，而当在平面（例如地图或计算机屏幕）上滚平时它们是弯曲的。

\n

因为一个例子超过 1000 个字，考虑这段代码：

\n

library(sf)\nlibrary(dplyr)\n\n# make polygons\npoly_1 <- c(xmin = -124.75961, ymin = 49.53330, xmax = -113.77328, ymax = 56.15249) %>% \n  st_bbox() %>%\n  st_as_sfc()\nst_crs(poly_1) <- 4326\n\npoly_2 <- c(xmin = -124.73214, ymin = 25.11625, xmax = -66.94889, ymax = 49.38330) %>% \n  st_bbox() %>%\n  st_as_sfc()\nst_crs(poly_2) <- 4326\n\n# this is what you *think* you see (and what GEOS sees, being planar) \n# = four corners connected by straight lines\n# & no intersecton\nmapview::mapview(list(poly_1, poly_2))\n

\n

\n

# this is what you *should* see (and what {s2} sees, being spherical)\n# = four corners connected by great circles\n# & an obvious intersection around the lower right corner of the small polygon\npoly_1alt <- st_segmentize(poly_1, units::set_units(1, degree))\npoly_2alt <- st_segmentize(poly_2, units::set_units(1, degree))\n\nmapview::mapview(list(poly_1alt, poly_2alt))\n

\n

\n

您有两个选择：

\n

\n
• 接受你对多边形的想法是错误的，并拥抱球形，即{s2}逻辑。
\n

\n

理论上这应该是正确的方法，但有些违反直觉。

\n

\n
• 放弃{sf}对多边形的球形方法，并强制其应用平面方法（例如GEOS使用）。
\n

\n

{sf}从理论上讲，这将是一种错误的方法，但与您的平面直觉和大多数 GIS 工具之前的行为（包括1.0.0 版本之前的版本）一致

\n

# remedy (of a kind...) = force planar geometry operations\n\nsf::sf_use_s2(FALSE)\n\nst_intersects(poly_1, poly_2, sparse = F)\n#      [,1]\n# [1,] FALSE\n

\n