Kortbaserede punkt- og tætheds-diagrammer med ggmap i R

Jeg har allerede givet eksempler på, hvordan man opretter varmekort (dvs. kortbaserede tæthedsdiagrammer) ved hjælp af deckgl og Leaflet i R. I dette indlæg vil jeg give et eksempel på, hvordan man visualiserer rumlige attributter for et datasæt ved hjælp af ggmap-pakken i R.

Jeg starter begynder med at indlæse de pakker, som jeg vil bruge til min analyse: ggmap, ggplot2, dplyr og gridExtra

#devtools::install_github("dkahle/ggmap", ref = "tidyup") 
# da ggmap i øjeblikket ikke er på CRAN
library(ggmap)
library(ggplot2)
library(dplyr)
library(gridExtra)

Jeg vil gerne bruge standarddatabasen, der er tilgængelig i R, til demonstrationsformål i dette kodeeksempel. Derfor giver jeg et glimt af datasættet ved at vise dets øverste poster i datarammen “crime”:

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

Datarammen “crime” indeholder lokationer for kriminelle handlinger. Du vil bemærke: Datasættet indeholder allerede længde- og breddegradskoordinater for alle dataindgange. Dette er de rumlige egenskaber i vores datasæt.

Som næste skridt i analysen giver jeg nu et eksempel på, hvordan basemap-fliser kan “trækkes” fra ggmap-pakken. I nedenstående kodestykke konstruerer jeg basemap-fliserne til USA.

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

“get_stamenmap”-funktionen er fra ggmap-pakken.

Vi er nu klar til at oprette et første plot baseret på vores datasets rumlige egenskaber. Nedenfor viser jeg fordelingen af mordscener, baseret på data og koordinaterne i datasættet “crime”. Funktionen “qmplot” er fra ggmap-pakken.

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

Dernæst tegner jeg et varmekort (dvs. et tæthedsdiagram). Til dette bliver jeg nødt til at specificere “geom” -parameteren i “qmplot” -funktionen til “polygon”. Jeg har også brug for funktionen “stat_density_2d” og “scale_fill_gradient2”. Densitetsestimationen er baseret på 2D-kernetæthedsestimering. Det beregnes af funktionen “stat_density_2d”.

# opret andre typer plot med ggmap-pakken
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)

I dette eksempel er visualiseringen ikke perfekt endnu og kan forbedres yderligere. Måder at gøre det på vil f.eks. vøre ved at justere beregningen af densitetsestimering.

Sørg for at tjekke mine andre indlæg om geodatavisualisering i R.