Skip to content

Instantly share code, notes, and snippets.

@aquilo
Created November 21, 2017 22:48
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save aquilo/c2a8c412e9394bf384a9199539131415 to your computer and use it in GitHub Desktop.
Save aquilo/c2a8c412e9394bf384a9199539131415 to your computer and use it in GitHub Desktop.
Reconstruction of Switzerland using swisscom open data (trips between cantons)
kt pop
ZH 1466.4
BE 1017.5
LU 398.8
UR 36.0
SZ 154.1
OW 37.1
NW 42.4
GL 40.0
ZG 122.1
FR 307.5
SO 266.4
BS 191.8
BL 283.2
SH 79.8
AR 54.5
AI 16.0
SG 499.1
GR 196.6
AG 653.7
TG 267.4
TI 351.9
VD 773.4
VS 335.7
NE 178.1
GE 484.7
JU 72.8
library(plyr)
# read swisscom open data from
# https://opendata.swisscom.com/explore/dataset/trips-der-reisenden-zwischen-schweizer-kantonen/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true
# already downloaded:
m0 = read.delim("trips-der-reisenden-zwischen-schweizer-kantonen.csv", header = TRUE, sep = ";", quote = "\"",
dec = ".", fill = TRUE, comment.char = "")
# remove unnecessary connections
m0 = subset(m0, Ziel.Kanton != 'Undefined')
m0 = subset(m0, Start.Kanton != 'Undefined')
m0 = subset(m0, Start.Kanton != Ziel.Kanton)
# m0 = m0[m0$Start.Kanton != 'Undefined',]
m0$Start.Kanton = factor(m0$Start.Kanton)
m0$Ziel.Kanton = factor(m0$Ziel.Kanton)
m = m0
kt = c("Aargau" = "AG", "Appenzell Ausserrhoden" = "AR", "Appenzell Innerrhoden" = "AI","Basel-Landschaft" = "BL", "Basel-Stadt" = "BS",
"Bern" = "BE", "Fribourg" = "FR", "Geneva" = "GE", "Grisons" = "GR", "Jura" = "JU", "Lucerne" = "LU", "Neuchâtel" = "NE",
"Nidwalden" = "NW", "Obwalden" = "OW", "Schaffhausen" = "SH", "Schwyz" = "SZ", "Solothurn" = "SO", "St. Gallen" = "SG",
"Thurgau" = "TG", "Ticino" = "TI", "Uri" = "UR", "Valais" = "VS", "Vaud" = "VD", "Zug" = "ZG", "Zurich" = "ZH")
m$kt.start <- revalue(m$Start.Kanton, kt)
m$kt.ziel <- revalue(m$Ziel, kt)
m = subset(m, select = c("kt.start", "kt.ziel", "Anzahl.Reisende"))
# simple way to make the distance matrix symmetric
m2 = m
names(m2) = c("kt.ziel", "kt.start", "Anzahl.Reisende")
m = rbind(m,m2)
# we are just interested in the sums and not individual dates
ma = aggregate(m$Anzahl.Reisende, list(m$kt.start, m$kt.ziel), sum)
names(ma) = c("kt.start", "kt.ziel", "Anzahl.Reisende")
# population 2015 originally copied table at
# https://www.bfs.admin.ch/bfsstatic/dam/assets/1922812/master
pop = read.delim("pop.txt", header = TRUE, sep = "\t", quote = "\"",
dec = ".", fill = TRUE, comment.char = "")
ma = merge (ma, pop, by.x="kt.start", by.y="kt")
ma = rename(ma, c("pop" = "pop.start"))
ma = merge (ma, pop, by.x="kt.ziel", by.y="kt")
ma = rename(ma, c("pop" = "pop.ziel"))
# gravitation model: F = (P1 * P2) / (distance^2)
# therefore dist:
ma$dist = sqrt(ma$pop.start * ma$pop.ziel / ma$Anzahl.Reisende)
# matrix
ma = with(ma, tapply(dist, list(kt.start, kt.ziel),"[[",1))
d = ma
d[is.na(d)] = 0
fit <- isoMDS(d, k=2)
fit
# plot solution
x <- -fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="", ylab="",
main="Nonmetric MDS", type="n", asp=1, xaxt='n', yaxt='n', ann=FALSE)
text(x, y, labels = row.names(fit$points), cex=1.1)
@aquilo
Copy link
Author

aquilo commented Nov 21, 2017

switzerland

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment