HDR-UK Summer School: Introduction to Webmapping in R

Interactive Child Obesity Map

Webmaps are particularly useful when we want to publish detailed small-area or point data across wider regions, where zooming in and out is essential in order to explore patterns.

In this practical we will use the library leaflet to make a simple interactive map. In this scenario we will make a choropleth map, i.e. a thematic map in which small areas, depicted as polygons, are shaded in proportion to the value being displayed on the map.

We can start by loading the necessary libraries:

library(readODS)
library(sf)
library(leaflet)
library(RColorBrewer)

Get Child Obesity data

For this example we will use the Child obesity and excess weight: small area level data dataset, published by Public Health England. Unfortunately it comes in open office format (.ods), which can be a little tricky to load into R (and slow). We will explore the dataset relating to the percentage of children characterised as obese at reception.

# Read data
obs_msoa <- read_ods("NCMP_data_MSOA_update_2019.ods", sheet = 3, skip = 2)

## Parsed with column specification:
## cols(
##   .default = col_character()
## )

## See spec(...) for full column specifications.

# We specify periods by the last year
colnames(obs_msoa)[5:44] <- paste0("Y", rep(2011:2018, each = 5), "_",
                                   colnames(obs_msoa)[5:44])
# Sample
head(obs_msoa, 4)

##   MSOA code MSOA name   LA code LA name Y2011_Numerator Y2011_Denominator            Y2011_% Y2011_LCI          Y2011_UCI Y2012_Numerator Y2012_Denominator Y2012_%          Y2012_LCI Y2012_UCI Y2013_Numerator Y2013_Denominator            Y2013_% Y2013_LCI Y2013_UCI Y2014_Numerator Y2014_Denominator Y2014_% Y2014_LCI Y2014_UCI Y2015_Numerator Y2015_Denominator            Y2015_% Y2015_LCI Y2015_UCI Y2016_Numerator Y2016_Denominator            Y2016_% Y2016_LCI          Y2016_UCI Y2017_Numerator Y2017_Denominator            Y2017_%          Y2017_LCI          Y2017_UCI Y2018_Numerator Y2018_Denominator            Y2018_%          Y2018_LCI          Y2018_UCI
## 1 E02006534  Adur 001 E07000223    Adur              14               225                6.2       3.7 10.199999999999999              18               230     7.8                  5        12              22               238 9.1999999999999993       6.2      13.6              24               206    11.7         8      16.7              22               235                9.4       6.3      13.8              18               244                7.4       4.7               11.4              12               273 4.4000000000000004                2.5                7.5              16               285                5.6                3.5                8.9
## 2 E02006535  Adur 002 E07000223    Adur              10               204 4.9000000000000004       2.7 8.8000000000000007              12               188     6.4                3.7      10.8              17               205 8.3000000000000007       5.2      12.9              16               186     8.6       5.4      13.5              15               215                  7       4.3      11.2              12               222                5.4       3.1 9.1999999999999993              14               263                5.3                3.2 8.6999999999999993              18               283                6.4 4.0999999999999996 9.8000000000000007
## 3 E02006536  Adur 003 E07000223    Adur               9               198                4.5       2.4                8.4              14               174       8 4.9000000000000004      13.1              17               198                8.6       5.4      13.3              16               167     9.6         6        15              14               194                7.2       4.3      11.7              14               211                6.6         4               10.8              19               249                7.6 4.9000000000000004               11.6              22               240 9.1999999999999993                6.1               13.5
## 4 E02006537  Adur 004 E07000223    Adur              16               261                6.1       3.8 9.6999999999999993              21               276     7.6                  5      11.4              31               314                9.9         7      13.7              35               314    11.1       8.1      15.1              32               329 9.6999999999999993         7      13.4              34               346 9.8000000000000007       7.1               13.4              30               360 8.3000000000000007                5.9               11.6              25               364                6.9                4.7                9.9

# In this example we will just keep the percentages
obs_msoa <- obs_msoa[, c(1, 4, seq(7, 42, by = 5))]

# Some cleaning (Warnings are ok)
colnames(obs_msoa) <- gsub("%", "PC", colnames(obs_msoa))
colnames(obs_msoa) <- gsub(" ", "_", colnames(obs_msoa))
obs_msoa[, 3:10] <- as.data.frame(lapply(obs_msoa[, 3:10], as.numeric))

Get MSOA boundaries

In order to visualize the data into polygons we will the Middle Super Output Area Boundaries. These are part of the Census 2011 Geography. Boundary data can be found in Census Support Easy Download: English Boundary datasets section of the UK Data Service website, although a much detailed repository of spatial datasets and lookup tables can be found at the ONS Open Geography Portal.

Note that spatial data come in various formats. The most common one is arguably the ESRI Shapefile, which is a collection of files i.e. one for the spatial data (.shp), one for the dataset underneath (.dbf), one for the spatial projection (.prj), etc. However, most online repositories denote them with .SHP and offer them within a single .zip file.

There are a number of libraries that can be used to load vector data (shapefiles) into R. We will use the library sf a.k.a simple features here for simplicity and convenience.

# Read data
msoa_poly <- st_read("england_msoa_2011_sgen_clipped.shp")

## Reading layer `england_msoa_2011_sgen_clipped' from data source `C:\Users\alale\Desktop\HDRUK2019-Webmapping_in_R\england_msoa_2011_sgen_clipped.shp' using driver `ESRI Shapefile'
## Simple feature collection with 6791 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 82678 ymin: 5343 xmax: 655604.7 ymax: 657534.1
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs

# Fix projection
msoa_poly <- st_transform(msoa_poly, 4326)

# Make a simple plot of just the features
# Note that plot(msoa_poly) plots all features with all variables
plot(st_geometry(msoa_poly))

Join obesity data to spatial data

This step is perhaps the most common in analysing spatial data. We will use the ONS MSOA codes to join the obesity data to the polygon features. It is relatively straightforward:

# Merge data
obs_poly <- merge(msoa_poly, obs_msoa,  by.x = "code", by.y = "MSOA_code", all.x = T)

# Sf has its own merge method
str(obs_poly)

## Classes 'sf' and 'data.frame':   6791 obs. of  12 variables:
##  $ code    : Factor w/ 6791 levels "E02000001","E02000002",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ name    : Factor w/ 6790 levels "Adur 001","Adur 002",..: 1385 121 122 123 124 125 126 127 128 129 ...
##  $ LA_name : chr  "City of London" "Barking and Dagenham" "Barking and Dagenham" "Barking and Dagenham" ...
##  $ Y2011_PC: num  NA 13.6 10.7 16.6 13.9 16.4 14.4 15.6 12.3 11.9 ...
##  $ Y2012_PC: num  NA 11.3 12.4 16.2 12.2 15.9 12.9 14.2 12.7 12.9 ...
##  $ Y2013_PC: num  NA 12.7 11.1 9.5 11.2 16 12.7 13.1 9.4 12.5 ...
##  $ Y2014_PC: num  NA 13.4 12.2 14.4 11.2 16.1 13.9 13.1 10 11.5 ...
##  $ Y2015_PC: num  9.4 15.7 12.4 13.4 12.2 16.4 15.6 13.1 11.5 11.8 ...
##  $ Y2016_PC: num  13.3 14.5 13 13.8 12.4 16.5 16 12.9 13.6 12.3 ...
##  $ Y2017_PC: num  9.2 14.2 13.6 9 12.5 13.6 13.5 13.3 14.3 13.7 ...
##  $ Y2018_PC: num  7.1 13.6 13.7 10.8 11.2 12.5 12.1 14.1 13.3 13 ...
##  $ geometry:sfc_MULTIPOLYGON of length 6791; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:26, 1:2] -0.0928 -0.0881 -0.0852 -0.0785 -0.0774 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "code" "name" "LA_name" "Y2011_PC" ...

# Plot the latest values (2015/16 to 2017/18)
plot(obs_poly["Y2018_PC"], border = NA)

Webmap Basics

Now, we can’t see a lot of details in that scale, as there are over 6,000 MSOAs in England. We will try to make an interactive map, or webmap, where we can zoom in/out, get pop-up information by clicking at features, dynamic labels, etc.

We can create a map “widget’ object using the leaflet() function. There is a ton of features / functionality we can add in a webmap, like adding a backdrop (background map) automatically. By default this is an OpenStreetMap, but there are also other visualisations available.

An issue to look out for here is projections. Projections are kind of complicated and out of the scope of this practical, but you should know that the Leaflet package expects all point, line, and shape data to be specified in latitude and longitude using WGS84 (a.k.a. EPSG:4326).

There is a function to return the projection, within the sf package we can call st_crs. Our MSOA data have a projection of:

st_crs(obs_poly)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Which is fine, but remember to always test!. Note that many UK spatial datasets come with the OSGB 1936 / British National Grid (BNG) projection, or EPSG:27700. A common issue is that you will not get an error about this, it will be just an empty map; your actual i.e. GP data might be somewhere in the Pacific Ocean.

Now we can make a basic webmap - without any data - as a start:

obs_map <- leaflet(obs_poly) %>% addTiles() 
obs_map %>% fitBounds(-3.195992, 53.338868, -2.762375, 53.462912) # zoom to Liverpool

There is a variety of graphics that can be passed to the map. For instance, we can pass a tile layer with graphics using addTiles() or addProviderTiles(). Simply run addProviderTiles("Stamen.Toner") to change the backdrop. You can find a complete list of available backdrops by visiting http://leaflet-extras.github.io/leaflet-providers/preview/.

Child Obesity Map

As seen in the example above, adding elements to the map can be done effectively with the the pipe operator %>%, which “pipes” objects down multiple functions. Otherwise, you can just add elements to the basic map the conventional way, i.e. by overwriting the original map object (both methods are demonstrated here).

Now we can visualise the polygon data using the fillColor argument, specifying a pre-defined palette, Yellow-Orange-Red, and the value to map, in this case Y2018_PC. colorQuantile maps the value using quantiles, but we can use other methods, such as jenks or custom breaks. We can access more palettes with the RColorBrewer library.

obs_map %>% addPolygons(color = "black", 
                        weight = 1, 
                        smoothFactor = 0.5,
                        opacity = 1.0, 
                        fillOpacity = 0.5, 
                        fillColor = ~colorQuantile("YlOrRd", Y2018_PC)(Y2018_PC))

This is a basic map but we can make a few adjustments. We can remove the black borders of the areas since they are clogging up our image - we can specify a grey colour with just a lot of tranpsarency with opacity=0.1, and add a legend with addLegend(). We will not use the pipe this time, we will just add each element to the map object obs_map.

# Custom breaks
col_breaks <- c(0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 20.8)
# Make a palette
col_pal <- colorBin("YlOrRd", domain = obs_poly$Y2018_PC, bins = col_breaks)

# Make initial map
obs_map <- addPolygons(map = obs_map,
                       color = "grey",
                       weight = 1, 
                       smoothFactor = 0.5,
                       opacity = 0.1, 
                       fillOpacity = 0.75,
                       fillColor = ~col_pal(Y2018_PC),
                       highlightOptions = highlightOptions(color = "white", 
                                                            opacity = 1,
                                                            weight = 2,
                                                            bringToFront = TRUE),
                       popup = ~paste0("MSOA: ", name, "(", Y2018_PC,"%)"))

# Now add the legend
obs_map <- addLegend(obs_map, 
                      pal = col_pal, 
                      values = ~Y2018_PC, 
                      opacity = 0.9, 
                      title = "Child Obesity at Reception (%)",
                      labFormat = labelFormat(suffix="%"), 
                      # labFormat: you can add transform = function(x) x*100 if x={0,1}
                      na.label = "No Data",
                      position = "bottomright")

# Plot the map
obs_map

Note that we can highlight areas in the map with highlightOptions and include a pop-up window giving more information on the area with popup. We can also add HTML integration in the map for some increased pop-up functionality (more details here).There are many more options to explore, the best way to do this is to look at the leaflet documentation.

It is also worth noting that there is a shiny integration with leaflet. One can add scroll bars, drop-down lists, buttons etc. within the webmap widget. For those feeling adventurous, take a look here.

Optional Tasks / Excercises:

• Make a map showing the percentage of change between 2011 and 2018.
• Add the Lower and Upper CI in the pop-up.
• Try to apply a palette that is neutral (white) around the mean value.

If you want to map points you can also check another practical, Mapping the refugee crisis in Greece using leaflet, which is available here, and might be of interest.