How to Resolve Spherical Geometry Failures When Joining Spatial Data

How to resolve spherical geometry failures when joining spatial data

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

How to fix spherical geometry errors caused by conversion from GEOS to s2

This is not a problem with code, the issue is with data. S2 is just more strict about polygon conformance, and throws this errors when it encounters an invalid polygon.

The polygons here seem to have self-intersection, like

A--B
|  |
D--C--E
   |  |
   G--F

This shape should be described as two polygons, ABCDA and CEFGC. But often it is described as a single polygon ABCEFGCDA (note C is repeated twice) - which probably happened here too. Some libraries would happily accept this, but S2 complains about duplicate vertex C in non-consecutive edges BC and GC.

If you can, fix the data, before loading it. I know PostgreSQL/PostGIS can usually fix these - it would accept the input WKB, and has ST_MakeValid that will split this into multiple polygons. R seem to have st_make_valid https://rdrr.io/cran/sf/man/valid.html too (I don't have personal experience with it).

How to solve this problem with the st_distance function from sf giving an error when inspecting the data?

Hard to do with the sample data provided. It needs to be in an sf style data frame, and needs a crs as well. Assuming you have those, but had difficulty posting them, the solution below should work. Your sf object needs to have two geometry columns, and it looks like yours should.

Using st_distance() should work, with the by_element = T argument. Examples below to either use st_distance() directly, or in dplyr::mutate to add a column for distance to the sf data frame.

library(sf)
library(tidyverse)

#### Making reproducible data
# get the nc data, make the geometry column a point with st_centroid
nc = st_read(system.file("shape/nc.shp", package="sf")) %>%
  select(NAME) %>% st_centroid()

# jitter the centroid point and add (cbind) as a second geometry column
geo2 <- st_geometry(nc) %>% st_jitter()
nc <- cbind(nc, geo2)
####

# Find the distance between the points, row-by-row
st_distance(nc$geometry,nc$geometry.1, by_element = T) %>% head()
#> Units: [m]
#> [1]  965.8162 2030.5782 1833.3081 1909.5538 1408.7908  820.0569

# or use mutate to add a column to the sf df.
nc %>% mutate(dist = st_distance(geometry, geometry.1, by_element = T))
#> Simple feature collection with 100 features and 2 fields
#> Active geometry column: geometry
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -84.05986 ymin: 34.07671 xmax: -75.8095 ymax: 36.49111
#> Geodetic CRS:  NAD27
#> First 10 features:
#>           NAME                   geometry                 geometry.1
#> 1         Ashe  POINT (-81.49823 36.4314) POINT (-81.49685 36.44001)
#> 2    Alleghany POINT (-81.12513 36.49111) POINT (-81.13681 36.47545)
#> 3        Surry POINT (-80.68573 36.41252) POINT (-80.69163 36.39673)
#> 4    Currituck POINT (-76.02719 36.40714) POINT (-76.02305 36.39029)
#> 5  Northampton POINT (-77.41046 36.42236) POINT (-77.40909 36.40974)
#> 6     Hertford POINT (-76.99472 36.36142) POINT (-76.98777 36.36623)
#> 7       Camden POINT (-76.23402 36.40122)  POINT (-76.23969 36.4181)
#> 8        Gates POINT (-76.70446 36.44428) POINT (-76.70953 36.45603)
#> 9       Warren POINT (-78.11042 36.39693) POINT (-78.11619 36.38541)
#> 10      Stokes POINT (-80.23429 36.40042) POINT (-80.24365 36.39904)
#>             dist
#> 1   965.8162 [m]
#> 2  2030.5782 [m]
#> 3  1833.3081 [m]
#> 4  1909.5538 [m]
#> 5  1408.7908 [m]
#> 6   820.0569 [m]
#> 7  1944.8192 [m]
#> 8  1382.9058 [m]
#> 9  1381.7946 [m]
#> 10  851.1106 [m]

^{Created on 2022-04-13 by the reprex package (v2.0.1)}

How to tell if a shape file (.shp) contains planar geometry or not

You are looking for this https://r-spatial.github.io/sf/reference/st_is_longlat.html EPSG:4326 is non planar, just st_transform(obj, 3857) or any other projection should work.

---

Related Topics