R中sf与rgeos多边形相交判断结果不一致，如何用sf获正确结果？

我在R中判断两个大多边形是否相交时遇到问题：可视化显示二者无交集，但rgeos::gIntersects返回FALSE，sf::st_intersects却返回TRUE。推测原因是多边形跨度大且距离接近，平面几何与球面几何的计算差异导致的。希望保持sf工作流，请问如何用sf函数得到正确的FALSE结果？

示例代码：

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

#### 创建多边形
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

#### 可视化（无交集）
leaflet() %>%
  addTiles() %>%
  addPolygons(data = poly_1) %>%
  addPolygons(data = poly_2)

#### 交集判断
# 返回 FALSE
gIntersects(poly_1 %>% as("Spatial"), 
            poly_2 %>% as("Spatial"))

# 返回 TRUE
st_intersects(poly_1,
              poly_2,
              sparse = F)

可视化效果：

解决方案 #

方法1：投影到平面坐标系

sf在处理WGS84（EPSG:4326）时默认用球面几何计算，大跨度多边形容易出现偏差。可以将多边形投影到合适的平面坐标系（比如Albers等面积投影、UTM分带投影），再进行交集判断：

# 选择覆盖北美地区的Albers等面积投影
aea_crs <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

poly_1_proj <- st_transform(poly_1, crs = aea_crs)
poly_2_proj <- st_transform(poly_2, crs = aea_crs)

# 此时计算交集返回 FALSE
st_intersects(poly_1_proj, poly_2_proj, sparse = FALSE)

方法2：强制使用平面几何计算

如果不需要考虑球面曲率，可将几何对象转换为无CRS的XY平面几何，强制sf用GEOS平面算法计算：

# 移除CRS，转换为XY平面几何
poly_1_xy <- st_geometry(poly_1) %>% st_set_crs(NA)
poly_2_xy <- st_geometry(poly_2) %>% st_set_crs(NA)

# 计算交集返回 FALSE
st_intersects(poly_1_xy, poly_2_xy, sparse = FALSE)

原因解释 #

• rgeos::gIntersects始终基于平面几何计算，直接把经纬度当作平面坐标处理；
• sf::st_intersects在处理WGS84这类地理坐标系时，默认采用球面几何计算，大跨度多边形的边缘在球面模型中可能被判定为相交，和视觉上的平面投影效果产生差异。

内容的提问来源于stack exchange，提问作者Rob Marty