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]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
