rigraphr-sfspdep

Generate Neighbour List object for spatial lines in R

My goal is to generate an object of neighbourhood relationships between spatial lines in a road network. If my data were spatial polygons I could use spdep::poly2nb to do this, but I'm having trouble working out how to do this for spatial lines.

In the reprex below I try using igraph::as_adjacency_matrix create an adjacency matrix and then convert that to a neighbours list object using spdep::mat2listw. But is this the correct way to go?

Once I have the Neighours list I also want to label with the road_id attribute.

library(sfnetworks)
library(sf)
library(spdep)
library(igraph)

net <- as_sfnetwork(roxel, directed = FALSE) %>% 
  activate("edges") %>% 
  mutate(road_id = row_number()+1000)

# Make adjacency matrix
B_net <- igraph::as_adjacency_matrix(net, edges = TRUE, attr = names)

# Make neighbour list
nb_net <- mat2listw(B_net)$neighbours
# Can't use row.names in mat2listw because how do I set row.names in igraph::as_adjacency_matrix

EDIT: New approach using methof in this sfnetworks issue https://github.com/luukvdmeer/sfnetworks/issues/10 gives a neighbour list.

net_sf <- st_as_sf(net)

net_neigh <- st_touches(net_sf)

# define ad-hoc function to translate sgbp into nb (as documented in 
# https://r-spatial.github.io/spdep/articles/nb_sf.html#creating-neighbours-using-sf-objects)
as.nb.sgbp <- function(x) {
  attrs <- attributes(x)
  x <- lapply(x, function(i) { if(length(i) == 0L) 0L else i } )
  attributes(x) <- attrs
  class(x) <- "nb"
  x
}

net_nb <- as.nb.sgbp(net_neigh)
net_lw <- nb2listw(net_nb)
Solution

  • IMO there are a few approaches for creating a neighborhood list object for spatial lines similar to the output of spdep::poly2nb. Unfortunately there are also a few assumptions (analogous to the queen parameter of spdep::poly2nb) that I will try to explain here.

    First we load some packages and data

    library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sfnetworks)
library(igraph)
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union

    The first solution is based on the sf::st_touches function: 

    roxel_touches <- stplanr::geo_projected(roxel, st_touches)
roxel_touches
#> Sparse geometry binary predicate list of length 851, where the predicate was `touches'
#> first 10 elements:
#>  1: 101, 146, 493, 577, 742
#>  2: 149, 285, 294, 354, 521
#>  3: 268, 269, 270, 376, 431
#>  4: 8, 9, 276, 735
#>  5: 8, 94, 95, 108, 590
#>  6: 272, 728, 732
#>  7: 274, 275, 276, 734
#>  8: 4, 5, 95, 108, 735
#>  9: 4, 276, 281, 728
#>  10: 273, 281, 729, 733

    The output is an sgbp object that summarizes the relationships between the spatial lines according to a touches predicate. For example the first line of the output means that the first line of roxel object touches the lines 101, 146, 493, 577 and 742. The stplanr::geo_projected function is used to apply a geometric binary predicate function without reprojecting the data (otherwise we would have get some warning message).

    The touches predicate is defined in the corresponding help page (?sf::st_touches) and also here. Briefly, the sf::st_touches function matches the LINESTRING geometries that share one common point when that point lies in the union of their boundaries (i.e. the first and last point of the line). For example, the following lines intersect each other but they are not touching since the common point do not lie in their boundaries.

    es_1 <- st_sfc(
  st_linestring(rbind(c(0, -1), c(0, 1))), 
  st_linestring(rbind(c(-1, 0), c(1, 0)))
)
par(mar = rep(1, 4))
plot(es_1, col = c("red", "blue"), lwd = 2)
legend("topright", legend = c(1, 2), col = c("red", "blue"), lty = 1, lwd = 2)


    st_intersects(es_1)
#> Sparse geometry binary predicate list of length 2, where the predicate was `intersects'
#>  1: 1, 2
#>  2: 1, 2
st_touches(es_1)
#> Sparse geometry binary predicate list of length 2, where the predicate was `touches'
#>  1: (empty)
#>  2: (empty)

    Another, slightly different, solution is created using the sf::st_relate function:

    roxel_relate <- stplanr::geo_projected(roxel, st_relate, pattern = "FF*FT****")
roxel_relate
#> Sparse geometry binary predicate list of length 851, where the predicate was `relate_pattern'
#> first 10 elements:
#>  1: 101, 146, 493, 577, 742
#>  2: 149, 285, 294, 354, 521
#>  3: 268, 269, 270, 376, 431
#>  4: 8, 9, 276, 735
#>  5: 8, 94, 95, 108, 590
#>  6: 272, 732
#>  7: 274, 275, 276, 734
#>  8: 4, 5, 95, 108, 735
#>  9: 4, 276, 281, 728
#>  10: 273, 281, 729, 733

    The output is (again) an sgbp object that summarizes the relationships between the spatial lines according to a "relate" predicate with pattern "FF*FT****". That type of pattern is used to match LINESTRINGS if and only if they share at least one point in their boundaries, their interiors are not intersecting and there is no shared point between interiors and boundaries. For example the following two lines are touching but they are not related according to the pattern "FF*FT****".

    es_2 <- st_sfc(
  st_linestring(rbind(c(0, -1), c(0, 1)), "red"), 
  st_linestring(rbind(c(0, 0), c(1, 0)), "blue"), 
  crs = 32632
)
plot(es_2, col = c("red", "blue"), lwd = 2)
legend(x = "topright", legend = c(1, 2), lty = 1, col = c("red", "blue"), lwd = 2)


    st_touches(es_2)
#> Sparse geometry binary predicate list of length 2, where the predicate was `touches'
#>  1: 2
#>  2: 1
st_relate(es_2, pattern = "FF*FT****")
#> Sparse geometry binary predicate list of length 2, where the predicate was `relate_pattern'
#>  1: (empty)
#>  2: (empty)

    You can check more details about sf::st_relate function and how to build patterns in the help page of sf::st_relate. We can understand why this distinction between st_touches and st_relate is important when we compare the output of the two predicates. 

    roxel_touches
#> Sparse geometry binary predicate list of length 851, where the predicate was `touches'
#> first 10 elements:
#>  1: 101, 146, 493, 577, 742
#>  2: 149, 285, 294, 354, 521
#>  3: 268, 269, 270, 376, 431
#>  4: 8, 9, 276, 735
#>  5: 8, 94, 95, 108, 590
#>  6: 272, 728, 732
#>  7: 274, 275, 276, 734
#>  8: 4, 5, 95, 108, 735
#>  9: 4, 276, 281, 728
#>  10: 273, 281, 729, 733
roxel_relate
#> Sparse geometry binary predicate list of length 851, where the predicate was `relate_pattern'
#> first 10 elements:
#>  1: 101, 146, 493, 577, 742
#>  2: 149, 285, 294, 354, 521
#>  3: 268, 269, 270, 376, 431
#>  4: 8, 9, 276, 735
#>  5: 8, 94, 95, 108, 590
#>  6: 272, 732
#>  7: 274, 275, 276, 734
#>  8: 4, 5, 95, 108, 735
#>  9: 4, 276, 281, 728
#>  10: 273, 281, 729, 733

    We can see that they are almost identical but for a few cases (see, for example, the 6th row). We can understand the situation even better if we plot the corresponding lines:

    plot(st_geometry(roxel[c(6, 272, 728, 732), ]), col = c("red", "orange", "blue", "darkgreen"), lwd = 2)
legend(x = "topright", legend = c(6, 272, 728, 732), lty = 1, col = c("red", "orange", "blue", "darkgreen"), lwd = 2)

    We can see that the 6th line is touching all the other lines but it's not "related" to line 728 according to the pattern "FF*FT****". This implies that if you create the neighbours list starting from st_relate, then you are missing (at least) one existing link. I explicitly chose that type of pattern since the default algorithm used by sfnetworks for inferring a graph structure from an sf object is very similar (maybe identical) to that type of pattern. We can check this fact creating the line graph associated to the graph created by as_sfnetwork:

    roxel_sfnetworks <- as_sfnetwork(roxel, directed = FALSE)
roxel_line_graph <- make_line_graph(roxel_sfnetworks)
roxel_adj_list <- as_adj_list(roxel_line_graph)
str(roxel_adj_list[1:10], give.attr = FALSE)
#> List of 10
#>  $ : 'igraph.vs' int [1:5] 101 146 493 577 742
#>  $ : 'igraph.vs' int [1:5] 149 285 294 354 521
#>  $ : 'igraph.vs' int [1:5] 268 269 270 376 431
#>  $ : 'igraph.vs' int [1:4] 8 9 276 735
#>  $ : 'igraph.vs' int [1:5] 8 94 95 108 590
#>  $ : 'igraph.vs' int [1:2] 272 732
#>  $ : 'igraph.vs' int [1:4] 274 275 276 734
#>  $ : 'igraph.vs' int [1:5] 4 5 95 108 735
#>  $ : 'igraph.vs' int [1:4] 4 276 281 728
#>  $ : 'igraph.vs' int [1:4] 273 281 729 733

    and test their equality: 

    all(unlist(mapply(function(x, y) unique(x) == unique(y), roxel_adj_list, roxel_relate)))
#> [1] TRUE

    I had to use the unique function since the igraph output may return duplicated indexes when there are multipled edges between two vertices, i.e.

    is.simple(roxel_sfnetworks)
#> [1] FALSE

    Anyway, we can see that the output is identical to roxel_relate, which means that we cannot always safely use st_relate or make_line_graph functions without particular considerations on the underlying street network object. I don't want to add too many details here since the examples are too complicated, but we cannot either always use the st_touches function for inferring the neighbour list object (but if you want to read a little bit more on this topic I recently wrote a paper on this topic, which is currently under review).

    Anyway, the summary of that paper is that I think we can safely use st_touches or st_relate functions to generate the neighbour list only after transforming the roxel object using the stplanr::rnet_breakup_vertices function. (Actually there are still some really minor differences between the two approaches, but I still have to properly figure out how to solve them. For the moment I would use st_touches).

    # 1 - Transform the input roxel object
roxel2 <- stplanr::rnet_breakup_vertices(roxel)
#> Splitting rnet object at the intersection points between nodes and internal vertexes

# 2- Apply st_touches function
roxel2_touches <- stplanr::geo_projected(roxel2, st_touches)
roxel2_touches
#> Sparse geometry binary predicate list of length 876, where the predicate was `touches'
#> first 10 elements:
#>  1: 105, 150, 504, 591, 765
#>  2: 153, 291, 300, 363, 533
#>  3: 274, 275, 276, 386, 441
#>  4: 8, 9, 282, 758
#>  5: 8, 98, 99, 112, 604
#>  6: 278, 750, 751, 755
#>  7: 280, 281, 282, 757
#>  8: 4, 5, 99, 112, 758
#>  9: 4, 282, 287, 750
#>  10: 279, 287, 752, 756

    Then, if you want, you can transform the result into an nb object: 

    as.nb.sgbp <- function(x) {
  attrs <- attributes(x)
  x <- lapply(x, function(i) { if(length(i) == 0L) 0L else i } )
  attributes(x) <- attrs
  class(x) <- "nb"
  x
}

roxel2_nb <- as.nb.sgbp(roxel2_touches)
roxel2_lw <- spdep::nb2listw(roxel2_nb)
roxel2_lw
#> Characteristics of weights list object:
#> Neighbour list object:
#> Number of regions: 876 
#> Number of nonzero links: 3318 
#> Percentage nonzero weights: 0.4323826 
#> Average number of links: 3.787671 
#> 
#> Weights style: W 
#> Weights constants summary:
#>     n     nn  S0       S1       S2
#> W 876 767376 876 494.4606 3580.933

    We can also use the "line-graph" option: 

    str(as_adj_list(make_line_graph(as_sfnetwork(roxel2, directed = FALSE)))[1:10], give.attr = FALSE)
#> List of 10
#>  $ : 'igraph.vs' int [1:5] 105 150 504 591 765
#>  $ : 'igraph.vs' int [1:5] 153 291 300 363 533
#>  $ : 'igraph.vs' int [1:5] 274 275 276 386 441
#>  $ : 'igraph.vs' int [1:4] 8 9 282 758
#>  $ : 'igraph.vs' int [1:5] 8 98 99 112 604
#>  $ : 'igraph.vs' int [1:4] 278 750 751 755
#>  $ : 'igraph.vs' int [1:4] 280 281 282 757
#>  $ : 'igraph.vs' int [1:5] 4 5 99 112 758
#>  $ : 'igraph.vs' int [1:4] 4 282 287 750
#>  $ : 'igraph.vs' int [1:4] 279 287 752 756

    Created on 2020-06-01 by the reprex package (v0.3.0)

    I know that I wrote a lot of (maybe unnecessary) information but I think they are useful to properly give context about the suggested solution. If it's not clear please add more questions.