I really liked Ramiro Gómez’ rendering of the density of pub locations in the UK and Ireland at http://ramiro.org/notebook/mapping-pubs/, and just thought I would replicate his code from Matplotlib to R – also to point that very few lines of code are needed in R to produce a decent print quality map. Moreover the good folks at RStudio are also making considerable progress developing bindings between R and Leaflet Javascript library, so here is an example of how to turn Ramiro’s map into a fully interactive version.

First, we need to load a couple of workhorse R libraries.

library(data.table)
library(raster)
## Loading required package: sp
library(tmap)
library(leaflet)

I always use data.table for data wrangling operation (it’s incredibly fast with very concise syntax). data.table will load osm-x-tractor 800MB CSV file in under a second! raster is not strictly needed here (but useful to get access to GADM admnistrative boundary maps). raster will also load all necessary libraries to create and manipuate spatial features.

tmap is Martijn Tennekes’ excellent convenience methods to create quick and beautiful print maps. I highly recommend his introduction tmap in a Nutshell. Yes, you can use R base plot() and spplot() methods to print nice choropleth maps as well, but tmap goes a step further in simplicity.

We start by loading a complete CSV file in our R workspace. Note that osm-x-tractor requires that you download the file by hand (there’s no permanent URL). No need to use csvgrep though. I saved the archive in ./temp/POIEurope_2015.06.20.7z and then had to uncompress it.

# Uncompress 7z archive
system("7za e ./temp/POIEurope_2015.06.20.7z")  

# Read it in, note fread() would be much faster, but these CSV files have skipped 
# field delimiters, so we need to use fill=T
# poi <- fread("./temp/POIEurope.csv")
poi <- read.csv("./temp/POIEurope.csv", fill=T)
poi <- data.table(poi)
# Explore
kable(data.frame(colName=names(poi), colClass=sapply(poi, class)), row.names=F)
colName colClass
osmid numeric
name factor
amenity factor
emergency factor
geological factor
historic factor
leisure factor
man_made factor
office factor
shop factor
sport factor
tourism factor
craft factor
Longitude numeric
Latitude numeric 
# Sumarize the number of locations across amenity classes (show only top 20 rows)
kable(poi[, .N, keyby=amenity][order(-N)][1:20], row.names=F)
amenity N
2824010
bench 442592
restaurant 234048
parking 205450
post_box 156604
recycling 129263
waste_basket 104173
cafe 103737
place_of_worship 101117
fuel 95373
bank 86800
pharmacy 84147
bicycle_parking 83492
pub 75770
fast_food 74551
drinking_water 71761
telephone 70070
school 66235
hunting_stand 65506
atm 55929 
# Keep only pubs, and only the columns we need
poi <- poi[amenity=="pub", list(name, amenity, Longitude, Latitude)]

# Also remove pubs with no coordinate
poi <- poi[!is.na(Longitude) & !is.na(Latitude)]

There are 75770 pubs to be mapped across Europe. Next we load province boundaries for the UK and Ireland, and convert the list of pubs to a list of spatial points.

# Use raster::getData() to conveniently return admin boundaries
gbr <- getData("GADM", country="GBR", level=1)
irl <- getData("GADM", country="IRL", level=1)

# We can merge these 2 SpatialPolygonsDataframe objects into 1 for convenience
# Good practice is to ensure they have similar projection before merging
# and then generate unique polygon IDs with spChFIDs()
gbr <- spTransform(gbr, CRS("+init=epsg:4326"))
irl <- spTransform(irl, proj4string(gbr))

gbr <- spChFIDs(gbr, paste0("gbr", row.names(gbr)))
irl <- spChFIDs(irl, paste0("irl", row.names(irl)))

# Merge features
adm1 <- rbind(gbr, irl)

# Let's map it
tm_shape(adm1) +
  tm_grid(n.x=5, n.y=5, col="grey90", labels.col="grey50", on.top=F) +
  tm_polygons(col="NAME_1", border.col="white", palette="Paired",
    title="Provinces", labels=sort(unique(adm1$NAME_1)), max.categories=30) +
    tm_layout(draw.frame=F,
      legend.position=c(-0.1,.05), legend.bg.color="white")

We’re ready to convert the list of pubs to a SpatialPointsDataframe, and then limit its extent to Ireland and UK.

poi.sp <- SpatialPointsDataFrame(poi[, list(Longitude, Latitude)], 
  data.frame(poi), proj4string=CRS("+init=epsg:4326"))

# Simply crop the points to the map
poi.sp <- crop(poi.sp, extent(adm1))

Now using tmap package we can generate a bubble map for print.

# Let's overlay the pub locations
tm_shape(adm1, bbox=bbox(adm1)-c(-1,1)) +
  tm_polygons(col="#f8f5de", border.col=NA) +
  tm_shape(poi.sp) + tm_bubbles(col="#e5662b", alpha=.8, size=0.004) +
  tm_credits("Author: Melanie Bacou. Data at OpenStreetMap - openstreetmap.org", size=.8) +
  tm_layout(title="United Kingdom \nand Ireland... \nPub-wise", 
    bg.color="#73B6E6", draw.frame=F)

Else we can also draw an interactive map using leaflet.

  # Create the map
  m <- leaflet(data=adm1) %>%
  setView(0, 52, 7) %>%
  # Add basemap
  addTiles(
    urlTemplate="http://{s}.tiles.mapbox.com/v3/jcheng.map-5ebohr46/{z}/{x}/{y}.png",
    attribution="Mapbox", group="Basemap") %>%
  # Add pub locations clustered
  addCircleMarkers(data=poi.sp, layerId=~name, group="Pub location",
    radius=10, fillColor="steelblue", stroke=F,
    clusterOptions=markerClusterOptions(showCoverageOnHover=T)) %>%
  # Add layer control
  addLayersControl(
    overlayGroups="Pub location",
    options=layersControlOptions(collapsed=F))  

# Show
m