inner_edges.R

## topology for simple features
## - must be added in
## - is ever needed

## rmapshaper has a very specific function to get the "inner lines" of a
## polygon mesh, but it's just the geometry no identity of which object
## or part its from
library(sf)

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

xsf <- sf::read_sf(system.file("shape/nc.shp", package="sf"))
ib <- rmapshaper::ms_innerlines(xsf)
plot(ib, col = sample(viridis::viridis(length(ib))))

## here we use the silicate PATH model and the new functions sc_segment
## with its edge classification to reconstruct an sf object where all
## individual line segments know their parent/s from the PATH model

library(tidyr)
library(silicate)
path <- PATH(xsf)
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

## arcs are all rings, not just shared
#arc <- sc_arc(path)
## segments are all instances, classified by edge groups
## and these are what we want
segs <- sc_segment(path) %>% 
        group_by(edge_) %>% mutate(n_edge = n()) %>% 
  dplyr::filter(n_edge > 1) %>% ungroup()
segs

## # A tibble: 2,128 x 6
##    .vertex0 .vertex1    path_ segment_    edge_ n_edge
##       <chr>    <chr>    <chr>    <chr>    <chr>  <int>
##  1 0e80d3b8 2963b00b e7b725bd 4bb9779e c3acfdf5      2
##  2 2963b00b 552702d3 e7b725bd f22bc478 f91cd15f      2
##  3 552702d3 74bfb0d4 e7b725bd 936dc4f2 a7e0f1cc      2
##  4 74bfb0d4 c7ad8967 e7b725bd 2dcc6bd0 7000645d      2
##  5 2177fc2e b79c6bd9 e7b725bd a28b18b7 552247ef      2
##  6 b79c6bd9 0262e8b7 e7b725bd b48a65c9 93aa65b7      2
##  7 0262e8b7 ee600ba4 e7b725bd 980fdcc1 be3e5124      2
##  8 ee600ba4 7bfe328b e7b725bd f851430c 5bc426b9      2
##  9 7bfe328b e79c8669 e7b725bd 25413041 84d5f403      2
## 10 e79c8669 e742eddc e7b725bd 662e8e60 d6ff9363      2
## # ... with 2,118 more rows

## we only want to draw them once, but knowing which of two paths we belong to
## is key
shared <- segs %>% select(path_, edge_) %>% 
  group_by(edge_) %>% 
  mutate(i = sprintf("path%i", row_number())) %>% spread(i, path_)

build_sf_from_path <- function(x) {
  g <- tibble::tibble(vertex_ = unlist(x[c(".vertex0", ".vertex1")])) %>% 
    inner_join(path$vertex, "vertex_") %>% dplyr::select(x_, y_) %>% as.matrix() %>% sf::st_linestring()
  sf::st_sf(edge = x$edge_, path1 = x$path1, path2 = x$path2, geometry = st_sfc(g))
}
## reconstruct the inner edges
d <- do.call(rbind, shared %>% inner_join(segs, "edge_") %>% 
               purrr::transpose() %>% purrr::map(build_sf_from_path))

## consider now that we know path1 and path2 for every edge, but they don't 
## necessarily identify closed rings (since there are two choices for that)
library(ggplot2)
st_crs(d) <- st_crs(xsf)
ggplot(d) + geom_sf(aes(colour = path1)) + coord_sf(datum = NA) + guides(colour = FALSE)

ggplot(d) + geom_sf(aes(colour = path2)) + coord_sf(datum = NA) + guides(colour = FALSE)

## a second plot with the original boundaries first, because it looks neat
ggplot(xsf) + geom_sf() +  geom_sf(data = d, aes(colour = path1)) + coord_sf(datum = NA) + guides(colour = FALSE)