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

Hereunto

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

Limpid

you are right! I tried sf::sf_use_s2(FALSE) and it is working now! – Hereunto

Super, glad to be of service! :) – Limpid

I stumbled across this error too and it was driving me crazy because the operation is running smooth in python using the same inputs. @jiladata Would you be so kind to explain what changed in release 1.0 to cause the unexpected behavior? Thank you – Thereafter

Hi @IvanP! The change in sf v1.0 was move from backend engine for unprojected coordinates (the geographic ones, i.e. lat - long as in EPSG 4326) from GEOS to s2 from Google. GEOS treats projected coordinates as planar (i.e. two points lie on a line of infinite max lenght) while s2 is more "correct" (two points lie on a great circle of circumference of 40 075 kilometers). The change of default backend had implications, as both GEOS and s2 are making shortcuts and taking (different) assumptions. Have a look at r-spatial.github.io/sf/articles/sf7.html for more information – Limpid

Thank you for the clear explanation! I had a quick look at the documentation of sf v1.0 and ended up at the same page you linked. – Thereafter

:) as was to be expected, since I linked an official communication by the {sf} maintainers :) – Limpid

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

Linn