Introduction to Webmapping using Leaflet: Child Obesity Map
Alexandros Alexiou
18/10/2020
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(sf)
library(readODS)
library(leaflet)
library(classInt)
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)
# 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\OneDrive - The University of Liverpool\Other\practical\Practical_WebMapping\england_msoa_2011_sgen_clipped.shp' using driver `ESRI Shapefile'
## Simple feature collection with 6791 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 82678 ymin: 5343 xmax: 655604.7 ymax: 657534.1
## Projected CRS: OSGB 1936 / British National Grid# 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 : chr "E02000001" "E02000002" "E02000003" "E02000004" ...
## $ name : chr "City of London 001" "Barking and Dagenham 001" "Barking and Dagenham 002" "Barking and Dagenham 003" ...
## $ 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 [1:3] "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 [1:11] "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; so we will make an interactive map, i.e. where we can zoom in and out, get more data by clicking at features, 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, the most basic one being including a backdrop (background map) automatically (by default it’s OpenStreetMap but there are other sources 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:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]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 LiverpoolThere 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 %>% which pipe 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 use 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().
# 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 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 = ~paste("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_mapNote 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 can be a shiny integration within the leaflet webmaps, and we could add scroll bars, drop-down lists, buttons etc. 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.