SPUR R Lecture - Lecture 3

J. Kavanagh

2022-07-14

Introduction

In this series of slides we will be exploring how to re-create historical visualisations using the ggplot2 package and basic mapping using tmap

load('SPUR_2.RData')

First we are going to re-create Minard’s map

# Load the following data via your console
data(Minard.troops)
data(Minard.cities)
data(Minard.temp)
glimpse(Minard.troops)
## Rows: 51
## Columns: 5
## $ long      <dbl> 24.0, 24.5, 25.5, 26.0, 27.0, 28.0, 28.5, 29.0, 30.0, 30.3, …
## $ lat       <dbl> 54.9, 55.0, 54.5, 54.7, 54.8, 54.9, 55.0, 55.1, 55.2, 55.3, …
## $ survivors <int> 340000, 340000, 340000, 320000, 300000, 280000, 240000, 2100…
## $ direction <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, R, R, R, R, …
## $ group     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
glimpse(Minard.troops)
## Rows: 51
## Columns: 5
## $ long      <dbl> 24.0, 24.5, 25.5, 26.0, 27.0, 28.0, 28.5, 29.0, 30.0, 30.3, …
## $ lat       <dbl> 54.9, 55.0, 54.5, 54.7, 54.8, 54.9, 55.0, 55.1, 55.2, 55.3, …
## $ survivors <int> 340000, 340000, 340000, 320000, 300000, 280000, 240000, 2100…
## $ direction <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, R, R, R, R, …
## $ group     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
glimpse(Minard.troops)
## Rows: 51
## Columns: 5
## $ long      <dbl> 24.0, 24.5, 25.5, 26.0, 27.0, 28.0, 28.5, 29.0, 30.0, 30.3, …
## $ lat       <dbl> 54.9, 55.0, 54.5, 54.7, 54.8, 54.9, 55.0, 55.1, 55.2, 55.3, …
## $ survivors <int> 340000, 340000, 340000, 320000, 300000, 280000, 240000, 2100…
## $ direction <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, R, R, R, R, …
## $ group     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …

Start creating the first Troops layer for the map

# These are required, be sure to run this code at this time
require(scales)
## Loading required package: scales
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
require(gridExtra)
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

The following code will create a number of new objects.

# This creates the base layer of troop movements
Minard.troops %>% ggplot(aes(long, lat)) +
geom_path(aes(size = survivors, colour = direction, group = group),
              lineend = "round", linejoin = "round") -> plot_troops
# Check the result
plot_troops

Next create the Cities layer

This is more for the labels and some of the breaks used, it seems simple but is very important overall

geom_text(aes(label = city), size = 4, data = Minard.cities) -> plot_cities
breaks <- c(1, 2, 3) * 10^5

Create the joined Troops and Cities layer

This creates a joined map showing the Troops and Cities layers

plot_troops + plot_cities +
scale_size("Survivors", range = c(1, 10), 
                breaks = breaks, labels = scales::comma(breaks)) +
scale_color_manual("Direction", 
                    values = c("grey50", "red"), 
                    labels=c("Advance", "Retreat")) +
coord_cartesian(xlim = c(24, 38)) +
xlab(NULL) + 
ylab("Latitude") + 
ggtitle("Napoleon's March on Moscow") +
theme_bw() +
theme(legend.position=c(.8, .2), legend.box="horizontal") -> plot_minard
# Check the result
plot_minard

Create the Temp layer

Minard.temp %>% ggplot(aes(long, temp)) +
geom_path(color="grey", size=1.5) +
geom_point(size=2) +
geom_text(aes(label=date)) +
xlab("Longitude") + ylab("Temperature") +
coord_cartesian(xlim = c(24, 38)) + theme_bw() -> plot_temp
# Check the result
plot_temp
## Warning: Removed 1 rows containing missing values (geom_text).

Display all the layers with an adjusted re-scaled plot window

grid.arrange(plot_minard, plot_temp, nrow=2, heights=c(3,1))
## Warning: Removed 1 rows containing missing values (geom_text).

Introduction to tmap & sf

# load a map of Ireland 
st_read('./counties/counties.shp') -> cty
## Reading layer `counties' from data source 
##   `/Users/jackkavanagh/Dropbox/R_SPUR/counties/counties.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -10.66221 ymin: 51.41943 xmax: -5.995289 ymax: 55.44661
## Geodetic CRS:  WGS 84
# Use the qtm() function to check that it imported correctly
cty %>% qtm

Adding a variable to the map

cty %>% qtm(fill="POPDENSITY")

Changing the background

cty %>% qtm(fill="POPDENSITY") -> cty_dens

# Add new features
cty_dens + tm_compass() + tm_style("natural")
## Note that tm_style("natural") resets all options set with tm_layout, tm_view, tm_format, or tm_legend. It is therefore recommended to place the tm_style element prior to the other tm_layout/tm_view/tm_format/tm_legend elements.

Adjusting for colour blind

cty %>% qtm(fill="POPDENSITY") -> cty_dens

# Add new features
cty_dens + tm_scale_bar() + tm_compass() + tm_style("col_blind")
## Note that tm_style("col_blind") resets all options set with tm_layout, tm_view, tm_format, or tm_legend. It is therefore recommended to place the tm_style element prior to the other tm_layout/tm_view/tm_format/tm_legend elements.

Creating a new shapefile

# You can change the type of projection 
qtm(cty,projection="+init=epsg:29900",fill="POPDENSITY") + tm_scale_bar()
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.

# Create a new shapefile with this projection
cty %>% set_projection("+init=epsg:29900") -> cty2

Labelling and finalising the map

tm_shape(cty2) + 
tm_style("classic") + 
tm_fill(col="POPDENSITY",title='Population Density') + 
tm_borders() + 
tm_compass(size=4) 

Cartograms

# This creates a new cartogram of Ireland using the population from the census
cty2 %>% cartogram(weight="TOTAL2011") -> cart
## 
## Please use cartogram_cont() instead of cartogram().
## Mean size error for iteration 1: 7.22625035681413
## Mean size error for iteration 2: 3.90812705063481
## Mean size error for iteration 3: 2.42846635451752
## Mean size error for iteration 4: 1.76512515998386
## Mean size error for iteration 5: 1.43479584410048
## Mean size error for iteration 6: 1.25671605744483
## Mean size error for iteration 7: 1.15695321419564
## Mean size error for iteration 8: 1.09928846404667
## Mean size error for iteration 9: 1.06498099166659
## Mean size error for iteration 10: 1.04385824473235
## Mean size error for iteration 11: 1.03016175221771
## Mean size error for iteration 12: 1.02100187189164
## Mean size error for iteration 13: 1.01481601671164
## Mean size error for iteration 14: 1.01056748842113
## Mean size error for iteration 15: 1.00761925463118
# Check the results
cart %>% qtm

Using a cartogram to visualise two variables

tm_shape(cart) + 
tm_fill(col="PCVAC20111",title='% Vacant Houses') + 
tm_borders()

Facet or Comparative Maps

tm_shape(cty2) + 
tm_fill(col=c("PCVAC20111","POPDENSITY"),
        title=c('% Vacant Houses','Population Density')) + 
tm_borders() 

Creating a ‘live’ map

tm_shape(cty) + 
tm_fill(col="PCVAC20111",
        title='% Vacant Houses',
        style='jenks') + 
tm_borders() -> cmap

tmap_leaflet(cmap)

Live mapping with greater adjustments

# First set tmap to live view
tmap_mode('view')
## tmap mode set to interactive viewing
# Recreate the cmap 
tm_shape(cty) + 
tm_fill(col="PCVAC20111",
        title='% Vacant Houses',
        style='jenks',
        alpha = 0.5) + 
tm_borders() -> cmap

# The results
cmap

Adding two maps using tmap

# Import the map of Gaeltacht areas
st_read('./Census2011_Gaeltacht/Census2011_Gaeltacht.shp') -> gael
## Reading layer `Census2011_Gaeltacht' from data source 
##   `/Users/jackkavanagh/Dropbox/R_SPUR/Census2011_Gaeltacht/Census2011_Gaeltacht.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 17481.2 ymin: 19589.93 xmax: 285663.3 ymax: 448059.1
## Projected CRS: TM65 / Irish Grid
# Create a new object with both maps
tm_shape(cty) + 
tm_fill(col="PCVAC20111",
        title='% Vacant Houses',
        style='jenks') +
tm_borders() + 
tm_shape(gael) + 
tm_fill(col='darkblue',alpha=0.6) -> cmap

cmap + tm_view(alpha = 0.5)