Mapping with style

Robert Schlegel

Problem

A blank map is rarely enough; how do we add data layers?
How do we manage these multiple different colour palettes?
Are there any tricks for making complex maps?

Solution

We may import data, or use built-in options
There are best practices in ggplot2 for this
Yes, heaps!

Setup

library(tidyverse) # The base

library(marmap) # For downloading bathymetry data

Base map

Did anyone notice the x-axis in the default global map?

ggplot() +
  borders(fill = "grey70", colour = "black") +
  coord_equal()

Base map

This happens primarily because the Kamchatka peninsula goes over the date line
We can fix this with the following code chunk

map_global_fix <- map_data('world') %>% 
  rename(lon = long) %>% 
   # Why +2000?
  mutate(group = ifelse(lon > 180, group+2000, group),
         lon = ifelse(lon > 180, lon-360, lon))

Base map

Note that Antarctica still doesn’t look great

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # The default coordinate system, with specific limits
  coord_cartesian(xlim = c(-180, 180), ylim = c(-90, 90), expand = FALSE)

Bathymetry data

We don’t want to download every time we run
Better to download it once and save it

# Download bathy data
bathy_WA <-  getNOAA.bathy(lon1 = 111, lon2 = 117, 
                           # NB: smaller value first, i.e. more negative
                           lat1 = -36, lat2 = -19, 
                            # In degree minutes
                           resolution = 4)

# Convert to data.frame for use with ggplot2
bathy_WA_df <- fortify.bathy(bathy_WA) %>% 
  # Remove altimetry data
  filter(z <= 0) 

# Save
save(bathy_WA_df, file = "course_material/data/bathy_WA_df.RData")

Bathymetry data

But it is good to keep a record of the download
One option is ‘comment out’ lines of code
The keyboard short cut is ctrl+shift+c

# Download bathy data
# bathy_WA <-  getNOAA.bathy(lon1 = 111, lon2 = 117, 
#                            lat1 = -36, lat2 = -19,
#                            resolution = 4)

# Convert to data.frame for use with ggplot2
# bathy_WA_df <- fortify.bathy(bathy_WA)

# Save
# save(bathy_WA_df, file = "../data/bathy_WA_df.RData")

# Load
load("course_material/data/bathy_WA_df.RData")

Mapping bathymetry

A basic bathymetry map

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # Add 200 m contour
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z),
               breaks = c(-200), 
               linewidth = c(0.3), colour = "grey") +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Mapping bathymetry

Choose multiple contour lines

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # Add 200 and 2000 m contours
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z),
               breaks = c(-200, -2000), 
               linewidth = c(0.3), colour = "grey") +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Mapping bathymetry

Manually change contour colours

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # Assign colour per depth
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z),
               breaks = c(-200), linewidth = c(0.3), colour = "black") +
  # Assign colour per depth
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z),
               breaks = c(-2000), 
               linewidth = c(0.3), colour = "blue") +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Mapping bathymetry

Change contour colours with aes()

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # Rather use `aes()`
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z, colour = after_stat(level)),
               linewidth = c(0.3)) +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Mapping bathymetry

Discrete colour palettes for contours

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # Combine `aes()` and `breaks = c()` for more control
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z, colour = after_stat(level)),
               breaks = c(-50, -200, -1000, -2000), 
               linewidth = c(0.3)) +
  # Also change colour palette
  scale_colour_distiller(palette = "BuPu") + 
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Mapping bathymetry

Tidy up the appearance

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # create discrete factors
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z, colour = as.factor(after_stat(level))),
               breaks = c(-50, -200, -1000, -2000), 
               linewidth = c(0.3)) +
  # Use discrete palette
  scale_colour_brewer("Depth [m]", palette = "Set1", direction = -1) +  
  # Reverse legend order and make symbols thicker
  guides(color = guide_legend(reverse = TRUE, 
                              override.aes = list(linewidth = 5))) +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Ocean temperature

We will cover the downloading of data more on Thursday (Day 7)
For now we will use pre-downloaded data from NOAA

# Load sea surface temperatures for 2000-01-01
load("course_material/data/OISST_2000.RData")

Ocean temperature

Plot Western Australia

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_polygon(aes(group = group), 
               colour = "black", fill = "grey60") +
  geom_raster(data = OISST_2000, aes(fill = temp)) +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Combining layers

Remember that the order of the code is how the layers are plotted

ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  # First layer
  geom_raster(data = OISST_2000, aes(fill = temp)) + 
  # Second layer
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  # Third layer
  geom_contour(data = bathy_WA_df, 
               aes(x = x, y = y, z = z, colour = as.factor(after_stat(level))), 
               breaks = c(-50, -200, -1000, -2000), 
               linewidth = c(0.3)) +
  guides(color = guide_legend(reverse = TRUE, 
                              override.aes = list(linewidth = 5))) + 
  scale_fill_viridis_c("Temperature [°C]") +
  scale_colour_brewer("Depth [m]", palette = "BuPu") +
  coord_cartesian(xlim = c(111, 117), 
                  ylim = c(-36, -19), expand = FALSE)

Combining layers

Finishing touches

An example of a complete map code chunk

final_map <- ggplot(data = map_global_fix, aes(x = lon, y = lat)) +
  geom_raster(data = OISST_2000, aes(fill = temp)) +
  geom_polygon(aes(group = group), colour = "black", fill = "grey60") +
  geom_contour(data = bathy_WA_df,
               aes(x = x, y = y, z = z, 
                   colour = as.factor(after_stat(level))), 
               breaks = c(-50, -200, -1000, -2000), 
               linewidth = c(0.3)) +
  guides(color = guide_legend(reverse = TRUE, 
                              override.aes = list(linewidth = 5))) + 
  scale_fill_viridis_c("Temperature [°C]") +
  scale_colour_brewer("Depth [m]", palette = "BuPu") +
  coord_cartesian(xlim = c(111, 117), ylim = c(-36, -19), expand = FALSE) +
  # Put x axis labels on top of figure and assign °E
  scale_x_continuous(position = "top", 
                     breaks = c(112, 114, 116), 
                     labels = c("112°E", "114°E", "116°E")) + 
  # Put y axis labels on right of figure and assign °S
  scale_y_continuous(position = "right",
                     breaks = c(-34, -28, -22), 
                     labels = c("34°S", "28°C", "22°C")) +
   # Remove the axis label text
  theme(axis.title = element_blank(),
        # Add black border
        panel.border = element_rect(fill = NA, colour = "black"), 
        # Change text size in legend
        legend.text = element_text(size = 7), 
        # Change legend title text size
        legend.title = element_text(size = 7), 
        # Change size of legend
        legend.key.height = unit(0.5, "cm"),
         # Add legend background
        legend.background = element_rect(fill = "white", colour = "black"),
        # Change position of legend
        legend.position = c(0.9, 0.5)
        )

The end result

The cleaned up map of Western Ausralia. Resplendent with ocean temperatures (°C) for 2000-01-01.

Don’t forget to save!

ggsave(plot = final_map, "figures/map_complete.pdf", height = 6, width = 9)