How to resolve spherical geometry failures when joining spatial data

Question

I have a shapefile (with several polygons) and a dataframe with coordinates. I want to assign each coordinate in a dataframe to a polygon in a shapefile. So to add a column in a data frame with a polygon name or id Here is the link to the data

library(sf)
library(readr)
shape <- read_sf("data/Provinces_v1_2017.shp")
data<- read_csv("data/data.csv")

But when I try to join them, I always get the error

pts = st_as_sf(data, coords = c("dec_lon", "dec_lat"), crs= 4326)

st_join(pts, shape)

i tried over() functions, and other tricks like st_make_valid() but I always get this error: Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : Evaluation error: Found 30 features with invalid spherical geometry.

It is a recent issue (before my code worked), but now I am unable to use sf package to do this task, I always end up with this error. I updated the libraries to see whether it would help, but I could not make it work.

I would really appreciate your help on this matter

Answer 1

You have two options:

1. turn off the s2 processing via sf::sf_use_s2(FALSE) in your script; in theory the behaviour should revert to the one before release 1.0
2. repair the spherical geometry of your polygons object; this will depend on the actual nature of your errors.

I can't access your file & make certain, but this piece of code has helped me in the past:

yer_object$geometry <- yer_object$geometry %>%
  s2::s2_rebuild() %>%
  sf::st_as_sfc()

Answer 2

I find that this 'invalid spherical geometry' does keep on popping up. If the s2::s2_rebuild() solution above doesn't work, a solution which usually works for me involves projecting and simplifying (reducing the map resolution a little). If your application can work with less resolution try this.

library(tidyverse)
library(sf)
crs_N = 3995 #northern polar projection

# example of FAILING map - with bad spherical geometry.
m_RU <- rnaturalearthdata::countries50 %>% 
  st_as_sf() %>% 
  filter((admin %in% c("Russia") )) |> 
  st_as_s2()

In the example, I chose Russia because it crosses the dateline, which can be one of the challenges. I switch to an Arctic polar projection, and reduce the map to 10km resolution (5km is not enough in this case!).

# with 2 extra lines the problem is gone
m_RU <- rnaturalearthdata::countries50 %>% 
  st_as_sf() %>% 
  filter((admin %in% c("Russia") )) |> 
  st_transform(crs = crs_N)   |>  
  st_simplify(dTolerance = 10000) |> # to get rid of duplicate vertex (reduce to 10km steps)
  st_as_s2()