## https://CRAN.R-project.org/package=spdep/vignettes/nb_igraph.html
## this function uniquifies the segments, very much WIP
u_edges <- function(x, ...) UseMethod("u_edges")
u_edges.PRIMITIVE <- function(x, ...) {
u_edges(x[["segment"]])
}
u_edges.data.frame <- function(x, ...) {
u2 <- x %>% mutate(uu = paste(pmin(.vertex0, .vertex1), pmax(.vertex0, .vertex1), sep = "_"))
#dplyr::distinct(select(u2, uu, segment_), uu, .keep_all = TRUE)
select(u2, uu, segment_)
}
find neighbours
Provides a nested tibble of row-indexes for both vertices and edges. No topological fixing or fudging is done here.
find_nb <- function(x) {
etab <- vtab <- vector("list", nrow(x$object))
for (ith in seq_along(etab)) {
xith <- x$object[ith, "object_"] %>% inner_join(x$path, "object_") %>%
inner_join(x$path_link_vertex, "path_") %>%
inner_join(x$vertex, "vertex_")
## join by vertex
idx_vertex <- xith %>% dplyr::select(vertex_) %>%
inner_join(x$path_link_vertex, "vertex_") %>%
inner_join(x$path, "path_") %>%
distinct(object_) ##%>% inner_join(nc %>% mutate(object_ = x$object_))
v_idx <- match(idx_vertex$object_, x$object$object_)
## join the subset data to the main on unique segment
## (unique as in order of vertices is irrelevant)
idx_edge <- u_edges.data.frame(x$object[ith, "object_"] %>%
inner_join(x$path, "object_") %>%
inner_join(x$segment, "path_")) %>%
select(-segment_) %>% inner_join(u_edges.PRIMITIVE(x), "uu") %>%
inner_join(x$segment, "segment_") %>% inner_join(x$path, "path_") %>%
distinct(object_)
e_idx <- match(idx_edge$object_, x$object$object_)
## remove the self link
v_idx <- setdiff( v_idx, ith)
e_idx <- setdiff(e_idx, ith)
vtab[[ith]] <- tibble(neighbour = v_idx, object_ = ith)
etab[[ith]] <- tibble(neighbour = e_idx, object_ = ith)
}
tibble(vnb = vtab, enb = etab)
}
build edge-based igraph from nested tibbles (see find_nb)
g_from_nb <- function(x, layout = NULL) {
g <- as_tibble(x) %>% dplyr::select(enb) %>% unnest() %>% select(object_, neighbour) %>% graph_from_data_frame()
if (!is.null(layout)) {
igraph::V(g)$x <- layout[[1]]
igraph::V(g)$y <- layout[[2]]
}
g
}
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
library(tibble)
library(igraph)
##
## Attaching package: 'igraph'
## The following object is masked from 'package:tibble':
##
## as_data_frame
## The following objects are masked from 'package:dplyr':
##
## %>%, as_data_frame, groups, union
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
library(scsf)
## Loading required package: sc
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:igraph':
##
## %>%, crossing
pp <- PRIMITIVE(minimal_mesh)
library(sf)
## Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3
##
## Attaching package: 'sf'
## The following object is masked from 'package:igraph':
##
## %>%
columbus <- sf::read_sf(system.file("etc/shapes/columbus.shp", package="spdep"))
nb_q <- spdep::poly2nb(as(columbus, "Spatial"))
str(nb_q)
## List of 49
## $ : int [1:2] 2 3
## $ : int [1:3] 1 3 4
## $ : int [1:4] 1 2 4 5
## $ : int [1:4] 2 3 5 8
## $ : int [1:8] 3 4 6 8 9 11 15 16
## $ : int [1:2] 5 9
## $ : int [1:4] 8 12 13 14
## $ : int [1:6] 4 5 7 11 12 13
## $ : int [1:8] 5 6 10 15 20 22 25 26
## $ : int [1:4] 9 17 20 22
## $ : int [1:5] 5 8 12 15 16
## $ : int [1:6] 7 8 11 13 14 16
## $ : int [1:4] 7 8 12 14
## $ : int [1:6] 7 12 13 16 18 19
## $ : int [1:6] 5 9 11 16 25 26
## $ : int [1:8] 5 11 12 14 15 18 24 25
## $ : int [1:3] 10 20 23
## $ : int [1:4] 14 16 19 24
## $ : int [1:3] 14 18 24
## $ : int [1:10] 9 10 17 22 23 27 32 33 35 40
## $ : int [1:3] 24 30 34
## $ : int [1:6] 9 10 20 26 27 28
## $ : int [1:3] 17 20 32
## $ : int [1:7] 16 18 19 21 25 29 30
## $ : int [1:8] 9 15 16 24 26 28 29 30
## $ : int [1:6] 9 15 22 25 28 29
## $ : int [1:4] 20 22 28 33
## $ : int [1:9] 22 25 26 27 29 33 35 37 38
## $ : int [1:7] 24 25 26 28 30 37 38
## $ : int [1:5] 21 24 25 29 37
## $ : int [1:3] 34 36 39
## $ : int [1:4] 20 23 40 41
## $ : int [1:4] 20 27 28 35
## $ : int [1:4] 21 31 36 42
## $ : int [1:7] 20 28 33 38 40 43 44
## $ : int [1:5] 31 34 39 42 46
## $ : int [1:6] 28 29 30 38 43 45
## $ : int [1:6] 28 29 35 37 43 44
## $ : int [1:3] 31 36 46
## $ : int [1:5] 20 32 35 41 47
## $ : int [1:3] 32 40 47
## $ : int [1:2] 34 36
## $ : int [1:6] 35 37 38 44 45 48
## $ : int [1:5] 35 38 43 48 49
## $ : int [1:4] 37 43 48 49
## $ : int [1:2] 36 39
## $ : int [1:2] 40 41
## $ : int [1:4] 43 44 45 49
## $ : int [1:3] 44 45 48
## - attr(*, "class")= chr "nb"
## - attr(*, "region.id")= chr [1:49] "ID1" "ID2" "ID3" "ID4" ...
## - attr(*, "call")= language spdep::poly2nb(pl = as(columbus, "Spatial"))
## - attr(*, "type")= chr "queen"
## - attr(*, "sym")= logi TRUE
prim <- PRIMITIVE(columbus)
d <- bind_cols(columbus, find_nb(prim) )
## the tidy print method
as_tibble(d) %>% select(geometry, vnb, enb)
## # A tibble: 49 x 3
## geometry vnb enb
## <simple_feature> <list> <list>
## 1 <POLYGON((8.6...> <tibble [2 x 2]> <tibble [2 x 2]>
## 2 <POLYGON((8.2...> <tibble [3 x 2]> <tibble [3 x 2]>
## 3 <POLYGON((8.6...> <tibble [4 x 2]> <tibble [4 x 2]>
## 4 <POLYGON((8.4...> <tibble [4 x 2]> <tibble [4 x 2]>
## 5 <POLYGON((8.6...> <tibble [8 x 2]> <tibble [7 x 2]>
## 6 <POLYGON((9.4...> <tibble [2 x 2]> <tibble [2 x 2]>
## 7 <POLYGON((8.0...> <tibble [4 x 2]> <tibble [3 x 2]>
## 8 <POLYGON((8.2...> <tibble [6 x 2]> <tibble [5 x 2]>
## 9 <POLYGON((9.3...> <tibble [8 x 2]> <tibble [6 x 2]>
## 10 <POLYGON((10....> <tibble [4 x 2]> <tibble [3 x 2]>
## # ... with 39 more rows
## the coordinates of a centroid of each polygon
lyout <- as.data.frame(do.call(rbind, lapply(st_geometry(st_centroid(columbus)), function(x) unclass(x))))
as_tibble(d) %>% g_from_nb(layout = lyout) %>% plot(vertex.size = 0, edge.arrow.size = 0, asp = 0)
plot(lyout, type = "n", axes = FALSE, xlab = "", ylab = "", asp = "")
plot(nb_q, lyout, col="grey", add = TRUE)