Since the update of the sf
package (V1.0) simple operations such as st_join()
or st_intersection()
may be hampered by features with invalid spherical geometry
such as:
[XY] Loop XY is not valid: Edge XY crosses edge XY
[XY] Loop XY is not valid: Edge XY has duplicate vertex with edge XY.
A quick fix and second best solution is to turn off s2 by sf::sf_use_s2(FALSE)
, but ideally would be to fix the features with invalid spherical geometry, so that s2 can process it.
How to do in R?
Related questions can be found here and here, but offer only work-arounds to the features with invalid spherical geometry
problem.
To reproduce the problem download the following shapefile here, and run:
path_to_shape_file_NUTS = "NUTS_RG_01M_2016_4326_LEVL_3.shp"
shape_file_NUTS_WGS84 = st_make_valid(st_read(path_to_shape_file_NUTS, quiet=TRUE))
y_coord <- c(46.976, 46.948)
x_coord <- c(7.483, 7.45)
coordinates_as_df <- data.frame(x_coord, y_coord)
coordinates_WGS84 <- st_as_sf(x = coordinates_as_df,
coords = c("x_coord", "y_coord"),
crs = 4326)
overlay <- st_join(coordinates_WGS84, shape_file_NUTS_WGS84, join=st_intersects)
print(overlay)
Best Answer
I realize it is a workaround, and not a generic solution. But issues with invalid geometries are by nature tricky (All happy geometries are alike; each unhappy geometry is unhappy in its own way, to misquote Anna Karenina) and a truly generic solution may be not feasible.
If 2016 level 3 NUTS from GISCO is the only shapefile you are after you may consider accessing it via
{giscoR}
- it should in theory be the same thing, but accessed via API and known to be valid.Consider this code: