-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathnw.Rmd
178 lines (171 loc) · 6.28 KB
/
nw.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
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:
```{r}
library(sp)
library(igraph)
```
and define the igraph class as S4, in order to reuse it in a new S4 class
```{r}
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:
```{r}
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.
```{r}
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)
}
```
```{r}
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:
```{r fig.width=7, fig.height=6}
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:
```{r}
plot(sln@g)
```
A larger example: shortest path through a Delauny triangulations network
------------------------------------------
We'll use Delauny triangulation from package `deldir`:
```{r}
require(deldir)
```
The following function converts a set of points into a `SpatialPolygons` or into a `SpatialLines` object:
```{r}
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:
```{r}
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:
```{r}
sl = dd(pts, to = "lines")
```
and plot it:
```{r}
plot(sl)
```
From this object, we can create a `SpatialLinesNetwork` by
```{r}
ln = SpatialLinesNetwork(sl)
```
when plotted, we can again use colour to denote the number of vertices connected to each edge:
```{r fig.width=7, fig.height=6}
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):
```{r}
sp = get.shortest.paths(ln@g, 1, 100)[[1]]
```
and plot it
```{r fig.width=7, fig.height=7}
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
```{r}
sp
```
which gives the node sequence, and match each pair, `r sp[1]` - `r sp[2]`, `r sp[2]` - `r sp[3]` etc., to the start and end fields (which are now un-sorted!) -- some challenges ahead.