Mit Hilfe von R lassen sich flächenbezogenene Daten auf Karten visualisieren. Nachfolgend wird ein mögliches Vorgehen für die Erstellung einer Chorophletenkarte kurz dargestellt. Eine Choroplethenkarte ist eine thematische Karte, bei der die Gebiete gemäß der flächenbezogenen Daten eingefärbt ist. Als Beispiel dient die Bevölkerungsveränderung der 402 deutschen kreisefreien Städte und Landkreise zwischen dem Jahr 2012 und 2013.

Das hier dargestellte Vorgehen orientiert sich an einem Blogeintrag von Kevin Johnson und verwendet maßgeblich das R-Package ggplot2.

library(ggplot2)
library(sp)
library(rgdal)
library(knitr)
library(dplyr)
library(rmarkdown)
library(scales)

Für die Visualisierung von geodaten-bezogenenen Informationen auf einer Karte müssen die Kartendaten (hier die Kreisgrenzen) den zu visualisierenden Daten (hier der Bevölkerungsveränderung) verheiratet werden.

Die Kartendaten liegen als Shapefile vor. Ein Shapefile ist keine einzelne Datei, sondern besteht aus verschiedenen Dateien, die anhand eines gemeinsamen layers miteinander verbunden sind. Die hier verwendeten Kartendaten zu den Kreisgrenzen in Deutschland stammen von dem Dienstleistungszentrum des Bundes für Geoinformation und Geodäsie.

Die Daten für die Bevölkerungsveränderung stammen von der Regionaldatenbank Deutschland. Sie sind in der csv-Datei name.data gespeichert.

name.data <- "EWZ_Kreis_2012_2013.csv"
layer.shp <- "VG250_KRS"

Im ersten Schritt muss das Shapefile eingelesen werden. Hierzu steht der readOGR-Befehl zur Verfügung. Nach dem erfolgreichen Import wird durch den Befehl names eine Übersicht über die vorhandenen Spalten des Shapefiles angezeigt. Hierdurch ist es möglich einen eindeutigen Identifikationsschlüssel in dem Shapefile zu identifizieren.

kreise.shp <- readOGR(dsn="SpatialData", layer=layer.shp)
names(kreise.shp)
##  [1] "ADE"      "GF"       "BSG"      "RS"       "AGS"      "SDV_RS"  
##  [7] "GEN"      "BEZ"      "IBZ"      "BEM"      "NBD"      "SN_L"    
## [13] "SN_R"     "SN_K"     "SN_V1"    "SN_V2"    "SN_G"     "FK_S3"   
## [19] "NUTS"     "RS_0"     "AGS_0"    "WSK"      "EWZ"      "DEBKG_ID"

Das Koordinatenreferenzsystem (CRS) des Shapefiles lässt sich mit dem Befehl proj4string anzeigen und kann nach dem Einlesen verändert werden. In diesem Beispiel wird die in Deutschland häufig verwendete Projektion EPSG:4839 verwendet.

proj4string(kreise.shp)
## [1] "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"
kreise.shp <- spTransform(kreise.shp, CRS=CRS("+init=epsg:4839"))

Die Übersicht über die Struktur des Shapefiles hat gezeigt, dass verschiedene Identifikationsschlüssel vorliegen. In diesem Beispiel wird als Identifikationsschlüssel der fünfstellige allgemeine Gemeindeschlüssel (AGS) verwendet. Mit Hilfe des AGS werden 1. die Umformung des Shapefiles in ein Dataframe mit Hilfe des Befehls fortify 2. die zu visualisierenden Daten mit den Kartendaten in dem Dataframe verheiratet

Die Umformung des Shapefiles in ein Dataframe ist notwendig, da in dem R-Package ggplot2 die zu visualisierenden Daten in einem Dataframe vorgehalten werden müssen. Der Befehl str zeigt die Struktur des transformierten Dataframes an. Die Informationen innerhalb des Dataframes stellen einzelne Koordinaten der Kreisgrenzen dar (long und lat). Diese Koordinaten werden später mit ggplot und der zugehörigen Funktion geom_polygon miteinander verbunden. Durch id=AGS ist sichergestellt, dass die Polygon-Züge die jeweiligen Kreisgrenzen darstellen. Hierfür ist auch eine Reihenfolge order angegeben, in welcher Reihenfolge die Punkte miteinander verbunden werden sollen.

kreise.df <- fortify(kreise.shp, region="AGS")
str(kreise.df)
## 'data.frame':    525831 obs. of  7 variables:
##  $ long : num  -69965 -69934 -69922 -69829 -69747 ...
##  $ lat  : num  425790 425633 425634 425641 425684 ...
##  $ order: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ hole : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ piece: Factor w/ 36 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ group: Factor w/ 657 levels "01001.1","01002.1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ id   : chr  "01001" "01001" "01001" "01001" ...

Im zweiten Schritt müssen die flächenbezogenen Daten importiert werden. Die flächenbezogenen Daten liegen in name.data vor und werden per read.csv importiert. Nach dem Import wird der fünfstellige Identifkationsschlüssel AGS als String transformiert und mit einer führenden 0 versehen, falls diese beim Import gelöscht worden ist.

kreise.data <- read.csv(file=paste("Data",name.data,sep="/"), header=TRUE, sep=";", dec=",", stringsAsFactors=TRUE)
kreise.data$AGS <- as.character(kreise.data$AGS)
kreise.data$AGS <- ifelse(nchar(kreise.data$AGS)<5, paste0("0", kreise.data$AGS), kreise.data$AGS)
str(kreise.data)
## 'data.frame':    402 obs. of  13 variables:
##  $ AGS        : chr  "01001" "01002" "01003" "01004" ...
##  $ SDV_RS     : num  1.00e+10 1.00e+10 1.00e+10 1.00e+10 1.05e+10 ...
##  $ GEN        : Factor w/ 380 levels "Ahrweiler","Aichach-Friedberg",..: 93 159 182 219 66 140 226 251 258 260 ...
##  $ IBZ        : int  40 40 40 40 42 42 42 42 42 42 ...
##  $ NBD        : Factor w/ 2 levels "ja","nein": 1 1 1 1 1 1 1 1 1 1 ...
##  $ RS_0       : num  1.00e+10 1.00e+10 1.00e+10 1.00e+10 1.05e+10 ...
##  $ AGS_0      : int  1001000 1002000 1003000 1004000 1051000 1053000 1054000 1055000 1056000 1057000 ...
##  $ Zuwanderung: int  7701 16354 12136 7518 11012 14922 13928 14925 23155 9311 ...
##  $ Fortzüge   : int  6980 14555 9951 7125 10571 13315 13509 13571 20205 8758 ...
##  $ Gestorbene : int  1057 2483 2796 951 1758 2135 1952 2712 3327 1503 ...
##  $ Geborene   : int  764 2255 1766 636 948 1475 1119 1200 2453 818 ...
##  $ EWZ_2013   : int  83971 241533 212958 77058 132665 189043 161923 197835 301223 126643 ...
##  $ EWZ_2012   : int  83462 239866 211713 76951 132965 187905 162237 197882 298826 126721 ...

Die Bevölkerungsveränderung ist noch nicht Bestandteil des Dataframes. Entsprechend wird die Veränderung der Einwohnerzahl in Prozent aus den vorhanden Informationen erzeugt. Hierzu wird das dplyr-Package verwendet.

kreise.data <- mutate(kreise.data,
                      id = AGS,
                      DifferenzEWZ = EWZ_2013 - EWZ_2012,
                      DifferenzEWZ.proz = (DifferenzEWZ/EWZ_2012)*100
)

Anschließend werden die Kartendaten mit den flächenbezogenen Daten verheiratet. Hierfür wird der Befehl left_join aus dem dplyr-Package verwendet. Nachdem dies erfolgt wird zur Sicherheit die gewünschte Sortierung des Dataframe hergestellt.

plot.data <- left_join(x=kreise.df, y=kreise.data, by="id")
plot.data <- arrange(plot.data, id, order)

Das Dataframe plot.data enthält alle Informationen, die für die Choroplethenkarte notwendig sind. Nachfolgend beispielhaft dargestellt für die ersten 5 Koordinaten der Stadt Weimar mit dem allgemeinen Gemeindeschlüssel 16055.

##      id   AGS  long  lat  order   group EWZ_2013 EWZ_2012
## 1 16055 16055 52463 4200 503867 16055.1    63315    63236
## 2 16055 16055 52463 4262 503868 16055.1    63315    63236
## 3 16055 16055 52468 4324 503869 16055.1    63315    63236
## 4 16055 16055 52614 4350 503870 16055.1    63315    63236
## 5 16055 16055 52695 4187 503871 16055.1    63315    63236
##   DifferenzEWZ.proz
## 1         0.1249288
## 2         0.1249288
## 3         0.1249288
## 4         0.1249288
## 5         0.1249288

Für die Erstellung der Choroplethenkarte wird ggplot2 verwendet. Die durch ggplot2 erstellen Grafiken sind nach dem Prinzip “layer by layer” aufgebaut. Das geometrische Objekt (geom), das die Kreisgrenzen beinhaltet, wird durch geom_polygon erzeugt. Die farbliche Kennzeichnung wird über die Einstellung fill gesteuert. In diesem Beispiel soll die farbliche Kennzeichnung in Abhängigkeit der Variable DifferenzEWZ.proz erfolgen.

g <- ggplot(data=plot.data) 
g <- g+ geom_polygon(aes(x=long,y=lat, group=group, fill=DifferenzEWZ.proz), alpha=0.9)
g <- g + scale_fill_gradient2(low = muted("green"), mid="white", high= muted("blue"), midpoint=0, na.value="red")
g <- g + coord_equal()
g <- g + guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10, title="Veränderung in %"))
g <- g + myMapTheme
g

Für eine bessere Orientierung soll als weiterer layer die Grenzen der Bundesländer der Choroplethenkarte hinzugefügt werden. Das Vorgehen ist identisch mit dem oben beschriebenen Vorgen für die Aufbereitung des Shapefiles. Da die Bundesländergrenzen des Geodatenzentrum den territorialen Grenzen Deutschlands entsprechen und nicht den Küstenregionen folgen, wird ein Shapefile vom GADM verwendet. Die transformierten Kartendaten für die Bundesländergrenzen liegen anschließend in dem Dataframe laender.df vor und kann der Choroplethenkarte als neuer layer hinzugefügt werden.

g2 <- g + geom_polygon(data=laender.df, aes(x=long,y=lat, group=group), fill=NA, colour="grey50")

Im Ergebnis ergibt sich folgende Choroplethenkarte, die die gewünschten Informationen enthält und weiter optimiert werden kann.