09-mapping.Rmd

# (PART) Advanced methods {-}

# Making maps with R {#adv-map}

## Prerequisites {-}

- This chapter requires the following packages that we have already been using:

```{r, message = FALSE} 
library(sf)
library(spData)
library(tidyverse)
```

- In addition it will make use of the following visualisation packages:

```{r}
library(leaflet)
library(tmap)
# library(mapview)
# library(shiny)
```

## Introduction

The most satisfying part of many projects involving geocomputation is producing and communicating the results in the form of maps.
Map making --- the art of Cartography --- is an ancient skill that involves precision, consideration of the map-reader and often an element of creativity.

Maps have been used for several thousand years for a wide variety of purposes.
From asserting control real estate in ancient Mesapotamia to illustrating a scientific worldview in Ptolomy's Greek masterpiece *Geography* nearly 2,000 years ago [@talbert_ancient_2014], they have historically been out of reach for everyday people.
The advent of modern computing has changed this, however.
Open data and software has played a major role in the democratization of map making.
This chapter teaches a range of methods to amplify the clarity, reach and attractiveness of your maps.
From a research perspective maps are often the best way to present the results of geocomputational research.
From the perspective of public engagement and policy impact, a map may be the best way to get a message across.
Map making is therefore a critical part of geocomputation and its emphasis on changing the world outlined in Chapter \@ref(intro).

We have already seen how to create maps using methods building on the substantial graphical capabilities of base R in section \@ref(basic-map).
It is possible to go far with these methods, as illustrated in Chapter [14](https://www.stat.auckland.ac.nz/~paul/RG2e/chapter14.html) of *R Graphics* [@murrell_r_2016].
The focus of this chapter, however, is creating maps using packages that are dedicated to this purpose.
In many cases the simplicity and speed of maps created using `plot()` commands will be sufficient, especially when using R interactively to understand data quickly.
In some case you may benefit from the flexibility of low-level control over graphical parameters offered by `plot()` and functions provided by the **grid** package.
If your aim is to produce publication-quality or high-impact maps in minimal time and code, this chapter is for you.

<!-- info about relation between efficiency and editability -->
<!-- intro to the chapter structure -->

## Static maps

### Basics

<!-- I'm not sure if we can fill it all in the book, but it could be worth to try these three pacakges on the same problems/data -->
<!-- base plots (one example) -->
<!-- ggplots  (geom_sf, coord_sf; one example)-->
<!-- tmap (more than one example) -->

### Map styling

<!-- 
- lines widths, polygon lines vs polygon areas, etc.
- alpha colors
- joining points/lines/polygons
- joining vector/raster data
- titles, labels
- satellite tiles in the background
- choropleth maps
- north arrow
- scale bar
-->

### Class intervals

<!--
- color palettes (viridis, etc.)
- class intervals ("cat", "fixed", "sd", "equal", "pretty", "quantile", "kmeans", "hclust", "bclust", "fisher", and "jenks")
- bicolor maps?
- categorical rasters?
-->

### Facetted maps

Facetted maps are a common and potentially effective way of visualizing spatial relationships that are more complex that a single relationship.
The population of cities at one moment in time can be represented easily on a single map, for example by making the size of symbols variable for each city depending on population.
However, to represent the populations of cities at multiple moments in time requires an *extra dimension*.
This could be added by an additional *aesthetic* such as colour but this risks cluttering the map because it will involve multiple overlapping points (cities do not tend to move over time!).

Facetted maps overcome this issue by creating 'small multiples' of spatial information.
Typically, this will involve plotting the same geometry data multiple times, once for each variable present in the attribute data (this is the default plotting method for `sf` objects, as we saw in Chapter 2).
However, facets can also represent shifting geometries, e.g. as the spatial distribution of a point pattern evolves over time.
This use case of facetted plot is illustrated in Figure \@ref(fig:urban-facet).

```{r urban-facet, fig.cap="Facetted map showing the top 30 largest 'urban agglomerations' from 1950 to 2030 based on population projects by the United Nations"}
tm_shape(world) +
  tm_polygons() + 
  tm_shape(urban_agglomerations) +
  tm_dots(size = "population_millions") +
  tm_facets(by = "year")
```

### Inset maps

The role of spatial visualizations is to provide an information in the most accessible form. 
This is often done with the use of inset maps.
They could be used to focus on a smaller area in more detail<!--fig ref1--> or to bring some non-contiguous regions closer to ease a comparision<!--fig ref2-->.

<!-- intro  -->
<!-- used when: a/ non-contigous regions, b/ focus on small region while showing the large one -->
<!-- find the best distribution of insets -->

<!-- example1: classic inset map -->

<!-- example2: us states + hawaii + alaska -->


<!-- find the best projection for each individual inset. -->

<!-- reproj -->
```{r}
us_states_map = tm_shape(us_states, projection = 2163) +
  tm_polygons() + 
  tm_layout(frame = FALSE)
```

<!-- create separate maps -->
```{r}
hawaii_map = tm_shape(hawaii) +
  tm_polygons() + 
  tm_layout(title = "Hawaii", frame = FALSE, bg.color = NA, 
            title.position = c("left", "bottom"))
alaska_map = tm_shape(alaska) +
  tm_polygons() + 
  tm_layout(title = "Alaska", frame = FALSE, bg.color = NA)
```

<!-- create an inset map -->
```{r}
library(grid)
us_states_map
print(hawaii_map, vp = viewport(x = 0.4, y = 0.1, width = 0.2, height = 0.1))
print(alaska_map, vp = viewport(x = 0.15, y = 0.15, width = 0.3, height = 0.3))
```

<!-- extended info about using tm_layout to show legend in main plot and remove it in the others -->

<!-- example3: mix of maps and plots? -->
<!-- of maybe only mention this -->
<!-- https://gist.github.com/walkerke/93bfe80bb7735aa6265a61013eaed3fa -->

## Animations

Animated maps can be useful for communicating how spatial phenomena shift over time.
An advantage of facetted plots are that they can be printed, but the approach has disadvantages:
faceted maps can become very small with more than ~9 maps in one figure, and it can be hard to see the spatial relationships between each facet when each map is on a different part of the page!
Furthermore, with the increasing proportion of communication that happens via digital screens, the disadvantage that animations cannot be printed is diminished.
You can always link readers to a web-page containing an animated (or interactive) version of a printed map to help make it come alive.

Figure \@ref(fig:urban-animated) is a simple example of the benefits of an animated map.
Unlike the facetted plot presented in the previous section, it does not squeeze all 17 for them all to be displayed simultaneously (see the book's website for the animated version).

```{r urban-animated, fig.cap="Animated map showing the top 30 largest 'urban agglomerations' from 1950 to 2030 based on population projects by the United Nations."}
knitr::include_graphics("figures/urban-animated.gif")
```

```{r}
m = tm_shape(world) + 
  tm_polygons() +
  tm_shape(urban_agglomerations) +
  tm_dots(size = "population_millions") +
  tm_facets(by = "year", nrow = 1, ncol = 1) 
```

```{r, echo=FALSE, eval=FALSE}
source("code/09-urban-animation.R")
```

<!-- Robin, check the differences between by="year" and along="year" -->

```{r ani_data_prep, eval=FALSE}
options(scipen = 999)

wb_data_create = function(indicator, our_name, year, ...){
  df = wb(indicator = indicator, startdate = year, enddate = year, ...) %>%
    as_data_frame() %>%
    select(iso_a2=iso2c, value) %>%
    mutate(indicator = our_name) %>%
    spread(indicator, value)
  return(df)
}

data_lifeExp = seq(1963, 2013, by=5) %>%
  set_names(.) %>%
  map_df(~wb_data_create(.x, indicator = "SP.DYN.LE00.IN",
                   our_name = "lifeExp",
                   country = "countries_only"), .id='year') %>%
  spread(year, lifeExp)

world_sf_temporal = ne_countries(returnclass = 'sf') %>%
  left_join(., data_lifeExp, by = c('iso_a2')) %>%
  mutate(area_km2 = set_units(st_area(.), km^2)) %>%
  select(iso_a2, name_long, continent, region_un, subregion, type, area_km2, `1963`:`2013`) %>%
  gather(year, lifeExp, `1963`:`2013`)
```

```{r animation1, eval=FALSE}
m1 = tm_shape(world_sf_temporal) + 
  tm_polygons("lifeExp") +
  tm_facets(by = "year", nrow = 1, ncol = 1, drop.units = TRUE)

animation_tmap(m1, filename = "figures/11-lifeExp_animation.gif", width = 2000, height = 600, delay = 40)
```

```{r animation2, eval=FALSE}
world_sf_temporal2 = filter(world_sf_temporal, continent == "South America")
m2 = tm_shape(world_sf_temporal2) +
  tm_polygons("lifeExp", n = 12) +
  tm_facets(by = "name_long", along = "year", drop.units = TRUE, free.coords = TRUE)
animation_tmap(m2, filename = "figures/11-lifeExp_sa_animation.gif", width = 1600, height = 1000, delay = 40)
```

## Interactive maps

<!-- leaflet -->
<!-- leaflet plugins -->
<!-- mapview -->
<!-- mapview raster support! -->
<!-- are we going to focus on leaflet or mapview? -->

## Pseudo (unusual) maps 
<!--we need better name -->

<!-- geofacet -->
<!-- e.g. population lines -->
<!-- point maps -->
<!-- hexbin -->
<!-- cartogram - https://github.com/sjewo/cartogram -->
<!-- shiny? -->