© Mariusz Rafało

Spatial data visualization

There are several ways to visualize spatial data in R, using common map engines like Google Maps or OpenStreetMap. Google Maps offers the functionality with the use of ggmap package. However, two issues emerged, that made the use of Google Maps difficult. These are:

  1. Google Maps API key is now required (it wasn’t, some time ago). Still you can obtain the key for free (for limited requests) but you need to reguster to Google Cloud Platform, create application, etc.

  2. There are some compatibility isses with ggmap package. There seem to be two versions aviabiable: CRAN version and github version (known as dkahle/ggmap).

This documents intends to use OpenStreetMap to make spatial analyses as an alternative approach.

Solution

Threre are multiple documents and tutorials showing how to use OpenStreetMap package with common graphics package. However, ggplot2 offers much more functionality. In this document we focus on using ggplot to plot on OpenStreetMap map fragment.

Libraries needed

Prepare environment: clear workspace

rm(list=ls())

Key package needed is OpenStreetMap. CRAN specification can be accessed here.

Package ggmap is needed only to use crime dataset.

library(OpenStreetMap)
library(ggplot2)
library(ggmap)

Load data

This example bases on crime dataset. However, thefts subdataset has been created as crime contains many records, which can be difficult to visualize.

data("crime")
thefts <- subset(crime, offense == "auto theft")

Setup map area

Bottom left and top right points coorinates are specified.

p1.lat <- 29.7
p1.lon <- -95.5
p2.lat <- 30 
p2.lon <- -95.1

Download map

OSM Comes with variety of types, e.g.: “osm”, “stamen-toner”, “stamen-terrain”. Notice that in order to plot Large OpenStreetMap map, we cannot use standard ggplot functions as it expects dataframe. In order to plot downloaded map, we use autoplot function (it also comes from ggplot2 package).

stamen-terrain
map <- openmap(c(p1.lat,p1.lon),c(p2.lat,p2.lon), zoom = 11,
               type = "stamen-terrain",
               mergeTiles = TRUE)

autoplot(map)

stamen-toner
map <- openmap(c(p1.lat,p1.lon),c(p2.lat,p2.lon), zoom = NULL,
               type = "stamen-toner",
               mergeTiles = TRUE)

autoplot(map)

osm

Create custom, rawm map, based on coordinates in theftd dataframe

map <- openmap(c(p1.lat,p1.lon),c(p2.lat,p2.lon),zoom=11,'osm')
autoplot(map)

Adjust coordinates

Both, Google Maps and OSM use Mercator projection. We nedd to transform Mercator coordinates to WGS84 coordinates. Notice, that both scales changed.

map.latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84")

autoplot(map.latlon)

Sample plots using ggplot package

autoplot(map.latlon) + xlab("lon") + ylab("lat") +
  geom_point(data=thefts, aes(x=lon, y=lat), size=3) 

h <- autoplot(map.latlon) + xlab("lon") + ylab("lat") + ylim(c(p1.lat, p2.lat)) + xlim(c(p1.lon,p2.lon))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
h + 
  geom_point(data=thefts, aes(x=lon, y=lat,color=hour), size=2) + scale_color_gradient(low = "red", high = "yellow") 
## Warning: Removed 4113 rows containing missing values (geom_point).

h +
  stat_density2d(
    aes(x = lon, y = lat, fill = number, alpha = 0.2), 
    data = thefts,
    geom = "polygon") + 
  theme(legend.position="none")
## Warning: Removed 4113 rows containing non-finite values (stat_density2d).

h +
  stat_density2d(
    aes(x = lon, y = lat, fill = ..level..), alpha=0.2,
    size = 2, bins = 6, data = thefts,
    geom = "polygon") + scale_fill_continuous(low="green", high="red")
## Warning: Removed 4113 rows containing non-finite values (stat_density2d).