10 minutes to mapping

Mapping in R can be quite complicated and there are multiple packages to help you do it. This short tutorial based on x, y, z uses the tmap package to create an interactive map of self-harm rates by ward.

There are 5 steps.

1. Get the data you want to plot.
2. Find a shape file.
3. Do some data processing and linkage.
4. Plot the map
5. Make it interactive.

1 Step 1: Setup

For this exercise we will use a number of R packages:

• rgdal to load Shape files
• readr to load data
• dplyr for data wrangling
• janitor for cleaning up data lables

knitr::opts_chunk$set(echo = FALSE, include = FALSE)

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

suppressMessages(library(readr))
library(rgdal)

rgdal: version: 1.2-5, (SVN revision 648)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
 Path to GDAL shared files: /Users/julianflowers/Library/R/3.3/library/rgdal/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
 Path to PROJ.4 shared files: /Users/julianflowers/Library/R/3.3/library/rgdal/proj
 Linking to sp version: 1.2-4 

library(janitor)
library(tmap)
options(digits = 3)

2 Step 2: Get the data

The latest Mental Health JSNA tool has some indicators at ward level. Currently Fingertips doesn’t offer a faclity to map data at this level so this walk through shows how to create these maps in R.

We’ll download the data, convert it to .csv format and import into R.

testward <- read_csv("~/Downloads/testward.csv", col_types = c("ciccccdcccccccc")) %>%
  janitor::clean_names() %>%
  select(x_indicator, area_code, area_name, time_period, value)

15036 parsing failures.
  row         col               expected actual
30055 Time Period no trailing characters    /14
30056 Time Period no trailing characters    /14
30057 Time Period no trailing characters    /14
30058 Time Period no trailing characters    /14
30059 Time Period no trailing characters    /14
..... ........... ...................... ......
See problems(...) for more details.

testward <- testward %>%
  filter(x_indicator == "Hospital admissions for self-harm: standardised emergency admission ratio (all ages)") 

3 Step 3: Get a relevant shape file

Shape files for ward boundaries can be obtained from PHE’s ESRI server or from the ONS Open geography portal. In this example we’ll use the super generalised ward boundaries for 2015 from the Open Geography portal. We’ll use the rgdal package to read in the shape file as a Spatial Data Frame.

shp <- tempfile() ## create a temp file
download.file('http://geoportal.statistics.gov.uk/datasets/5fb8813978cc4e4892da4b57bcf4491f_3.zip', shp) ## download the file

trying URL 'http://geoportal.statistics.gov.uk/datasets/5fb8813978cc4e4892da4b57bcf4491f_3.zip'
downloaded 2.5 MB

unzip(shp, exdir = getwd()) ## unzip
s <- readOGR(dsn = "Wards_December_2015_Super_Generalised_Clipped_Boundaries_in_Great_Britain.shp", layer = "Wards_December_2015_Super_Generalised_Clipped_Boundaries_in_Great_Britain", verbose = FALSE) ## read in th shape file

4 Data linkage and pre-processing

seng <- subset(s, substr(wd15cd, 1, 1) == "E") ## restrict to England wards
seng@data <- left_join(seng@data, testward, by=c("wd15cd" = "area_code")) ## Join indicator data

joining character vector and factor, coercing into character vector

seng@data$country <- "Eng" ## restrict to England
library(viridis) ## colour palette
palette <- rev(plasma(15)) ## colour scheme
indicator <- stringr::str_wrap(seng@data$x_indicator, width = 30) ## labels
style <- "kmeans" ## categorisation
n <-  10 ## categories
credits <- "Contains ordnance survey data (c) \nCrown copyright and database right 2016"

5 Draw map

t <- tm_shape(seng) +
  tm_fill("value", style = style , n = n,
          palette = palette, title = "Self-harm rate SAR:\n100 = national average", 
          popup.vars = c("x_indicator", "value", "area_name")) +
  tm_layout(paste(indicator, ":", style),
            title.size = 1,
            legend.outside = TRUE, 
            frame = FALSE) +
  tm_scale_bar() +
  tm_compass()+
   tm_credits( credits, size = 0.7, align = "right") 
tmap_mode("plot")
t

6 Make it interactive

tl <- tmap_leaflet(t)

Credits not supported in view mode.Compass not supported in view mode.Scale bar not yet supported in view mode, it will be in the next version.

tl

That’s it. There are a lot more tweaks you can make of course.