library(sf)  
library(terra)  
library(tidyverse)  
library(spData)

#install.packages("spDataLarge", repos = "https://nowosad.r-universe.dev")
library(spDataLarge)
## Warning: package 'spDataLarge' was built under R version 4.3.3
library(magrittr)

library(tmap)    # for static and interactive maps
library(leaflet) # for interactive maps
library(ggplot2) # tidyverse data visualization package

Making Maps with R

Time to make some maps. We’ll focus on the ‘tmap’ package here, but there are a number of packages dedicated to making maps. We’ll also introduce ‘leaflet’ and the mapping capabilites of ‘ggplot2’. First, we will go over static maps and then move on to dynamic and interactive maps.

Static Maps

Static maps are the most common result when visualizing spatial data. Base R’s plot() can be used for static maps as we’ve demonstrated throughout the workshop. However, ‘tmap’ offers more options for aesthetics and follows a similar syntax with ggplot.

tmap

Like ggplot, tmap follows a ‘grammar of graphics’ (the gg in ggplot), where input data is separated from the aesthetics. Input data is defined using tm_shape() and aesthetics are layered on using a wide variety of functions such as tm_fill() and tm_dots().

# Add fill layer to nz shape
tm_shape(nz) +
  tm_fill() 

# Add border layer to nz shape
tm_shape(nz) +
  tm_borders() 

# Add fill and border layers to nz shape
tm_shape(nz) +
  tm_fill() +
  tm_borders()

Play around with tmap. Type help(“tmap-element”) into your console for a full list of functions.

Map objects

A useful feature of tmap is its ability to store objects representing maps.

map_nz <- tm_shape(nz) + tm_polygons()
class(map_nz)
## [1] "tmap"

Storing maps as map objects is useful for adding aesthetic layers and new spatial objects on top of each other. You can now add layers to ‘map_nz’ by simply using the + operator. If you add + tm_shape() then you can add another spatial object.

nz_elev = rast(system.file("raster/nz_elev.tif", package = "spDataLarge"))

map_nz1 = map_nz +
  tm_shape(nz_elev) + tm_raster(alpha = 0.7)
nz_water = st_union(nz) %>% st_buffer(22200) %>% 
  st_cast(to = "LINESTRING")
map_nz2 = map_nz1 +
  tm_shape(nz_water) + tm_lines()
map_nz3 = map_nz2 +
  tm_shape(nz_height) + tm_dots()
tmap_arrange(map_nz1, map_nz2, map_nz3)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## Some legend labels were too wide. These labels have been resized to 0.63, 0.55, 0.55, 0.55, 0.55, 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## Some legend labels were too wide. These labels have been resized to 0.63, 0.55, 0.55, 0.55, 0.55, 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
## Some legend labels were too wide. These labels have been resized to 0.63, 0.55, 0.55, 0.55, 0.55, 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Aesthetics

Like ggplot, the default aesthetics of each layer can be changed. Unlike ggplot, tmap does not use the aes() function but layers can accept arguments that are fixed values or variable fields (based on column names).

Here is an example of fixed values

ma1 = tm_shape(nz) + tm_fill(col = "red")
ma2 = tm_shape(nz) + tm_fill(col = "red", alpha = 0.3)
ma3 = tm_shape(nz) + tm_borders(col = "blue")
ma4 = tm_shape(nz) + tm_borders(lwd = 3)
ma5 = tm_shape(nz) + tm_borders(lty = 2)
ma6 = tm_shape(nz) + tm_fill(col = "red", alpha = 0.3) +
  tm_borders(col = "blue", lwd = 3, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)

And here is an example of how you would use values from a variable. The argument will take the character string associated with the column you are interested in.

tm_shape(nz) + tm_fill(col = "Land_area")

Finally, you can add nicer titles. expression() allows for superscripts

legend_title = expression("Area (km"^2*")")
map_nza = tm_shape(nz) +
  tm_fill(col = "Land_area", title = legend_title) + tm_borders()
map_nza

Color settings

Colors are useful. Here are a few ways to manipulate color in your maps. Play around and figure out what each of these arguments (breaks, n, palette) do.

col1 <- tm_shape(nz) + tm_polygons(col = "Median_income")
breaks <- c(0, 3, 4, 5) * 10000
col2 <- tm_shape(nz) + tm_polygons(col = "Median_income", breaks = breaks)
col3 <- tm_shape(nz) + tm_polygons(col = "Median_income", n = 10)
col4 <- tm_shape(nz) + tm_polygons(col = "Median_income", palette = "BuGn")
tmap_arrange(col1, col2, col3, col4)

tmap also allows for setting breaks using the style argument rather than manual setting breaks. Try out a few of the different style options: “pretty”, “equal”, “quantile”, “jenks”, “cont”, “cat”

tm_shape(nz) + tm_polygons(col = "Median_income", style = "pretty")

Palettes define the color range used and can be critical for interpretability of your figures. Palettes can be categorical, sequential, or diverging and each of these serve different purposes. To adjust the palette, you can use the palette argument in tmap.

tm_shape(nz) + tm_polygons("Population", palette = "Blues")
## Some legend labels were too wide. These labels have been resized to 0.64, 0.59, 0.59. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

tm_shape(nz) + tm_polygons("Population", palette = "YlOrBr")
## Some legend labels were too wide. These labels have been resized to 0.64, 0.59, 0.59. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Please reference section 8.2.4 of the online textbook to learn more about good habits for choosing palettes and when it is appropriate to use each.

Layouts

The map layout refers to the combination of all map elements into a cohesive map. Map elements include among others the objects to be mapped, the title, the scale bar, margins and aspect ratios.

These additional elements can have their own functions for easy layering onto your maps.

map_nz + 
  tm_compass(type = "8star", position = c("left", "top")) +
  tm_scale_bar(breaks = c(0, 100, 200), text.size = 1)

tm_layout allows for more customization of layout elements

layout1 <- map_nz + tm_layout(title = "New Zealand")
layout2 <- map_nz + tm_layout(scale = 5)
layout3 <- map_nz + tm_layout(bg.color = "lightblue")
layout4 <- map_nz + tm_layout(frame = FALSE)
tmap_arrange(layout1, layout2, layout3, layout4)

tmap also provides predetermined style options for layout that you can access using the style argument.

map_nza + tm_style("bw")

map_nza + tm_style("classic")

map_nza + tm_style("cobalt")

map_nza + tm_style("col_blind")

Please reference section 8.2.5 of the online textbook for more information on what you can do to the layout.

Faceted maps

Faceted maps are useful for representing how data change in relationship to another variable, such as time. Here is an example of population size of large cities changing over time (including predictions of the future).

Play around with the number of years included in this figure and the arguments for the tm_facets() function.

urb_1970_2030 = urban_agglomerations %>% 
  filter(year %in% c(1970, 1990, 2010, 2030))

tm_shape(world) +
  tm_polygons() +
  tm_shape(urb_1970_2030) +
  tm_symbols(col = "black", border.col = "white", size = "population_millions") +
  tm_facets(by = "year", nrow = 2, free.coords = FALSE)

Inset maps

An inset map is a smaller map rendered within or next to the main map and can be useful for adding more detail to your maps.

As an example, we can start by making a new spatial object

nz_region = st_bbox(c(xmin = 1340000, xmax = 1450000,
                      ymin = 5130000, ymax = 5210000),
                    crs = st_crs(nz_height)) %>% 
  st_as_sfc()

Next, we can make the elevation map

nz_height_map = tm_shape(nz_elev, bbox = nz_region) +
  tm_raster(style = "cont", palette = "YlGn", legend.show = TRUE) +
  tm_shape(nz_height) + tm_symbols(shape = 2, col = "red", size = 1) +
  tm_scale_bar(position = c("left", "bottom"))

Now we can make the inset map

nz_map = tm_shape(nz) + tm_polygons() +
  tm_shape(nz_height) + tm_symbols(shape = 2, col = "red", size = 0.1) + 
  tm_shape(nz_region) + tm_borders(lwd = 3)

And put it all together

library(grid)
## 
## Attaching package: 'grid'
## The following object is masked from 'package:terra':
## 
##     depth
nz_height_map
## stars object downsampled to 877 by 1140 cells. See tm_shape manual (argument raster.downsample)
print(nz_map, vp = viewport(0.8, 0.27, width = 0.5, height = 0.5))

Another purpose of inset maps is to include create a single map of non-contiguous areas such as a map of the US.

First, we can get a map of the contiguous US

us_states_map = tm_shape(us_states, projection = 2163) + tm_polygons() + 
  tm_layout(frame = FALSE)

Next, we can define maps for hawaii and alaska

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)

Finally, we can print these maps and define where we want them to go on the map

us_states_map
print(hawaii_map, vp = grid::viewport(0.35, 0.1, width = 0.2, height = 0.1))
print(alaska_map, vp = grid::viewport(0.15, 0.15, width = 0.3, height = 0.3))

Animated maps

Rather than using faceted plots, animated maps can help visualize the same data using less space and in a more intuitive way.

urb_anim = tm_shape(world) + tm_polygons() + 
  tm_shape(urban_agglomerations) + tm_dots(size = "population_millions") +
  tm_facets(along = "year", free.coords = FALSE)

This will create the animation and save it onto your computer. You’ll want to designate the filepath so that you can easily find it. This will default into your current working directory. If you are using MacOS then you can view it as a gif by opening it using cmd+y

#tmap_animation(urb_anim, filename = "urb_anim.gif", delay = 25)

Interactive maps

The development of packages for interactive maps in R has made interactivity extremely accessible. ‘tmap’ makes this easy by allowing most static maps easily transferrable to interactive maps using tmap_mode()

See below for an example of turning ‘map_nz’ into a simple interactive map using tmap_mode(“view”)

tmap_mode("view")
## tmap mode set to interactive viewing
map_nz

Now that tmap_mode() has been set to “view”, all of the maps generated with tmap will show up as interactive maps

map_nza

You can return to static maps by using tmap_mode(“plot”)

tmap_mode("plot")
## tmap mode set to plotting
map_nza

Another package that allows for quick generation of interactive maps is ‘mapview’

library(mapview)
trails %>%
  st_transform(st_crs(franconia)) %>%
  st_intersection(franconia[franconia$district == "Oberfranken", ]) %>%
  st_collection_extract("LINE") %>%
  mapview(color = "red", lwd = 3, layer.name = "trails") +
  mapview(franconia, zcol = "district", burst = TRUE) +
  breweries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

leaflet is another option that is more low-level allowing for more customization

pal = colorNumeric("RdYlBu", domain = cycle_hire$nbikes)
leaflet(data = cycle_hire) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(col = ~pal(nbikes), opacity = 0.9) %>% 
  addPolygons(data = lnd, fill = FALSE) %>% 
  addLegend(pal = pal, values = ~nbikes) %>% 
  setView(lng = -0.1, 51.5, zoom = 12) %>% 
  addMiniMap()

More mapping options

To see more mapping options, you can check out the rest of the material in Chapter 8 of the online textbook. For example, this section goes into using leaflet and shiny to make interactive web apps.

Exercises

These exercises rely on a new object, africa. Create it using the world and worldbank_df datasets from the spData package as follows:

africa = world %>% 
  filter(continent == "Africa", !is.na(iso_a2)) %>% 
  left_join(worldbank_df, by = "iso_a2") %>% 
  dplyr::select(name, subregion, gdpPercap, HDI, pop_growth) %>% 
  st_transform("+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25")

We will also use zion and nlcd datasets from spDataLarge:

zion = read_sf((system.file("vector/zion.gpkg", package = "spDataLarge")))
nlcd = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))

  1. Create a map showing the geographic distribution of the Human Development Index (HDI) across Africa with base graphics (hint: use plot()) and tmap packages (hint: use tm_shape(africa) + ...).

    • Name two advantages of each based on the experience.
    • Name three other mapping packages and an advantage of each.
    • Bonus: create three more maps of Africa using these three packages.

  2. Extend the tmap created for the previous exercise so the legend has three bins: “High” (HDI above 0.7), “Medium” (HDI between 0.55 and 0.7) and “Low” (HDI below 0.55).

    • Bonus: improve the map aesthetics, for example by changing the legend title, class labels and color palette.

  3. Represent africa’s subregions on the map. Change the default color palette and legend title. Next, combine this map and the map created in the previous exercise into a single plot.

  4. Create a land cover map of the Zion National Park.

    • Change the default colors to match your perception of the land cover categories
    • Add a scale bar and north arrow and change the position of both to improve the map’s aesthetic appeal
    • Bonus: Add an inset map of Zion National Park’s location in the context of the Utah state. (Hint: an object representing Utah can be subset from the us_states dataset.)

  5. Create facet maps of countries in Eastern Africa:

    • With one facet showing HDI and the other representing population growth (hint: using variables HDI and pop_growth, respectively)
    • With a ‘small multiple’ per country

  6. Building on the previous facet map examples, create animated maps of East Africa:

    • Showing first the spatial distribution of HDI scores then population growth
    • Showing each country in order

  7. Create an interactive map of Africa:

    • With tmap
    • With mapview
    • With leaflet
    • Bonus: For each approach, add a legend (if not automatically provided) and a scale bar