Mapping and Modelling Geographic Data in R

Session 2: Visualisation in R

Introduction

In their book R for Data Science, Wickham and Grolemund introduce the various stages of ‘data wrangling’ that go into analysis. A graphical portrayal of those stages, taken from their book, is below.

source: Wickham & Grolemund

Amongst those stages is visualisation. Visualisation is used for communication – to communicate important features of the data such as their distribution or the relationship between variables. Often we think of that communication as being to an audience – the readers of an academic paper or a policy report, for example. But equally important is communication to ourselves as analysts, as we use visual tools (as well as other methods such as descriptive statistics) to help us to explore and to understand the data.

We already have learned that there is often more than one way of doing things in R and that is true of producing graphics too. R comes with the core graphics package. Open the index for that package by running the snippet of code below and have a look at some of the plots available:

help(package = "graphics")

Now modify (add to) the following block of code to produce a box plot of the variable \(y\). See if you can get the box plot to span horizontally instead of vertically across the screen. you might find the help file, help(boxplot), useful.

y <- c(1:10, 12, 15, 20)
# Replace this comment with your own code

ggplot2

A different approach – initially harder to use but good for producing publishable graphics - is to the use the ggplot2 package:

# If package is not installed then install it
if(!"ggplot2" %in% installed.packages()) install.packages("ggplot2")
# Load / require the package
require(ggplot2)

The most comprehensive book for learning ggplot2 is (again) Hadley Wickham’s. However, a number of shorter tutorials can be found online, including this one.

Another notes that there are basically three steps to producing a plot using ggplot2, that last of which is optional. These are:

• specify the aesthetic, which is what you want to graph
• choose one of more geoms, which is how you want to graph it
• add some options, such as labelling or the theme (style), and so forth.

So, for a box plot of \(y\), \(y\) is the aesthetic, boxplot is the geom. The only complication is that ggplot2 will want \(y\) to be part of a data frame:

mydata <- data.frame(y = y)
ggplot(mydata, aes(x = "", y = y)) +
  geom_boxplot()

Actually, if you look again at the code, you will find there are two aesthetics, \(x\) and \(y\). This is because boxplots are often used to compare groups (categories) of data and \(x\) is the grouping variables. We can demonstrate this will some ‘in-built’ data that comes with ggplot2.

data(mpg)
ggplot(mpg, aes(x = class, y = hwy)) +
  geom_boxplot() +
  ylab("Miles per gallon") +
  coord_flip() +
  theme_minimal()

Note that in the example above, some options have been added: specifically, changing the axis title on one axis, flipping the coordinates (so the boxplots spread horizontally not vertically) and using a more minimal theme. A useful cheatsheet for ggplot2 is available. Use it to modify the following code chunk to produce a scatterplot with hwy as the \(Y\) variable and cty as the \(X\) variable, and add a regression line of best fit to the plot. See if you can get the regression line to be plotted without a confidence interval around it.

# Replace this comment with your own code

A useful feature of ggplot2 is its ability to facet. Say, for example, we wanted a dotplot showing the distribution of the fuel efficiency (hwy) statistics for each class of vehicle:

ggplot(data = mpg, aes(hwy)) +
  geom_dotplot() +
  facet_wrap(~ class) 

If we wanted, we could also colour each dot by the year of manufacture. Trying running the following code block first as a whole and then without the bottom line - what changes?

ggplot(data = mpg, aes(hwy)) +
  geom_dotplot(aes(fill = factor(year))) +
  facet_wrap(~ class) +
  scale_fill_manual(name = "year", values = c("white", "black"))

The function ggsave can be used to save the last plot, as in the example below. Take a look at the help menu (type ?ggsave in the Console) to learn how to change the file format, size of the graphic, resolution, and other options.

ggsave("myplot.pdf")

Spatial visualisation

From a geographical perspective, what often is of interest is geographical (spatial) variation - i.e. ‘what is where?’. This might focus on geographical clustering (similar things being found in close proximity to each other) and/or geographical heterogeneity (different things being found in different sub-regions). The principle tool for visualising these spatial patterns and variations is, of course, a map.

Introducing tmap

One way to draw maps in R is to use the tmap package where the t is short for thematic - i.e. it is a package for producing thematic maps.

# If package is not installed then install it
if(!"tmap" %in% installed.packages()) install.packages("tmap")

The following code chunk is from the vignette tmap: get started! and maps the Human Poverty Index (HPI) for countries around the World:

require(tmap)
data("World")
tm_shape(World) +
    tm_polygons("HPI")

Looking at the code chunk, the first line of code (data(“World”)) simply loads a data set that comes bundled with the package (the World data). The data set is of class,

class(World)

The “data.frame” part of it is obvious enough (it includes a table of data) but what of the sf? It stands for simple feature but you might like to think of it as spatial feature - World is a mappable object, which is why we can use it to produce the choropleth map, above.

More specifically, we need two things to produce the choropleth map - the geometry of the countries (i.e. a digital representation of their boundaries) and something about the attributes of those countries (something measured about them). The attributes contained in the World map are,

require(sf)
st_drop_geometry(World)

In the following code chunk, enter the code required to produce a map of the estimated population of each country (using the variable pop_est).

# Replace this comment with your own code

Look again at the map you have just produced of the population estimates. You will probably find it is dominated by the largest countries, which we can identify in the data:

require(tidyverse)
World %>%
  st_drop_geometry() %>%
  dplyr::select(name, pop_est) %>%
  arrange(desc(pop_est)) %>%
  slice(1:5)

Looking at the distribution of the data, they have a very strong positive skew, which means that many of the map classes have few or no countries contained within them under an equal interval classification.

require(classInt)
brks <- classIntervals(World$pop_est, n = 7, style = "equal")$brks
World %>%
  ggplot(aes(pop_est)) +
  geom_histogram() +
  geom_rug() +
  geom_vline(xintercept = brks, linetype = "dotted")

We do better using a ‘natural breaks’ classification that tries to identify break points that reflect the distribution of the data.

brks <- classIntervals(World$pop_est, n = 7, style = "jenks")$brks
World %>%
  ggplot(aes(pop_est)) +
  geom_histogram() +
  geom_rug() +
  geom_vline(xintercept = brks, linetype = "dotted")

Applying those to the map is straightforward,

tm_shape(World) +
    tm_polygons("pop_est", style = "jenks", n = 7)

Use the following code chunk to produce a map wherein the population estimates are grouped into quartiles.

# Replace this comment with your own code

Adding annotation

Once we have decided upon the basic design of the map, there are various ways to add to or modify it. For example,

data(metro)
tm_shape(World) +
    tm_polygons("pop_est", style = "jenks", n = 7, palette = "Blues",
                border.alpha = 0, title = "population estimate") +
    tm_compass(type = "radar", position = c("left", "top"), size = 4) +
    tm_scale_bar(position = c("right", "bottom")) +
    tm_text("iso_a3", size = "AREA") +
tm_shape(metro) +
  tm_symbols(col = "red", size = "pop2020", scale = .5, alpha = 1)

Or,

tmap_mode("plot")
tm_shape(World) +
  tm_polygons("HPI") +
  tm_layout(bg.color = "skyblue", inner.margins = c(0, .02, .02, .02))

The map can also be changed to an ‘interactive’ mode, as in the following example …

tmap_mode("view")
tm_basemap("Stamen.Watercolor") +
tm_shape(metro) +
  tm_symbols(col = "red", size = "pop2020", scale = 2, alpha = 1)

… and in this,

tmap_mode("view")
World %>%
  filter(grepl("America", continent)) %>%
  tm_shape(.) +
    tm_polygons("HPI")

Use the code block below to produce a map of the HPI scores for Africa to which you should add a scale bar and a north arrow. The first part of the code is completed for you.

World %>%
  filter(continent == "Africa") %>%
  # Replace this comment with your own code

Faceting in tmap

There are various ways to create facets in tmap. One is to combine the variables you wish to map using the c(…) [combine] function:

tmap_mode("plot")
World %>%
  filter(continent == "Asia") %>%
  tm_shape(.) +
  tm_polygons(c("well_being", "inequality"), style = "jenks") +
  tm_facets(ncol = 2) +
  tm_legend(position = c("left", "bottom"))

… although, mapping it this way implies that well-being rises with inequality. Is that true?

World %>%
  filter(continent == "Asia") %>%
  ggplot(aes(x = inequality, y = well_being)) +
  geom_point() +
  geom_smooth(method = "lm")

… no, it isn’t! It would therefore be better to create a new equality variable and map that with well-being:

tmap_mode("plot")
World %>%
  filter(continent == "Asia") %>%
  mutate(equality = 1 - inequality) %>%
  tm_shape(.) +
  tm_polygons(c("well_being", "equality"), style = "jenks") +
  tm_facets(ncol = 2) +
  tm_legend(position = c("left", "bottom"))

Or, in ‘interactive’ mode:

tmap_mode("view")
World %>%
  filter(continent == "Asia") %>%
  mutate(equality = 1 - inequality) %>%
  tm_shape(.) +
  tm_polygons(c("well_being", "equality"), style = "jenks") +
  tm_facets(sync = TRUE, ncol = 2) +
  tm_legend(position = c("left", "bottom"))

(try setting sync to FALSE in the code above and see what happens)

In the examples above, the geographical region is held constant whilst the attribute variable changes. This can be reversed, as in the following example, which is modified from tmap: get started!

data(NLD_muni)
tmap_mode("plot")
NLD_muni %>%
  mutate(perc_men = pop_men / population * 100) %>%
  tm_shape(.) +
  tm_polygons("perc_men", palette = "RdYlBu") +
  tm_facets(by = "province")

Maps can also be arranged side-by-side using tmap’s arrange function:

NLD_muni %>%
  tm_shape(.) +
  tm_polygons("population", convert2density = TRUE) -> tm1
NLD_muni %>%
  tm_shape(.) +
  tm_bubbles(size = "population") -> tm2
tmap_arrange(tm1, tm2)

Modify the following code block so that the HPI scores for North America are mapped alongside the HPI scores for South America. The beginning of the code is included for you.

tmap_mode("plot")
World %>%
  filter(grepl("America", continent)) %>%
  # Replace this comment with your own code 

Exporting the plots

Plots can be saved using the tmap_save(…) function - see ?tmap_save for a list of the possible file formats. For example,

World %>%
  filter(grepl("America", continent)) %>%
  tm_shape(.) +
  tm_polygons("HPI") -> mymap

tmap_mode("plot")
tmap_save(mymap, "mymap.pdf", height = 29.7, width = 21.0, units = "cm")

tmap_mode("view")
tmap_save(mymap, filename = "mymap.html")

browseURL("mymap.pdf")

browseURL("mymap.html")

To Do

To end this session, your task is to produce a map of Noise Pollution from roads in London, using data from the London Datastore

The first task if to download the data from the website where they are stored in a zipped format and then unzip them:

download.file("https://data.london.gov.uk/download/noise-pollution-in-london/1154a491-c68f-4451-9111-09bac975be99/Road_LAeq_16h_London.zip", "noise.zip")
unzip("noise.zip")

Next, we need to read in the data:

require(sf)
road_noise <- st_read("Road_LAeq_16h_London.shp")

Taking a quick look at the data, we can see it is of class

class(road_noise)

and contains the following attribute variable

names(st_drop_geometry(road_noise))

Because it is also a large data set,

nrow(road_noise)

we will just map the locations of the noisiest roads,

road_noise %>%
  filter(NoiseClass == ">=75.0") %>%
  filter(st_is_valid(.)) -> noise_subset

(The additional filter using st_is_valid is to avoid some geometric errors in the data)

Now produce a map of the noisiest roads to which you could add a North arrow and a scale bar, and save it in a suitable file format. The code block may take a little longer to run than before so please be patient.

noise_subset %>%
  # Replace this comment with your own code 

If you still have time…

If you have finished the exercises above and have some time left over, work through one of these (you can choose which):

1. (recommended, if you have a copy of the book)

Chapter 3 of An Introduction to R for Spatial Analysis & Mapping by Chris Brunsdon and Lex Comber. The code (but not text) for the chapter is here

Wickham H & Grolemund G, 2017, R for Data Science, Chapter 28: Graphics for Communication. https://r4ds.had.co.nz/graphics-for-communication.html

which is about ways to better communicate graphics.

browseURL("https://r4ds.had.co.nz/graphics-for-communication.html")

Leaflet for R. https://rstudio.github.io/leaflet/

browseURL("https://rstudio.github.io/leaflet/")

which offers another way to generate interactive maps.

Further Reading

Lovelace R, Nowosad J & Muenchow J, 2019, Geocomputation with R, Chapter 8: Making Maps with R. https://geocompr.robinlovelace.net/adv-map.html

Tennekes M, 2018, tmap: Thematic Maps in R. Journal of Statistical Software, 84 (6). https://www.jstatsoft.org/article/view/v084i06