general_inner_ring_rewinder.R

library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3

library(rnaturalearth)
##devtools::install_github("mdsumner/scsf")
library(scsf)

## Loading required package: sc

library(sc)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

spdf_france_country <- ne_states(country = 'france') %>% head(-5)
set.seed(4)
d2 <-  st_as_sf(spdf_france_country) %>% sample_frac(0.39)
plot(d2[1])

inner_ring_touching <- st_sf(geometry  = st_sfc(st_multipolygon(unlist(lapply(st_geometry(d2), unclass), recursive = FALSE))), name = "internal boundaried abomination")

prim <- PRIMITIVE(inner_ring_touching)

## which segments have to go? find where interchanging vertex pairs is not unique
prim$segment$useg <-  apply(cbind(prim$segment$.vertex0, prim$segment$.vertex1), 1, function(x) paste(sort(x), collapse = "-"))
prim$segment$bad <- duplicated(prim$segment$useg)



### aa is a matrix of vertex id
## this is slow and a bit bogus, note that there's a problem with 
## multiple matches - but they should be covered by filter out those
## segments in the first place as per "useg" anti join
cycles <- function(aa) {
  ii <- 1
  set0 <- ii
  visited <- logical(nrow(aa))
  while(!all(visited)) {
    i0 <- ii
    repeat {
      ii <- which(aa[,1] == aa[ii, 2])
      if (length(ii) < 1 | ii[1] == i0) {
        set0 <- c(set0, NA_integer_)
        break; 
      }
      set0 <- c(set0, ii)
    }
    visited <- seq(nrow(aa)) %in% na.omit(set0)
    ii <- which(!visited)[1L]
    if (!is.na(ii)) set0 <- c(set0, ii)
  }
  l <- split(set0, c(0, cumsum(abs(diff(is.na(set0))))))
  bind_rows(lapply(l[!unlist(lapply(l, function(x) all(is.na(x))))], function(x) tibble(row = x)), .id = "cycle")
}

## these are the ids of the vertex pairs to be kept, one row per segment
#bb <- prim$segment %>% dp %>% dplyr::select(.vertex0, .vertex1) %>% as.matrix()
## I'm failing to get anti_join to work today for some reason
anti <- prim$segment %>% filter(bad)  %>% inner_join(prim$segment %>% select(useg))

## Joining, by = "useg"

segs <- prim$segment %>% filter( !prim$segment$useg %in% anti$useg) 
bb <- prim$segment %>% filter( !prim$segment$useg %in% anti$useg) %>% dplyr::select(.vertex0, .vertex1) %>% as.matrix() 
## these are the indices from bb of each valid
ind <- cycles(bb)

#segs <- prim$segment %>% dplyr::filter(!bad) 

new_paths_list <- lapply(split(ind$row, ind$cycle), function(index) tibble(vertex_path = segs[index, c(".vertex0", ".vertex1")] %>% as.matrix() %>% t() %>% as.vector()))
plot(inner_ring_touching)

new_paths <- bind_rows(lapply(new_paths_list, 
                              function(vertex_tab) prim$vertex[match(vertex_tab$vertex_path[!duplicated(vertex_tab$vertex_path)], prim$vertex$vertex_), c("x_", "y_")]), .id = "newpath")
library(ggplot2)
#ggplot(new_paths, aes(x = x_, y = y_, group = newpath)) + geom_polygon()

library(ggpolypath)
ggplot(new_paths, aes(x = x_, y = y_, group = newpath, fill = newpath)) + ggpolypath::geom_polypath()

## what about the crazy buffer kludge? 

fix <- st_buffer(inner_ring_touching, dist = 0)
st_is_valid(inner_ring_touching)

## Warning in evalq((function (..., call. = TRUE, immediate. = FALSE,
## noBreaks. = FALSE, : Self-intersection at or near point -2.4610974491818638
## 47.452644464352034

## [1] FALSE

st_is_valid(fix)

## [1] TRUE

## it just amazes me how powerful this graph decomposition is, but
## seemingly no understanding of how close to our reach it is
plot(fix)