Kartenbasierte Punkt- und Dichtediagramme mit ggmap in R

Ich habe bereits Beispiele zur Erstellung von Wärmekarten (d.h. Dichtediagrammen) mit Deckgl und Leaflet in R bereitgestellt. In diesem Beitrag möchte ich ein Beispiel zur Visualisierung räumlicher Attribute in Datensatzen mithilfe des ggmap-Pakets in R teilen. Ich beginne mit dem Import der relevanten Pakete:

#devtools::install_github("dkahle/ggmap", ref = "tidyup") 
# da ggmap derzeit nicht auf CRAN ist
library(ggmap)
library(ggplot2)
library(dplyr)
library(gridExtra)

In diesem Codierungsbeispiel möchte ich die in R verfügbare Standarddatenbank zu Demonstrationszwecken verwenden. Daher gebe ich einen Einblick in diesen Datensatz indem ich die oversten Dateneinträge zeige.

head(crime)

##                      time     date hour premise            offense  beat
## 82729 2010-01-01 07:00:00 1/1/2010    0     18A             murder 15E30
## 82730 2010-01-01 07:00:00 1/1/2010    0     13R            robbery 13D10
## 82731 2010-01-01 07:00:00 1/1/2010    0     20R aggravated assault 16E20
## 82732 2010-01-01 07:00:00 1/1/2010    0     20R aggravated assault  2A30
## 82733 2010-01-01 07:00:00 1/1/2010    0     20A aggravated assault 14D20
## 82734 2010-01-01 07:00:00 1/1/2010    0     20R           burglary 18F60
##           block    street type suffix number   month    day
## 82729 9600-9699   marlive   ln      -      1 january friday
## 82730 4700-4799 telephone   rd      -      1 january friday
## 82731 5000-5099  wickview   ln      -      1 january friday
## 82732 1000-1099   ashland   st      -      1 january friday
## 82733 8300-8399    canyon           -      1 january friday
## 82734 9300-9399     rowan   ln      -      1 january friday
##                       location           address       lon      lat
## 82729    apartment parking lot   9650 marlive ln -95.43739 29.67790
## 82730 road / street / sidewalk 4750 telephone rd -95.29888 29.69171
## 82731        residence / house  5050 wickview ln -95.45586 29.59922
## 82732        residence / house   1050 ashland st -95.40334 29.79024
## 82733                apartment       8350 canyon -95.37791 29.67063
## 82734        residence / house     9350 rowan ln -95.54830 29.70223

Als nächstes zeige beispielhaft wie Grundkartenkacheln aus dem ggmap-Paket „gezogen“ werden können. Im folgenden Code erstelle ich die Grundkartenkacheln für die USA. Man stellt fest: Der Datensatz enthält bereits Längen- und Breitengradkoordinaten für alle Dateneinträge. Dies ist die räumliche Eigenschaft unseres Datensatzes. Jeder Dateneintrag stellt einen Tatort für eine bestimmte kriminelle Handlung dar.

# plotte eine ggmap-Grundkarte
us <- c(left = -125, bottom = 25.75, right = -67, top = 49)
map <- get_stamenmap(us, zoom = 5, maptype = "toner-lite",legend="none")
plot(map)

Die Funktion „get_stamenmap“ stammt aus dem ggmap-Paket.

Wir sind jetzt soweit, dass wir ein erstes Diagramm erstellen könnne. Dieses soll auf den räumlichen Eigenschaften unseres Datensatzes basieren. Im Folgenden zeige ich die Verteilung der Tatorte auf der Grundlage der Koordinaten, die im Datensatz „Verbrechen“ angegeben sind. Die Funktion „qmplot“ stammt aus dem ggmap-Paket.

scatterplot_murder <- qmplot(x=lon,y=lat,data=filter(crime,offense=="murder"),legend="none",color=I("darkred"))
plot(scatterplot_murder)

Als nächstes werde ich eine Wärmekarte (d.h. ein Dichtediagramm) erstellen. Dazu muss ich den „geom“ -Parameter in der „qmplot“ -Funktion auf „polygon“ setzen. Außerdem muss ich die Funktionen „stat_density_2d“ und „scale_fill_gradient2“ verwenden. Die Dichteschätzung basiert auf einer 2D-Kernel-Dichteschätzung. Sie wird mit der Funktion „stat_density_2d“ berechnet.

# Erstelle mit dem ggmap-Paket weitere plots (Dichtediagramm; also Wärmekarte)
densityplot_murder <- qmplot(x=lon, y=lat, 
                             data = filter(crime,offense=="murder"), 
                             geom = "blank",
                             maptype = "toner-background", 
                             darken = .7, 
                             legend = "topright") + stat_density_2d(aes(fill = ..level..), 
                  geom = "polygon", 
                  alpha = .5,
                  color = NA) + scale_fill_gradient2(low = "blue", 
                       mid = "green", 
                       high = "red")

plot(densityplot_murder)

In diesem Beispiel ist die Visualisierung noch nicht perfekt und könnte weiter verbessert werden. Eine Möglichkeit um dies zu tun wäre bspw. das Anpassen der Dichteschätzungsberechnung.

Schaue Dir auch meine anderen Beiträge zur Visualisierung von Geodaten in R an.