Spatial networks

A missing sp-compatible spatial data class is that for spatial networks. This document will do some simple things by combining the spatial (SpatialLines) classes of sp with the graph classes and methods from igraph. We will compute and plot the shortest path through a graph made of Delauney triangle edges from random points.

To start, we'll need to inport the spatial and graph packages:

library(sp)
library(igraph)

and define the igraph class as S4, in order to reuse it in a new S4 class

setClass("igraph")

Next, we define the SpatialLinesNetwork class, consisting of

setClass("SpatialLinesNetwork", representation(sl = "SpatialLines", g = "igraph", 
    nb = "list"), validity = function(object) {
    stopifnot(length(object@sl) == length(E(object@g)))
    stopifnot(length(object@nb) == length(V(object@g)))
})

The following section creates a SpatialLines object, and converts it into a SpatialLinesNetwork. As you can see, it only accepts “simple” objects, having only one Line per Lines set. Much of what it does is figure out the common points (igraph: vertices) where line segments (igraph: edges) meet.

SpatialLinesNetwork = function(sl) {
    stopifnot(is(sl, "SpatialLines"))
    if (!is(sl, "SpatialLinesDataFrame")) 
        sl = new("SpatialLinesDataFrame", sl, data = data.frame(id = 1:length(sl)))
    if (!all(sapply(sl@lines, length) == 1)) 
        stop("SpatialLines is not simple: each Lines element should have only a single Line")
    startEndPoints = function(x) {
        firstLast = function(L) {
            cc = coordinates(L)[[1]]
            rbind(cc[1, ], cc[nrow(cc), ])
        }
        do.call(rbind, lapply(x@lines, firstLast))
    }
    s = startEndPoints(sl)
    zd = zerodist(SpatialPoints(s))
    pts = 1:nrow(s)

    # the following can't be done vector-wise, there is a progressive effect:
    if (nrow(zd) > 0) {
        for (i in 1:nrow(zd)) pts[zd[i, 2]] = pts[zd[i, 1]]
    }
    stopifnot(identical(s, s[pts, ]))

    # map to 1:length(unique(pts))
    pts0 = match(pts, unique(pts))
    node = rep(1:length(sl), each = 2)
    nb = lapply(1:length(unique(pts)), function(x) node[which(pts0 == x)])
    g = graph(pts0, directed = FALSE)  # edges
    nodes = s[unique(pts), ]
    g$x = nodes[, 1]  # x-coordinate vertex
    g$y = nodes[, 2]  # y-coordinate vertex
    g$n = as.vector(table(pts0))  # nr of edges
    # line lengths:
    sl$length = sapply(sl@lines, function(x) LineLength(x@Lines[[1]]))
    E(g)$weight = sl$length
    # create list with vertices, starting/stopping for each edge?  add for
    # each SpatialLines, the start and stop vertex
    pts2 = matrix(pts0, ncol = 2, byrow = TRUE)
    sl$start = pts2[, 1]
    sl$end = pts2[, 2]
    new("SpatialLinesNetwork", sl = sl, g = g, nb = nb)
}
l0 = cbind(c(1, 2), c(0, 0))
l1 = cbind(c(0, 0, 0), c(0, 1, 2))
l2 = cbind(c(0, 0, 0), c(2, 3, 4))
l3 = cbind(c(0, 1, 2), c(2, 2, 3))
l4 = cbind(c(0, 1, 2), c(4, 4, 3))
l5 = cbind(c(2, 2), c(0, 3))
l6 = cbind(c(2, 3), c(3, 4))
l = list(Lines(list(Line(l0)), "e"), Lines(list(Line(l1)), "a"), Lines(list(Line(l2)), 
    "b"), Lines(list(Line(l3)), "c"), Lines(list(Line(l4)), "d"), Lines(list(Line(l5)), 
    "f"), Lines(list(Line(l6)), "g"))
sl = SpatialLines(l)
sln = SpatialLinesNetwork(sl)

We can plot these data in geographical space, using colour to denote the number of edges at each vertex:

plot(sln@g$x, sln@g$y, col = sln@g$n, pch = 16, cex = 2, asp = 1)
lines(sl)
text(sln@g$x, sln@g$y, E(sln@g), pos = 4)

plot of chunk unnamed-chunk-6

The plot in the graph space looks like this:

plot(sln@g)

plot of chunk unnamed-chunk-7

A larger example: shortest path through a Delauny triangulations network

We'll use Delauny triangulation from package deldir:

require(deldir)
## Loading required package: deldir
## deldir 0.0-22

The following function converts a set of points into a SpatialPolygons or into a SpatialLines object:

dd <- function(x, ..., to = "polygons") {
    stopifnot(is(x, "Spatial"))
    cc = coordinates(x)
    dd = deldir(list(x = cc[, 1], y = cc[, 2]), ...)
    if (to == "polygons") {
        tm = triMat(dd)
        fn = function(ix) {
            pts = tm[ix, ]
            pts = c(pts, pts[1])
            Polygons(list(Polygon(rbind(cc[pts, ]))), ix)
        }
        SpatialPolygons(lapply(1:nrow(tm), fn), proj4string = x@proj4string)
    } else if (to == "lines") {
        segs = dd$delsgs
        lst = vector("list", nrow(segs))
        for (i in 1:nrow(segs)) lst[i] = Lines(list(Line(cc[c(segs[i, 5], segs[i, 
            6]), ])), i)
        SpatialLines(lst, proj4string = x@proj4string)
    } else stop("argument to should be polygons or lines")
}

We'll generate a set of 100 points in a unit square:

set.seed(5432)
x = runif(100)
x = x[order(x)]
y = runif(100)
pts = SpatialPoints(cbind(x, y))

Next, we'll get the SpatialLines object from it:

sl = dd(pts, to = "lines")
## 
##      PLEASE NOTE:  The components "delsgs" and "summary" of the 
##      object returned by deldir() are now DATA FRAMES rather than 
##      matrices (as they were prior to release 0.0-18). 
##      See help("deldir").
##  
##      PLEASE NOTE: The process that deldir() uses for determining
##      duplicated points has changed from that used in version
##      0.0-9 of this package (and previously). See help("deldir").

and plot it:

plot(sl)

plot of chunk unnamed-chunk-12

From this object, we can create a SpatialLinesNetwork by

ln = SpatialLinesNetwork(sl)

when plotted, we can again use colour to denote the number of vertices connected to each edge:

plot(ln@g$x, ln@g$y, col = ln@g$n, pch = 16, cex = 2, asp = 1)
lines(sl, col = "grey")

plot of chunk unnamed-chunk-14

Now we compute the shortest path from the left-most point (1) to the right-most one (100):

sp = get.shortest.paths(ln@g, 1, 100)[[1]]

and plot it

plot(ln@g$x, ln@g$y, col = ln@g$n, pch = 16, cex = 1.5, asp = 1)
lines(sl, col = "grey")
points(ln@g$x[sp], ln@g$y[sp], col = "red", cex = 2)
text(ln@g$x[c(1, 100)], ln@g$y[c(1, 100)], c("start", "end"), pos = 4)

plot of chunk unnamed-chunk-16

As the edge weights are computed by Line lengths, this is the geographically shortest path.

I haven't figure out how to draw the shortest path by selecting the appropriate line segments; I guess I have to use

sp
##  [1]   1  16  18  34  44  49  63  74  90 100

which gives the node sequence, and match each pair, 1 - 16, 16 - 18 etc., to the start and end fields (which are now un-sorted!) – some challenges ahead.