9.1 Creating a thematic map

In this post, we go through all the steps required to produce a complete good-looking map. We will see how to add a title, a legend, a scale, axis and a North arrow and choose a good color palette. To do so, we will use the Quebec province as our sampled area.

So, first, we import 2 vector layers readily available in R, the Canadian provincial boundaries and USA country boundary.

library(sf)
library(raster)

can1 <- getData("GADM", country = "CAN", level = 1, path = "data")

#R>  Warning in getData("GADM", country = "CAN", level = 1, path = "data"): getData will be removed in a future version of raster
#R>  . Please use the geodata package instead

usa0 <- getData("GADM", country = "USA", level = 0, path = "data")

#R>  Warning in getData("GADM", country = "USA", level = 0, path = "data"): getData will be removed in a future version of raster
#R>  . Please use the geodata package instead

Transform to sf objects to facilitate manipulation

can1_sf <- st_as_sf(can1)
usa0_sf <- st_as_sf(usa0)

Retrieve Quebec polygon and surrounding provinces

qc <- can1_sf[can1_sf$NAME_1 == "Québec", ]

qc_neigh <- can1_sf[can1_sf$NAME_1 %in% c("Québec", "Ontario", "Nova Scotia", "New Brunswick", "Newfoundland and Labrador"), ]

It would take a while to plot because there is a lot of unnecessary details, so we can simplify the shape of the polygons using st_simplify().

usa0_simple <- st_simplify(usa0_sf, dTolerance = .05, preserveTopology = FALSE)

#R>  Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
#R>  argument preserveTopology cannot be set to FALSE when working with ellipsoidal
#R>  coordinates since the algorithm behind st_simplify always preserves topological
#R>  relationships

qc_simple <- st_simplify(qc, dTolerance = .01, preserveTopology = FALSE)

#R>  Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
#R>  argument preserveTopology cannot be set to FALSE when working with ellipsoidal
#R>  coordinates since the algorithm behind st_simplify always preserves topological
#R>  relationships

qc_neigh_simple <- st_simplify(qc_neigh, dTolerance = .01, preserveTopology = FALSE)

#R>  Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
#R>  argument preserveTopology cannot be set to FALSE when working with ellipsoidal
#R>  coordinates since the algorithm behind st_simplify always preserves topological
#R>  relationships

The warning stating that st_simplify does not correctly simplify longitude/latitude data and it should be in decimal degrees. It’s important in some spatial operations involving distances, but we can ignore it for the purpose of mapping.

Let’s look at what we have:

plot(st_geometry(qc_simple), col = "grey85")
plot(st_geometry(qc_neigh_simple),  col = "grey45", add = TRUE)
plot(st_geometry(usa0_simple), col = "grey25", add = TRUE)

Let’s say we sampled vegetation in 100 plots in Quebec; we now want to plot them with points proportional to their species richness. We will now create a data frame containing coordinates and random values from 5 to 50.

# Sample random points from our study area
sample_pts <- st_sample(x = qc_simple, size = 100)

# Create an attribute of fake species richness
sample_richness <- sample(x = 5:50, size = 100, replace = TRUE)