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

- a
`SpatialLines object`

, containing the spatial lines representing the edges - an
`igraph`

object, containing the network topology - a neighbourhood list
`nb`

, 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.