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
SpatialLines object, containing the spatial lines representing the edgesigraph object, containing the network topologynb, containing for each vertex which edges it connets to: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)
The plot in the graph space looks like this:
plot(sln@g)
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)
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")
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)
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.