R/sf: Some edges are crossing in a MULTIPOLYGON. How to make it valid when using S2 geometry?

Question

Goal: I would like to know the NUTS region (2016 version) of points in Europe.

Method: (1) Download the official geodata from Eurostat, (2) Use st_join from sf with R.

Problem: I've updated sf, now using v1.0. Now that it uses S2 geometry, I have an error message:

Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
  Evaluation error: Found 2 features with invalid spherical geometry.
[821] Loop 0 is not valid: Edge 299 crosses edge 301
[1102] Loop 0 is not valid: Edge 575 has duplicate vertex with edge 579.

Note that it was not generating an error before the update of sf.

Part of a solution: I've used st_make_valid. Now I have a shorter error message:

Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
  Evaluation error: Found 1 feature with invalid spherical geometry.
[821] Loop 0 is not valid: Edge 299 crosses edge 301.

Minimal example:

library(sf)
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)

You should be able to get the shapefile using this link, then opening the zip file "NUTS_RG_01M_2016_4326_LEVL_3.shp.zip" in the zip file.

My question: How to get rid of the error message and get the NUTS region?

Answer 1

The answer by Jindra Lacko is specific to NUTS regions. Just a quick fix if you have an invalid geometry after updating sf to v1.0:

sf_use_s2(FALSE)

It does not solve the problem, but it brings back the old way of working. I still don't have a general solution to make a geometry valid when using S2, but just if it could help someone...

Answer 2

While the specific dataset mentioned in the question works for me with the latest sf (1.0-2), I had basically the same issue with other EU datasets, namely those on local administrative units (LAU) imported via shapefiles for 2019 and 2020 (previous years seem to be fine).

sf::st_make_valid() does not help, but projecting them with sf::st_transform(crs = 3857) (and then back to sf::st_transform(crs = 4326) if needed) fixed this for me.

Answer 3

You can consider getting the NUTS regions via {giscoR} package - it serves formally valid NUTS regions.

nuts <- giscoR::gisco_get_nuts(year = '2016')
all(sf::st_is_valid(nuts))
[1] TRUE

The example is now:

library(sf)
library(giscoR)
NUTS_WGS84 = gisco_get_nuts(year = '2016')
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, NUTS_WGS84, join=st_intersects)
print(overlay[overlay$LEVL_CODE == 3, ])

And the output:

Simple feature collection with 2 features and 9 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.45 ymin: 46.948 xmax: 7.483 ymax: 46.976
Geodetic CRS:  WGS 84
  NAME_LATN MOUNT_TYPE CNTR_CODE URBN_TYPE COAST_TYPE   FID NUTS_NAME NUTS_ID LEVL_CODE             geometry
4      Bern          3        CH         2          0 CH021      Bern   CH021         3 POINT (7.483 46.976)
8      Bern          3        CH         2          0 CH021      Bern   CH021         3  POINT (7.45 46.948)