Schlagwort-Archive: R

Schoenere Karten mit R und ggplot2

Vor einiger Zeit hatte ich hier schon einmal beschrieben, wie sich mit R einfache Karten und Visualisierungen herstellen lassen. Das geht natuerlich noch deutlich komplexer — aber eben auch schoener. Edward Tufte sollte eigentlich allen ein Begriff sein, die sich mit Visualisierung und grafischer Aufbereitung von Information beschaeftigen; „The Visual Display of Quantitative Information“ empfehle ich immer wieder gerne (Amazon-Affilliate-Link; alternativ ISBN 0961392142).

Neben Tufte gibt es noch ein weiteres grosses Werk, auf das ich bei der ganzen Kartiererei in R gestossen bin: „The Grammar of Graphics“ (Amazon-Affilliate-Link, alternativ ISBN 0387245448). „Gestossen“ ist hier weit ausgelegt, ich habe das Buch nicht einmal gelesen. Ich kenne aber die auf ihm basierende R-Library ggplot2, die die Grammar of Graphics fuer ihre Plots verwendet — und sie dabei aesthetisch ansehlich werden laesst.

Ich verwende im folgenden Beispiel wieder das ESRI-Shapefile der Stadt mit den integrierten Altersquotienten pro Stadtviertel. Prinzipiell laesst sich hier noch viel mehr machen, naemlich auch mit Daten aus Drittquellen. Dazu komme ich in einem spaeteren Post.

Zuerst holen wir uns einige Bibliotheken, setzen das Arbeitsverzeichnis und lesen das Shapefile ein. Ich binde in diesem Tutorial die Libraries immer erst dann ein, wenn eine ihrer Funktionen zum Einsatz kommt. Ich hoffe, so ist halbwegs nachzuvollziehen, was wozu gehoert. Die rgdal-Library erfordert PROJ.4, das einige Zeit lang in Debian und seinen Derivaten in einer kaputten Version daher kam; unter Ubuntu lohnt hier ein Blick auf UbuntuGIS.

library(sp)
library(rgdal)
setwd("/media/home/opendata/altersquotient")
ulm <- readOGR(".", "G_Altersquotient_311210")

Mit einem Blick in die ersten Zeilen des Datenteils unseres Shapefiles faellt auf, dass die Codierung (wieder) etwas gelitten hat. Das laesst sich durch Konvertierung und Ueberschreiben der urspruenglichen Spalten beheben.

head(ulm@data)
ulm@data$STT_Name <- iconv(ulm@data$STT_Name, "ISO_8859-1", "UTF-8")
ulm@data$STTV_Name <- iconv(ulm@data$STTV_Name, "ISO_8859-1", "UTF-8")

Im Folgenden trennen wir den Datenteil vom eigentlichen Shape und wandeln letzteres in Koordinatenreihen um. Um hinterher die einzelnen Teile wieder sinnvoll zusammenfuehren zu koennen, verpassen wir dem ulm-Objekt eine Spalte id, die einfach der Zeilennummer entspricht. Anschliessend speichern wir den Datenteil von ulm nach ulm.df

ulm@data$id <- rownames(ulm@data)
ulm.df <- as.data.frame(ulm)

Die bekannte Library maptools hilft uns nun mit seiner fortify()-Methode, aus dem Shapefile Koordinatenreihen zu machen. Leider ist das nicht so ganz GNU-konform, so dass wir erst einmal mit gpclibPermit() die Erlaubnis geben muessen, diese Teilfunktion auch zu benutzen. (Danke an diese Tutorials und Beispiele, die mich auf die Sache mit dem fortify() gebracht haben)

library(maptools)
gpclibPermit()
ulm_fort <- fortify(ulm, region="id")

Nun muessen ulm.df und ulm_fort wieder zu einem gemeinsamen Data Frame kombiniert werden. Das uebernimmt die join()-Funktion aus der plyr-Library, die wie ein SQL-Join funktioniert. Bei komplexeren Verknuepfungen sollte man das im Kopf behalten und gegebenenfalls explizit einen left oder right oder inner join machen, um nicht zu verzweifeln 😉

Wichtig ist natuerlich, dass beide Teile eine gemeinsame Spalte haben, anhand derer sie verbunden werden koennen. Mit einem head() sehen wir, dass das wunderbar der Fall ist.

library(plyr)
head(ulm.df)
head(ulm_fort)

ulm_merged <- join(ulm_fort, ulm.df, by=“id“)

Nun wird geplottet — mit dem altbekannten RColorBrewer und der ggplot2-Library. Wir fuegen hinzu:

  • Eine Aesthetikangabe: Welche Werte sollen fuer x, y und die Fuellung der Flaechen verwendet werden? Wir faerben anhand des Altersquotienten ein.
  • Die Angabe, aus den Koordinatenreihen Polygone zu zeichnen
  • Zusaetzlich weisse Umrisse um die Stadtviertel
  • Eine Farbskala fuer die Altersquotientenfuellung. colours = brewer.pal(9,“Greens“) wuerde schon reichen, dann gaebe es aber nur fuenf Abstufungen. Stattdessen definieren wir die breaks in Zehnerabstaenden, indem wir einen geeigneten Vektor uebergeben.

library(ggplot2)
library(RColorBrewer)
ggplot(ulm_merged) +
aes(long, lat, group=group, fill=ALTQU) +
geom_polygon() +
geom_path (color="white") +
scale_fill_gradientn(colours = brewer.pal(9, "Greens"), breaks = c(10, 20, 30, 40, 50, 60, 70, 80, 90)) +
coord_equal()

Et voila — eine wunderbare Karte mit schoenem, unaufdringlichem Gitternetz und abgesetztem Hintergrund 🙂

R ist jedoch nicht „nur“ dazu da, bunte Karten zu plotten, zumal ohne raeumlichen Bezug und (in diesem Fall) im Gauss-Krueger-Koordinatensystem. Die hier verwendeten Libraries koennen jedoch noch viel mehr: Komplexe Kombinationen von Daten, und vor allem auch die Ausgabe in verschiedensten Formaten. Google Fusion Tables? GeoJSON? Kein Problem.

Mehr dazu irgendwann demnaechst 😉

Karten mit R

Ich haette mich in den letzten Tagen eigentlich mit etwas voellig anderem beschaeftigen muessen. Deswegen habe ich mich in R eingearbeitet. Prokrastination, level: 14. Semester. Und weil mir das Spass gemacht hat, moechte ich meine Erfahrungen festhalten und teilen.

Am Ende dieses kleinen Durchlaufs steht eine Choroplethenkarte — also eine in Subregionen unterteilte Karte, deren Flaechenfaerbung gewisse Kennzahlen anzeigen kann. So etwas geht heutzutage auch schoen als Overlay-Karte mit OpenStreetMap, mich interessierte aber der Prozess in R, und wie man ueberhaupt einmal zu den passenden Farben kommt 😉

0: R installieren

R gibt es fuer Win, MacOS und Linux; ich habe mich mit der letzteren Variante beschaeftigt. Eigentlich ginge das wunderbar mit sudo apt-get install R, dann wird aber bei einigen Distributionen erst einmal eine aeltere Version installiert, die viele Dinge nicht kann, die wir spaeter noch brauchen. Hier ist beschrieben, wie man die passenden Paketquellen einbindet, sofern man nicht ohnehin aus dem Source heraus kompilieren moechte.

Wenn man moechte, kann man auch eine GUI verwenden. Ich habe ein wenig mit RKward gespielt, mit dem man schoen die zu erledigenden Schritte in eine Skriptdatei schreiben und Zeile fuer Zeile abarbeiten kann — das hilft, die eigenen Schritte hinterher auch festzuhalten. Prinzipiell geht aber alles auch in der Konsole → R ausfuehren und anfangen. Mit q() geht’s wieder zurueck in die Shell.

(Mittlerweile bin ich auch auf Deducer gestossen, der recht vielversprechend aussieht, dazu vielleicht spaeter noch etwas)

1: Erste Schritte

Ich werde nicht grossartig auf die R-Grundlagen eingehen — die eignet man sich am besten dann bei, wenn man irgendwo auf ein Problem stoesst. Das kurze Grundlagenhandbuch von Thomas Petzold hat mir hierfuer gut getaugt.

Die Uebersicht aller Ulmer Stadtviertel liegt auf ulmapi.de als ZIP-Archiv mit einem ESRI-Shapefile vor. Und von dort zur ersten primitiven Karte sind es eigentlich nur vier Zeilen:

library(maptools)
setwd("/media/home/opendata/stadtviertel")
ulm <- readShapePoly("Stadtviertel_Gesamt_270209.shp")
plot(ulm)

Falls die maptools- (oder eine beliebige andere) Bibliothek nicht vorhanden ist, kann sie aus R heraus mit install.packages("maptools") nachinstalliert werden. Gegebenenfalls ist hierfuer noch der gcc- und gfortran-Compiler zu installieren. Das Argument in setwd() sollte natuerlich das eigene Arbeitsverzeichnis sein, in dem die Shapes aus dem ZIP-Archiv liegen.

Mit diesen wenigen Zeilen haben wir also eine erstens zum schreien haessliche und zweitens vollkommen informationsfreie Karte mit den Umrissen der Ulmer Stadtviertel hinbekommen. Dieser Schritt bringt uns also ueberhaupt nichts nuetzliches, aber zumindest einmal das gute Gefuehl, dass da etwas geklappt hat 🙂

2: Farbe!

Interessanter wird es, wenn tatsaechlich auch Farbe ins Spiel kommt. Die Shapefiles bringen (in der Regel) eine DBase-dbf-Datei mit, in denen mit den Umrissen verknuepfte Daten mitgebracht werden. Die liegen nun im Datenslot der Variable ulm, die wir vorhin mit dem Shape gefuellt haben. Ein einfaches

ulm

wirft uns den kompletten Inhalt dieser Variable um die Ohren — also auch alle Koordinaten. Mit

names(ulm)

sehen wir die Namen der Datenvektoren, und mit

ulm@data

bekommen wir diese als „Tabelle“ angezeigt. Einzelne Spalten (bzw. Vektoren) koennen wir mit einem angehaengten $Spaltenname selektieren; so gibt

ulm@data$ST_NAME (oder einfach ulm$ST_NAME)

alle Eintraege im Datenvektor der Stadtteilnamen aus. Nach diesen Eintraegen koennen wir nun arbeiten — zum Beispiel, indem wir einfach einmal die Karte nach Stadtteilen faerben. Das geht in zwei Zeilen:

col <- rainbow(length(levels(ulm@data$ST_NAME)))
spplot(ulm, "ST_NAME", col.regions=col, main="Stadtviertel Ulms", sub="cc-by-sa, Datensatz der Stadt Ulm", lwd=.8, col="white")

rainbow() erzeugt hier einfach einen Farbverlauf — oder besser gesagt, einen Zeichenvektor mit so vielen (length()) RGB-Farben, wie es eindeutige Werte (levels()) im Vektor ST_NAME gibt. spplot() plottet hier ulm, faerbt die durch ST_NAME bezeichneten Regionen mit der vorhin erstellten Palette ein, setzt Titel und Untertitel und die Grenzen zwischen den einzelnen Flaechen auf weisse Linien.

(Zusammengebastelt anhand dieser Tutorials)

Auch diese Darstellung ist in der Darstellung eher unspannend — dargestellt werden so nur Raumordnungsdaten, keine statistischen Werte. Dazu kommen wir nun.

3: Malen nach Zahlen

Eigentlich sollen Zahlenwerte dargestellt werden, wie im R Choropleth Challenge demonstriert. Hierzu brauchen wir aber erst einmal die darzustellenden Daten samt Schluesseln, um sie den passenden Regionen zuordnen zu koennen.

Die Tutorials und Anleitungen, die ich durchgeackert habe, gehen von zwei Praemissen aus — entweder werden RData-Objekte aus der GADM-Datenbank verwendet, wie in diesem Tutorial ueber die Einwohnerdichte indischer Bundesstaaten. Oder es geht um ESRI-Shapefiles, die in ihrem dbf-Anhaengsel deutlich mehr Nutzdaten mitbringen als unsere Stadtteilekarte. Schoene Beispiele hierfuer finden sich auf dieser Uebersichtsseite sowie in diesem weiterreichenden Beispiel, beides von der University of Oregon, an der gerade (Fruehling 2012) offenbar eine passende Vorlesung laeuft.

Die Zuordnung von Shape und geschluesselten CSV-Dateien greife ich (vielleicht ^^) in einem Folgeartikel noch einmal auf. Stattdessen gibt es nun was ganz exklusives, eine UlmAPI-Sneak-Preview quasi: Einen nagelneuen Datensatz der Stadt, bevor er auf ulmapi.de landet. Wow!

Es handelt sich um *trommelwirbel* ein Shapefile, das neben den Schluesseln und Namen auch den Altersquotienten fuer das jeweilige Viertel beinhaltet. Wow! 😀

Der Altersquotient gibt das Verhaeltnis von Personen ab 65 zu Personen von 14–64 an und kann somit als Indikator fuer die Altersstruktur der Bevoelkerung dienen. Und das laesst sich natuerlich gut farbig darstellen 🙂

Wir fangen wieder mit unseren bisherigen Bibliotheken und dem Einlesen an:

library(sp)
library(maptools)
setwd("/media/home/opendata/altersquotient")
ulm <- readShapePoly("G_Altersquotient_311210.shp")
# Inhalte kontrollieren
names(ulm)
ulm@data

Na hoppla, da ist was schiefgegangen, oder? Die Chancen sind gross, dass hier nun etwas wie „ M\xe4hringen“ steht. Das ist aber nichts, was man nicht schnell beheben koennte — wir konvertieren die betroffenen Vektoren korrekt nach UTF-8 und ersetzen die Originale damit:

ulm@data$STT_Name <- iconv(ulm@data$STT_Name, "ISO_8859-1", "UTF-8")
ulm@data$STTV_Name <- iconv(ulm@data$STTV_Name, "ISO_8859-1", "UTF-8")

Um nun die Daten aus der Spalte ALTQU in Farben umzusetzen, verwenden wir das Paket RColorBrewer, das einige fuer Diagrammvisualisierungen praktische Farbverlaeufe mitbringt, und die Bibliothek classInt, um die Werte unseres Altersquotientenvektors zu klassieren:

library(RColorBrewer)
library(classInt)

plotvar <- ulm@data$ALTQU
nclr <- 7
plotclr <- brewer.pal(nclr,"Greens")
class <- classIntervals(plotvar, nclr, style="equal")
# class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)

in diesem Beispiel sind die Klassenbreiten konstant; classIntervals() laesst jedoch natuerlich auch noch andere Optionen zu.

Diese Vorbereitung reicht schon, um die oben abgebildete Karte anzuzeigen:

spplot(ulm, "ALTQU", col.regions=plotclr, at=round(class$brks, digits=2), main = "Altersquotient nach Stadtvierteln in Ulm", sub="Stand 31.12.2010 nach Daten der Stadt Ulm, cc-by-sa", lwd=.8, col="white")

Weiteres Tuning waere beispielsweise durch die Anzeige der Stadtviertelnamen moeglich, was aber angesichts der teilweise dichten Haeufung nur zur groben Orientierung taugt:

spplot(ulm, "ALTQU", col.regions=plotclr, at=round(class$brks, digits=2), main = "Altersquotient nach Stadtvierteln in Ulm", sub="Stand 31.12.2010 nach Daten der Stadt Ulm, cc-by-sa", lwd=.8, col="white", sp.layout = list("sp.text", coordinates(ulm), as.character(ulm$STTV_Name), cex=0.6))

So weit so gut…

Dieser Einstieg kratzt leider wirklich nur an der Oberflaeche — und auch die Darstellung mit spplot laesst ein wenig zu wuenschen uebrig. Ebenfalls nicht behandelt habe ich das „Problem“, Shapefiles mit ganz anderen Daten zu verknuepfen. Die verlinkten Tutorials helfen hier teilweise schon weiter — ich hoffe aber, die naechsten Tage zwischen Barcamporga und Arbeiten noch ein wenig Zeit fuer ein weiteres Posting ueber ggplot2 und weitere Beispiele hinzubekommen 🙂

Viele weitere Datensaetze und vor allem auch Visualisierungsjunkies findest du ausserdem auf dem OpenCityCamp am 12. und 13. Mai 2012 an der Uni Ulm 😉

(edit: peinliche Typos gefixt)