[R-sig-Geo] self-intersecting polygons

[Dexter Locke]

Dear list,

I have a set of 15 and 30 minute drive-time polygons created
from mapboxapi::mb_isochrone. Ultimately I want to erase the smaller
polygons from the larger so that the polygons are non-overlapping,
concentric rings.

The whole dataset includes ~1600 polygons for each drive-time spread out
over the 50 US states plus Washington DC and Puerto Rico. 13 of them have
apparent geometry errors, but the example below includes just 4 polygons
with one self intersection in Hawaii.

script (also below) and data available here:
http://dexterlocke.com/wp-content/uploads/2021/11/reprex_files.zip (I
couldn't figure out how to get the script to load data into R from a URL)

library(sf)
library(tidyverse)

# read in data
(small <- st_read(paste0(getwd(), '/data/reprex_files/small_polys.shp')))
(large <- st_read(paste0(getwd(), '/data/reprex_files/large_polys.shp')))

# look on basemap (optional)
mapview::mapview(small, color = 'black') + mapview::mapview(large, color =
'red')

# create function to erase holes to make 'donuts'
# https://github.com/r-spatial/sf/issues/346
st_erase <- function(x, y) st_difference(x, st_union(st_combine(y)))

# apply the erasing function
donut <- large %>% st_erase(., small) # fails

# 1st option explored, turn off spherical geometry as suggested here:
# https://github.com/r-spatial/sf/issues/1762

# here
https://gis.stackexchange.com/questions/404385/r-sf-some-edges-are-crossing-in-a-multipolygon-how-to-make-it-valid-when-using/404454

# and here
#
https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data/68481205#68481205

sf_use_s2(FALSE)

donut <- large %>% st_erase(., small) # results in errors..

# are the polygons valid? Yes
small %>% st_is_valid()
large %>% st_is_valid()

# 2nd option explored, more intensive geometry cleaning, as suggested here:
#
https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data/68481205#68481205
and elsewhere
small$geometry <- small$geometry %>%
  s2::s2_rebuild() %>%
  sf::st_as_sfc()

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

# re-attempt
donut <- large %>% st_erase(., small) # results in errors

# 3rd option explored, zero-distance buffer to help reset geometries
small <- small %>% st_buffer(., 0)
large <- large %>% st_buffer(., 0)

# re-attempt
donut <- large %>% st_erase(., small) # more self-intersection errors

# 4th option explored, reproject before buffering
# suggested here:
https://gis.stackexchange.com/questions/404385/r-sf-some-edges-are-crossing-in-a-multipolygon-how-to-make-it-valid-when-using/404388#404388
small <- small %>% st_transform(crs = 5070) %>% st_buffer(0) %>%
st_transform(crs = 4269)
large <- large %>% st_transform(crs = 5070) %>% st_buffer(0) %>%
st_transform(crs = 4269)

# re-attempt
donut <- large %>% st_erase(., small) # more self-intersection errors

# 5th option
small_s2 <- small %>%
  s2::s2_rebuild(.,
                 options = s2::s2_options(edge_type = "undirected",
                 duplicate_edges = FALSE,
                 split_crossing_edges = TRUE,
                 validate = TRUE)) %>%
  st_as_sf()

large_s2 <- large %>%
  s2::s2_rebuild(.,
                 options = s2::s2_options(edge_type = "undirected",
                 duplicate_edges = FALSE,
                 split_crossing_edges = TRUE,
                 validate = TRUE)) %>%
  st_as_sf()

# re-attempt
donut <- large_s2 %>% st_erase(., small_s2) # more self-intersection errors


# 6th option explored using rgeos and sp, as suggested here:
#
https://gis.stackexchange.com/questions/163445/getting-topologyexception-input-geom-1-is-invalid-which-is-due-to-self-intersec/163480#163480

# same problems..

# so where is the double-vertex?
large %>%
  filter(id == 1017) %>% # based on zooming and panning around via mapview
  st_cast('POINT') %>%
  mutate(geom_text = as.character(geometry)) %>%
  group_by(geom_text) %>%
  count() %>%
  arrange(desc(n))

# as indicated from some errors and from query directly above
bad_point <- tibble(row = 1, x = -155.056700377001,  y = 19.5963974859705)
%>%
  st_as_sf(., coords = c('x', 'y')) %>%
  st_set_crs(4269)

# take a look
large %>% mapview::mapview() + mapview::mapview(bad_point, color = 'red')

# Ok, so there are two vertices there.. how can I remove the duplicates and
get on with the rest of the analyses?

Thank you, Dexter

	[[alternative HTML version deleted]]