Problems with spatial and graph operations on a Spatial Lines Network

[agila5]

This issue shows the problems with the current implementations of spatial and graph operations on a network.

The problem is illustrated using code from stplanr, the challenge will be to solve this problem in future versions of sfnetworks.

We start this example by creating a spatial lines network object with the function SpatialLinesNetwork of the R package stplanr.

# packages
library(sf)
library(stplanr)
library(tmap)
library(leaflet)
library(ggspatial)
library(tidygraph)

sln_sf <- SpatialLinesNetwork(route_network_sf)

First of all, it should be noted that, at the moment, the plotting of spatial lines network objects with ggplot2/tmap/leaflet is unnecessary complicated. For example, all the following approaches fail.

tm_shape(sln_sf) +
  tm_lines()
#> Error in get_projection(mshp_raw, output = "crs"): shp is neither a sf, sp, nor a raster object
leaflet(sln_sf) %>%
  addTiles() %>%
  addPolylines()
#> Error in polygonData.default(data): Don't know how to get path data from object of class sfNetwork
ggplot() +
  annotation_map_tile() +
  geom_sf(data = sln_sf, size = 0.75)
#> Error: `data` must be a data frame, or other object coercible by `fortify()`, not an S4 object with class sfNetwork

As we can see from the error messages, the user first has to extract the sf data and then plot it. For example,

ggplot() +
  annotation_map_tile() +
  geom_sf(data = sln_sf@sl, size = 0.75)
#> Zoom: 13

The same idea can be applied to igraph functions, like

igraph::edge_betweenness(sln_sf)
#> Error in igraph::edge_betweenness(sln_sf): Not a graph object
igraph::is_connected(sln_sf)
#> Error in igraph::is_connected(sln_sf): Not a graph object

which work only if we extract the graph component from the spatial network object,

igraph::edge_betweenness(sln_sf@g)
#>  [1]  51 100   0   0 147 108  84 129  56 116 151   0 220  16  49 102 121
#> [18] 182 100 161  95  57  48   6  37 135 375 386   0 289 106  51 135  10
#> [35] 245  43   4 184 200  97   0 162 121   0  78  33 198  51  57 120  77
#> [52] 117 150  85  51 127 118   0  79  44 135  98  34   0 119  37 101 120
#> [69]  99  25  41  17 153  36 226 151  40 164  20  11
igraph::is_connected(sln_sf@g)
#> [1] TRUE

This separation between the spatial and the graph components of a spatial network is counterintuituve and complicated.

Moreover, we cannot apply to that network object any spatial or graph function that implies a modification of the spatial or graph structure. For example, if we create a buffer of 400m around Chapteltown, which is one of the neighborhoods in the center of the network object,

chapeltown_geocode <- tmaptools::geocode_OSM("chapeltown leeds")
chapteltown_buffer <- geo_projected(
  st_sf(st_sfc(st_point(c(chapeltown_geocode$coords[1], chapeltown_geocode$coords[2])), crs = 4326)),
  st_buffer,
  dist = 400
)

and filter all the edges inside the buffer, we get an error:

sln_sf[chapteltown_buffer, ]
#> Error in sln_sf[chapteltown_buffer, ]: object of type 'S4' is not subsettable

This example is even more problematic than the others, since, in this case, we cannot simply work on the sf component of the spatial lines network object. In fact, even if we modify it,

sln_sf2 <- sln_sf
sln_sf2@sl <- sln_sf2@sl[chapteltown_buffer, ]
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

this does not change the graph component since they are not strictly linked:

nrow(sln_sf2@sl)
#> [1] 9
igraph::E(sln_sf@g)
#> + 80/80 edges from 6e154b9:
#>  [1]  1-- 2  1-- 3  1-- 4  3-- 5  3-- 6  6-- 7  6-- 8  5-- 8  8-- 9  5--10
#> [11]  7-- 9  7--11  9--12  9--11 10--13 10--14 13--14 12--14 14--15 13--16
#> [21] 16--17  4--16  4--18 18--19 17--18 17--20 12--21 21--22 21--22 22--23
#> [31] 22--24 19--20 15--20 19--25 15--26 15--25 25--27 23--28 23--29 24--29
#> [41] 24--30 24--30 30--31 31--32 31--32 32--33 29--34 11--33 33--35 28--36
#> [51] 28--37 34--38 34--39 37--38 38--40 37--41 41--42 41--43 36--44 26--36
#> [61] 39--45 35--39 39--46 35--46 45--47 42--45 43--44 42--43 47--48 47--49
#> [71] 48--49 46--49 44--50 26--27 26--50 48--51 50--52 50--51 51--52 27--52
igraph::E(sln_sf2@g)
#> + 80/80 edges from 6e154b9:
#>  [1]  1-- 2  1-- 3  1-- 4  3-- 5  3-- 6  6-- 7  6-- 8  5-- 8  8-- 9  5--10
#> [11]  7-- 9  7--11  9--12  9--11 10--13 10--14 13--14 12--14 14--15 13--16
#> [21] 16--17  4--16  4--18 18--19 17--18 17--20 12--21 21--22 21--22 22--23
#> [31] 22--24 19--20 15--20 19--25 15--26 15--25 25--27 23--28 23--29 24--29
#> [41] 24--30 24--30 30--31 31--32 31--32 32--33 29--34 11--33 33--35 28--36
#> [51] 28--37 34--38 34--39 37--38 38--40 37--41 41--42 41--43 36--44 26--36
#> [61] 39--45 35--39 39--46 35--46 45--47 42--45 43--44 42--43 47--48 47--49
#> [71] 48--49 46--49 44--50 26--27 26--50 48--51 50--52 50--51 51--52 27--52

This means that, after the spatial filtering, some edges in the graph component are linked to LINESTRINGS that do no exist in the spatial component. The same will happen if we modify the graph component. The only procedure that works is to filter the original sf data and, then, rebuild the network from scratch, i.e.

sln_sf3 <- SpatialLinesNetwork(route_network_sf[chapteltown_buffer, ])
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

nrow(sln_sf3@sl)
#> [1] 9
igraph::E(sln_sf3@g)
#> + 9/9 edges from 700d8f9:
#> [1] 1--2 1--3 1--4 4--5 4--5 5--6 5--7 6--8 6--9

^{Created on 2019-10-13 by the reprex package (v0.3.0)}

[luukvdmeer]

The "sfnetworks answer" to close this issue:

library(sfnetworks)
library(igraph)
library(tidygraph)
library(sf)
library(ggplot2)

net = as_sfnetwork(roxel, directed = FALSE)

Regarding ggplot, you will still need to extract nodes and edges as sf objects. However, this is fairly easy now since you can do that without activating first, and instead provide "nodes" or "edges" as argument to st_as_sf.

st_as_sf(net, "nodes")
#> Simple feature collection with 701 features and 0 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 7.522622 ymin: 51.94151 xmax: 7.546705 ymax: 51.9612
#> CRS:            EPSG:4326
#> # A tibble: 701 x 1
#>               geometry
#>            <POINT [°]>
#>  1 (7.533722 51.95556)
#>  2 (7.533461 51.95576)
#>  3 (7.532442 51.95422)
#>  4  (7.53209 51.95328)
#>  5 (7.532709 51.95209)
#>  6 (7.532869 51.95257)
#>  7 (7.540063 51.94468)
#>  8  (7.53822 51.94546)
#>  9  (7.537673 51.9475)
#> 10 (7.537614 51.94562)
#> # … with 691 more rows

st_as_sf(net, "edges")
#> Simple feature collection with 851 features and 4 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 7.522594 ymin: 51.94151 xmax: 7.546705 ymax: 51.9612
#> CRS:            EPSG:4326
#> # A tibble: 851 x 5
#>     from    to name         type                                        geometry
#>    <int> <int> <fct>        <fct>                               <LINESTRING [°]>
#>  1     1     2 Havixbecker… reside…       (7.533722 51.95556, 7.533461 51.95576)
#>  2     3     4 Pienersallee second… (7.532442 51.95422, 7.53236 51.95377, 7.532…
#>  3     5     6 Schulte-Ber… reside… (7.532709 51.95209, 7.532823 51.95239, 7.53…
#>  4     7     8 <NA>         path    (7.540063 51.94468, 7.539696 51.94479, 7.53…
#>  5     9    10 Welsingheide reside…        (7.537673 51.9475, 7.537614 51.94562)
#>  6    11    12 <NA>         footway (7.543791 51.94733, 7.54369 51.94686, 7.543…
#>  7    13    14 <NA>         footway        (7.54012 51.94478, 7.539931 51.94514)
#>  8     8    10 <NA>         path    (7.53822 51.94546, 7.538131 51.94549, 7.538…
#>  9     7    15 <NA>         track   (7.540063 51.94468, 7.540338 51.94468, 7.54…
#> 10    16    17 <NA>         track   (7.5424 51.94599, 7.54205 51.94629, 7.54196…
#> # … with 841 more rows

Then plotting gets as easy as.

ggplot() +
  geom_sf(data = st_as_sf(net, "edges")) +
  geom_sf(data = st_as_sf(net, "nodes"))

The algorithms from igraph can be directly applied to the network.

head(igraph::edge_betweenness(net))
#> [1] 10406.362 12118.803  6837.847  4467.540 32124.052   674.000
igraph::is_connected(net)
#> [1] FALSE

Spatial filtering can be done with st_filter.

p1 = st_point(c(7.53173, 51.95662))
p2 = st_point(c(7.53173, 51.95190))
p3 = st_point(c(7.53778, 51.95190))
p4 = st_point(c(7.53778, 51.95662))

rect = st_multipoint(c(p1, p2, p3, p4)) %>% 
  st_cast('POLYGON') %>% 
  st_sfc(crs = 4326)

net_f = net %>%
  st_filter(rect, .pred = st_intersects)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

par(mar = c(1, 1, 1, 1))
plot(net, col = "grey")
plot(rect, border = "Red", lwd = 2, add = TRUE)
plot(net_f, add = TRUE)

^{Created on 2020-06-16 by the reprex package (v0.3.0)}