sf::st_intersection：CPL_nary_intersection(x) 中的错误答案

【问题标题】：sf::st_intersection: error in CPL_nary_intersection(x)sf::st_intersection：CPL_nary_intersection(x) 中的错误
【发布时间】：2021-04-29 16:16:37
【问题描述】：

假设我有以下带有 lon、lat 对的数据：

# Data frame
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"
                             ))

然后我想围绕每个点创建一个圆圈：

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"))

创建情节：

plot(st_geometry(circle_df))

现在st_intersection() 返回错误：

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

但是，如果我只使用前四行，我不会收到错误。

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 package (v1.0.0) 于 2021 年 4 月 29 日创建}

【问题讨论】：

这应该开箱即用，有一个 GitHub 问题：github.com/r-spatial/sf/issues/1668。

标签： r sf

【解决方案1】：

# Data frame
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"
                             ))


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")) %>%
  lwgeom::st_snap_to_grid(50) # using snap to grid you can avoid the error

plot(st_geometry(circle_df))

points_within <- st_intersection(circle_df) %>%
  st_centroid()

【讨论】：