R语言sf包执行st_intersection时出现CPL_nary_intersection的GEOS报错排查

复现步骤 #

• 准备包含经纬度配对的样本点数据框，共包含5个样本点，分别存储纬度lat、经度lon字段：

# 样本数据框
df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548, 
                             53.2142084420887, 53.2051969857393), lon = c(6.57022923569802, 
                                                                          6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
                             )), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
                             ))

• 以每个样本点为中心，生成半径50公里的圆形缓冲区，实现代码如下：

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
circle_df <- df %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(3035) %>%
  st_buffer(dist = units::set_units(50, "kilometers"))

异常现象 #

生成的5个缓冲区可视化时可正常显示圆形分布效果，但对全量5个要素执行多要素交集计算时，直接抛出运行错误：

points_within <- st_intersection(circle_df) %>%
  st_centroid()
#> Error in CPL_nary_intersection(x): GEOS exception

对照测试结果 #

仅选取前4行缓冲区要素，执行完全一致的交集计算流程时，可正常输出计算结果，无报错：

st_intersection(circle_df[1:4,]) %>%
  st_centroid()
#> Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
#> geometries of x
#> Simple feature collection with 14 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 4043254 ymin: 3307392 xmax: 4139992 ymax: 3395156
#> Projected CRS: ETRS89-extended / LAEA Europe
#> # A tibble: 14 x 3
#>    n.overlaps origins            geometry
#>  *      <int> <list>          <POINT [m]>
#>  1          1 <int [1]> (4136096 3329095)
#>  2          2 <int [2]> (4133590 3378661)
#>  3          1 <int [1]> (4069474 3387378)
#>  4          2 <int [2]> (4139992 3355924)
#>  5          3 <int [3]> (4109810 3392076)
#>  6          2 <int [2]> (4092517 3395156)
#>  7          1 <int [1]> (4137083 3371188)
#>  8          2 <int [2]> (4090567 3307392)
#>  9          3 <int [3]> (4112906 3316653)
#> 10          3 <int [3]> (4046277 3332453)
#> 11          4 <int [4]> (4091271 3351759)
#> 12          2 <int [2]> (4043254 3347152)
#> 13          3 <int [3]> (4047643 3372812)
#> 14          1 <int [1]> (4072727 3310143)

当前诉求 #

目前暂不明确触发该报错的根本原因，需要对应的问题成因解释与可落地的排查思路。

注：本示例基于reprex包v1.0.0生成，创建于2021-04-29

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