CODEWITHSUNDEEP
rr-sf
mharinga
mharinga

Reputation: 1780

sf::st_intersection: error in CPL_nary_intersection(x)

Lets say I have the following data with lon, lat pairs:

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

Then I would like to create a circle around each point:

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

Create plot:

plot(st_geometry(circle_df))

Now st_intersection() returns an error:

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

However, if I only use the first four rows I do not get an error.

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)

I am not sure why I get the error?

Created on 2021-04-29 by the reprex package (v1.0.0)

Upvotes: 1

Views: 503

Answers (1)

Klaus
Klaus

Reputation: 484

# 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()

Upvotes: 1

Related Questions