Last updated: 2019-06-10

workflowr checks: (Click a bullet for more information)
  • R Markdown file: up-to-date

    Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

  • Environment: empty

    Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

  • Seed: set.seed(20190513)

    The command set.seed(20190513) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

  • Session information: recorded

    Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

  • Repository version: 028d3cc

    Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

    Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
    
    Ignored files:
        Ignored:    .Rhistory
        Ignored:    .Rproj.user/
        Ignored:    data/NAPA_clim_U.Rda
        Ignored:    data/NAPA_clim_V.Rda
        Ignored:    data/NAPA_clim_W.Rda
        Ignored:    data/NAPA_clim_emp_ice.Rda
        Ignored:    data/NAPA_clim_emp_oce.Rda
        Ignored:    data/NAPA_clim_fmmflx.Rda
        Ignored:    data/NAPA_clim_mldkz5.Rda
        Ignored:    data/NAPA_clim_mldr10_1.Rda
        Ignored:    data/NAPA_clim_qemp_oce.Rda
        Ignored:    data/NAPA_clim_qla_oce.Rda
        Ignored:    data/NAPA_clim_qns.Rda
        Ignored:    data/NAPA_clim_qsb_oce.Rda
        Ignored:    data/NAPA_clim_qt.Rda
        Ignored:    data/NAPA_clim_runoffs.Rda
        Ignored:    data/NAPA_clim_ssh.Rda
        Ignored:    data/NAPA_clim_sss.Rda
        Ignored:    data/NAPA_clim_sst.Rda
        Ignored:    data/NAPA_clim_taum.Rda
        Ignored:    data/NAPA_clim_vars.Rda
        Ignored:    data/NAPA_clim_vecs.Rda
        Ignored:    data/node_mean_all_anom.Rda
        Ignored:    data/packet_all_anom.Rda
        Ignored:    data/som_all_anom.Rda
        Ignored:    data/synoptic_states.Rda
        Ignored:    data/synoptic_vec_states.Rda
    
    Untracked files:
        Untracked:  output/som_plot_uoce_anom.pdf
        Untracked:  output/som_plot_voce_anom.pdf
    
    Unstaged changes:
        Modified:   analysis/MHWNWA.bib
        Modified:   analysis/_site.yml
        Modified:   code/workflow.R
        Modified:   output/som_plot_mldr10_1_anom.pdf
        Modified:   output/som_plot_qt_anom.pdf
        Modified:   output/som_plot_sst_anom.pdf
        Modified:   output/som_plot_taum_anom.pdf
    
    
    Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
Expand here to see past versions:
    File Version Author Date Message
    html 028d3cc robwschlegel 2019-06-10 Build site.
    Rmd c6b3c7b robwschlegel 2019-06-10 Re-publish entire site.

Introduction

In this final vignette we will go over the creation of the figures used in the publication for this research. These figures are largely adapted from the techniques seen in Oliver et al. (2018) and Schlegel et al. (2017).

# Insatll from GitHub
# .libPaths(c("~/R-packages", .libPaths()))
# devtools::install_github("fabrice-rossi/yasomi")

# Packages used in this vignette
library(jsonlite, lib.loc = "../R-packages/")
library(tidyverse) # Base suite of functions
library(ncdf4) # For opening and working with NetCDF files
library(lubridate) # For convenient date manipulation
# library(scales) # For scaling data before running SOM
library(yasomi, lib.loc = "../R-packages/") # The SOM package of choice due to PCI compliance
library(data.table) # For working with massive dataframes

# Set number of cores
doMC::registerDoMC(cores = 50)

# Disable scientific notation for numeric values
  # I just find it annoying
options(scipen = 999)

# Set number of cores
doMC::registerDoMC(cores = 50)

# Disable scientific notation for numeric values
  # I just find it annoying
options(scipen = 999)

# Individual regions
NWA_coords <- readRDS("data/NWA_coords_cabot.Rda")

# The NAPA variables
NAPA_vars <- readRDS("data/NAPA_vars.Rda")

# Corners of the study area
NWA_corners <- readRDS("data/NWA_corners.Rda")

# Create smaller corners to use less RAM
  # This also better matches the previous South African work
  # The Tasmania work had corners of roughly 2 degrees greater than the study area
NWA_corners_sub <- c(NWA_corners[1]+8, NWA_corners[2]-8, NWA_corners[3]+8, NWA_corners[4]-8)

# The base map
map_base <- ggplot2::fortify(maps::map(fill = TRUE, col = "grey80", plot = FALSE)) %>%
  dplyr::rename(lon = long) %>%
  mutate(group = ifelse(lon > 180, group+9999, group),
         lon = ifelse(lon > 180, lon-360, lon)) %>% 
  select(-region, -subregion)
###########################################################################
### "3.Figures.R"
## This script creates the figures for the paper and supplemental
# 1. Load all libraries and functions used in this script
# 2. Create synoptic figure for each event
# 3. Create synoptic figure showing SOM nodes
# 4. Create lolliplots for the SOM nodes
# 5. Create dendrogram for HCA results
# 6. Create ordiplot for MDS results
# 7. Create map of study area
#############################################################################


# 1. Load all libraries and functions used in this script -----------------
source("func/synoptic.func.R")
source("func/som.func.R")
source("func/scale.bar.func.R")
library(vegan)
library(ggdendro)
library(broom)


# 2. Create synoptic figure for each event  -------------------------------

# Load SACTN data
load("~/data/SACTN/AHW/SACTN_clims.Rdata")
load("data/SACTN/SACTN_events.Rdata")
load("setupParams/SACTN_site_list.Rdata")

# The files for loading
event_idx <- data.frame(event = dir("data/SOM", full.names = TRUE),
                        x = length(dir("data/SOM")))

# Create a synoptic atlas figure for each MHW
system.time(plyr::ddply(event_idx, c("event"), synoptic.fig, .progress = "text")) # 539 seconds


# 3. Create synoptic figure showing SOM nodes -----------------------------

load("data/node_means.Rdata")
load("data/node_all_anom.Rdata")
all.panels(node_means, node_all_anom)


# 4. Create lolliplots for the SOM nodes ----------------------------------

# Load data for figure
load("data/SACTN/SACTN_events.Rdata")
load("data/node_all_anom.Rdata")

# Merge into one dataframe
node_all <- merge(node_all_anom, SACTN_events, by = c("event", "site", "season", "event_no"))

# Calculate mean and median per node for plotting
node_h_lines <- node_all %>% 
  group_by(node) %>% 
  summarise(mean_int_cum = mean(int_cum, na.rm = T),
            median_int_cum = median(int_cum, na.rm = T))

# Create the figure
ggplot(data = node_all, aes(x = date_start, y = int_cum)) +
  geom_lolli() +
  geom_point(aes(colour = season)) +
  geom_label(aes(x = as.Date("2005-01-01"), y = 580, label = paste0("n = ", count,"/",length(node))), 
             size = 3, label.padding = unit(0.5, "lines")) +
  geom_hline(data = node_h_lines, aes(yintercept = mean_int_cum), linetype = "dashed") +
  geom_hline(data = node_h_lines, aes(yintercept = median_int_cum), linetype = "dotted") +
  facet_wrap(~node) +
  labs(x = "", y = "Cummulative intensity (°C·days)", colour = "Season") +
  theme_grey() +
  # scale_y_continuous(expand = c(0, 100)) +
  theme(strip.background = element_rect(fill = NA),
        panel.border = element_rect(fill = NA, colour = "black", size = 1),
        axis.text = element_text(size = 12, colour = "black"),
        axis.ticks = element_line(colour = "black"))
ggsave("graph/SOM_lolli.pdf", height = 9, width = 10)


# 5. Create dendrogram for HCA results ------------------------------------

# Load data
load("data/all_anom_hclust.Rdata")
load("data/all_anom_env.Rdata")

# Prep dendrogram
dhc <- as.dendrogram(all_anom_hclust)
ddata <- dendro_data(dhc, type = "rectangle")
label_data <- label(ddata)
label_data$season <- all_anom_env$season[as.numeric(as.character(label_data$label))]
label_data$type <- all_anom_env$type[as.numeric(as.character(label_data$label))]

# Plot the dendrogram
ggplot(segment(ddata)) + 
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_point(data = label_data, aes(x = x, y = y, shape = type, colour = season)) +
  scale_shape_manual(values = c(1, 16)) +
  theme_grey() +
  theme(panel.border = element_rect(fill = NA, colour = "black", size = 1),
        axis.text = element_text(size = 12, colour = "black"),
        axis.ticks = element_line(colour = "black"))
ggsave("graph/HCA.pdf", height = 9, width = 9)


# 6. Create ordiplot for MDS results --------------------------------------

# Load data
load("data/all_anom_MDS.Rdata")
load("data/all_anom_env.Rdata")
load("data/SACTN/SACTN_events.Rdata")
SACTN_events$type <- NULL
# load("data/node_means.Rdata")
# load("data/node_all_anom.Rdata")

## Fit environmental variables
ord_fit <- envfit(all_anom_MDS ~ season + type, data = all_anom_env)
# ord_fit
ord_fit_df <- as.data.frame(ord_fit$factors$centroids)
ord_fit_df$factors <- c("autumn", "spring", "summer", "winter", "clim", "MHW")

# Merge event values
# SACTN_events_sub <- SACTN_events[,c(1,2)]
all_anom_env <- left_join(all_anom_env, SACTN_events, by = c("event", "season"))

# Create MDS dataframe
# mds_df <- data.frame(all_anom_MDS$points, type = all_anom_env$type, 
#                      event = all_anom_env$event, season = all_anom_env$season)
mds_df <- data.frame(all_anom_MDS$points, all_anom_env)
mds_df$duration[is.na(mds_df$duration)] <- mean(mds_df$duration, na.rm = T)

# Plot the fits
ggplot(data = mds_df, aes(x = MDS1, y = MDS2)) +
  geom_point(aes(colour = season, shape = type, size = duration)) +
  geom_segment(data = ord_fit_df, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
               arrow = arrow(angle = 40, length = unit(0.2, "cm"), type = "open"), 
               alpha = 1, colour = "black", size = 0.5)  +
  geom_text(data = ord_fit_df, aes(label = factors, x = NMDS1, y = NMDS2), size = 8) +
  scale_shape_manual(name = "State", values = c(19, 15), labels = c("clim", "MHW")) +
  scale_colour_discrete(name = "Season") +
  scale_size_continuous(name = "Duration\n(days)", breaks = c(20, 70, 120, 170, 220)) +
  guides(colour = guide_legend(order = 1),
         shape = guide_legend(order = 2),
         size = guide_legend(override.aes = list(shape = 15), order = 3)) +
  # labs(size = "Duration") +
  theme_grey() +
  theme(strip.background = element_rect(fill = NA),
        panel.border = element_rect(fill = NA, colour = "black", size = 1),
        axis.text = element_text(size = 12, colour = "black"),
        axis.ticks = element_line(colour = "black"))
ggsave("graph/MDS.pdf", height = 9, width = 12)

# Linear model: duration vs. ordination from centre point
lm_all_results <- mds_df %>%
  na.omit() %>%
  mutate(MDS2 = abs(MDS2)) %>% 
  select(-MDS1, -event, -season, -type, -coast, -site) %>%
  gather(key = group, 
         value = measurement,
         -MDS2) %>%
  group_by(group) %>% 
  nest() %>%
  mutate(model = purrr::map(data, ~lm(measurement ~ MDS2, data = .))) %>% 
  unnest(model %>% purrr::map(glance)) %>% 
  select(-data, -model)
lm_all_results[lm_all_results$adj.r.squared == max(lm_all_results$adj.r.squared, na.rm = T),]

# Visualization of the best linear model
ggplot(data = mds_df) +
  geom_smooth(aes(x = duration, y = abs(MDS2)), method = "lm", se = F) +
  theme_grey()


# 7. Create map of study area ---------------------------------------------

## Load data
# International borders
load("graph/africa_borders.Rdata")

# Hi-res bathy
load("graph/bathy.Rdata")

# Reanalysis data
load("data/ERA/ERA_temp_clim.Rdata")
load("data/ERA/ERA_u_clim.Rdata")
colnames(ERA_u_clim)[4] <- "val"
ERA_u_clim$variable <- "u"
load("data/ERA/ERA_v_clim.Rdata")
colnames(ERA_v_clim)[4] <- "val"
ERA_v_clim$variable <- "v"

# Remote data
load("data/OISST/OISST_temp_clim.Rdata")
load("data/AVISO/AVISO_u_clim.Rdata")
colnames(AVISO_u_clim)[4] <- "val"
AVISO_u_clim$variable <- "u"
load("data/AVISO/AVISO_v_clim.Rdata")
colnames(AVISO_v_clim)[4] <- "val"
AVISO_v_clim$variable <- "v"

# In situ time series locations
load("setupParams/SACTN_site_list.Rdata")
SACTN_site_list$order <- 1:nrow(SACTN_site_list)

# Create annual mean air-sea state
# sea_temp <- filter(OISST_temp_clim, date == "01-15") %>% 
sea_temp <- data.table::data.table(OISST_temp_clim)
sea_temp <- sea_temp[, .(temp = mean(temp, na.rm = TRUE)),
           by = .(x, y)] %>% 
  rename(lon = x, lat = y)
currents <- data.table::data.table(rbind(AVISO_u_clim, AVISO_v_clim))
currents <- currents[, .(val = mean(val, na.rm = TRUE)),
                     by = .(x, y, variable)] %>% 
  spread(key = variable, value = val) %>% 
  rename(lon = x, lat = y) %>% 
  mutate(arrow_size = ((abs(u*v)/ max(abs(u*v)))+0.2)/6)
air_temp <- data.table::data.table(ERA_temp_clim)
air_temp <- air_temp[, .(temp = mean(temp, na.rm = TRUE)),
                     by = .(x, y)] %>% 
  rename(lon = x, lat = y)
winds <- data.table::data.table(rbind(ERA_u_clim, ERA_v_clim))
winds <- winds[, .(val = mean(val, na.rm = TRUE)),
                     by = .(x, y, variable)] %>% 
  spread(key = variable, value = val) %>% 
  rename(lon = x, lat = y) %>% 
  mutate(arrow_size = ((abs(u*v)/ max(abs(u*v)))+0.3)/6)

# Reduce wind/ current vectors
lon_sub <- seq(10, 40, by = 1)
lat_sub <- seq(-40, -15, by = 1)
# currents <- currents[(currents$lon %in% lon_sub & currents$lat %in% lat_sub),]
winds <- winds[(winds$lon %in% lon_sub & winds$lat %in% lat_sub),]

# Establish the vector scalar for the currents
current_uv_scalar <- 2

# Establish the vector scalar for the wind
wind_uv_scalar <- 0.5

# Wind feature vector coordinates
cyc_atlantic <- data.frame(x = c(14.0, 16.1, 16.0), y = c(-36.0, -34.4, -32.1), xend = c(16.0, 16.1, 14.0), yend = c(-34.5, -32.2, -30.6))
cyc_indian <- data.frame(x = c(36.0, 33.9, 34.0), y = c(-31.5, -33.1, -35.4), xend = c(34.0, 33.9, 36.0), yend = c(-33.0, -35.3, -36.9))
westerlies <- data.frame(x = c(18.0, 21.1, 24.2), y = c(-38.0, -37.8, -37.8), xend = c(21.0, 24.1, 27.2), yend = c(-37.8, -37.8, -38.0))

# The top figure (sea)
fig_1_top <- ggplot(data = southern_africa_coast, aes(x = lon, y = lat)) +
  # The ocean temperature
  geom_raster(data = sea_temp, aes(fill = temp)) +
  # The bathymetry
  stat_contour(data = bathy[bathy$depth < -100 & bathy$depth > -300,], 
               aes(x = lon, y = lat, z = depth), alpha = 0.5,
               colour = "ivory", size = 0.5, binwidth = 200, na.rm = TRUE, show.legend = FALSE) +
  # The current vectors
  geom_segment(data = currents, aes(xend = lon + u * current_uv_scalar, yend = lat + v * current_uv_scalar),
               arrow = arrow(angle = 40, length = unit(currents$arrow_size, "cm"), type = "open"),
                             linejoin = "mitre", size = 0.4) +
  # The land mass
  geom_polygon(aes(group = group), fill = "grey70", colour = "black", size = 0.5, show.legend = FALSE) +
  geom_path(data = africa_borders, aes(group = group)) +
  # The legend for the vector length
  geom_label(aes(x = 37.0, y = -38.0, label = "1.0 m/s\n"), size = 3, label.padding = unit(0.5, "lines")) +
  geom_segment(aes(x = 36.0, y = -38.3, xend = 38.0, yend = -38.3), linejoin = "mitre",
               arrow = arrow(angle = 40, length = unit(0.2, "cm"), type = "open")) +
  # The in situ sites
  geom_point(data = SACTN_site_list, shape = 19,  size = 2.8, colour = "ivory") +
  geom_text(data = SACTN_site_list[-c(3,4,7:9,18,21,23:24),], aes(label = order), size = 1.9, colour = "red") +
  # Oceans
  annotate("text", label = "INDIAN\nOCEAN", x = 37.00, y = -34.0, size = 4.0, angle = 0, colour = "ivory") +
  
  annotate("text", label = "ATLANTIC\nOCEAN", x = 13.10, y = -34.0, size = 4.0, angle = 0, colour = "ivory") +
  # Benguela
  geom_segment(aes(x = 17.2, y = -32.6, xend = 15.2, yend = -29.5),
               arrow = arrow(length = unit(0.3, "cm")), size = 0.5, colour = "ivory") +
  annotate("text", label = "Benguela", x = 16.0, y = -31.8, size = 3.5, angle = 298, colour = "ivory") +
  # Agulhas
  geom_segment(aes(x = 33, y = -29.5, xend = 29.8, yend = -33.0),
               arrow = arrow(length = unit(0.3, "cm")), size = 0.5, colour = "ivory") +
  annotate("text", label = "Agulhas", x = 31.7, y = -31.7, size = 3.5, angle = 53, colour = "ivory") +
  # Agulhas Bank
  annotate("text", label = "Agulhas\nBank", x = 22.5, y = -35.5, size = 3.0, angle = 0, colour = "ivory") +
  # Cape Peninsula
  annotate("text", label = "Cape\nPeninsula", x = 17.2, y = -35, size = 3.0, angle = 0, colour = "ivory") +
  # Improve on the x and y axis labels
  scale_x_continuous(breaks = seq(15, 35, 5),
                     labels = scales::unit_format("°E", sep = ""),
                     position = "top") +
  scale_y_continuous(breaks = seq(-35, -30, 5),
                     labels = c("35°S", "30°S")) +
  labs(x = NULL, y = NULL) +
  # Slightly shrink the plotting area
  coord_cartesian(xlim = c(10.5, 39.5), ylim = c(-39.5, -25.5), expand = F) +
  # Use viridis colour scheme
  scale_fill_viridis(name = "Temp.\n(°C)", option = "D", breaks = c(24, 20, 16)) +
  # Adjust the theme
  theme_bw() +
  theme(panel.border = element_rect(fill = NA, colour = "black", size = 1),
        axis.text = element_text(size = 12, colour = "black"),
        axis.ticks = element_line(colour = "black"))
# fig_1_top

# False Bay inset
fb_inset <- ggplot(data = sa_shore, aes(x = lon, y = lat)) +
  # The land mass
  geom_polygon(aes(group = PID),
               fill = "grey70", colour = NA, size = 0.5, show.legend = FALSE) +
  # The in situ sites
  geom_point(data = SACTN_site_list, shape = 1,  size = 3, colour = "black") +
  geom_text(data = SACTN_site_list[-6,], aes(label = order), size = 2.0, colour = "red") +
  # Text label
  geom_text(aes(x = 18.65, y = -34.25, label = "False\nBay"), size = 2.7) +
  # Control the x and y axes
  coord_cartesian(xlim = c(18.2, 19), ylim = c(-34.5, -33.8), expand = F) +
  scale_x_continuous(breaks = c(18.5), label = "18.5°E") +
  scale_y_continuous(breaks = c(-34.1), label = "34.1°S") +
  labs(x = NULL, y = NULL) +
  # Change the theme for cleaner over-plotting
  theme_bw() +
  theme(plot.background = element_blank(),
        axis.text = element_text(colour = "ivory"),
        axis.text.y = element_text(angle = 90, hjust = 0.5),
        axis.ticks = element_line(colour = "ivory"),
        panel.border = element_rect(colour = "ivory"),
        panel.grid = element_blank())
# fb_inset

# The bottom figure (air)
fig_1_bottom <- ggplot(data = southern_africa_coast, aes(x = lon, y = lat)) +
  # The ocean temperature
  geom_raster(data = air_temp, aes(fill = temp)) +
  # The land mass
  geom_polygon(aes(group = group), fill = NA, colour = "black", size = 0.5, show.legend = FALSE) +
  geom_path(data = africa_borders, aes(group = group)) +
  # The current vectors
  geom_segment(data = winds, aes(xend = lon + u * wind_uv_scalar, yend = lat + v * wind_uv_scalar),
               arrow = arrow(angle = 40, length = unit(winds$arrow_size, "cm"), type = "open"),
               linejoin = "mitre", size = 0.4) +
  # The legend for the vector length
  geom_label(aes(x = 37.0, y = -38.0, label = "4.0 m/s\n"), size = 3, label.padding = unit(0.5, "lines")) +
  geom_segment(aes(x = 36.0, y = -38.3, xend = 38.0, yend = -38.3), linejoin = "mitre",
               arrow = arrow(angle = 40, length = unit(0.2, "cm"), type = "open")) +
  # The coastal sections
  geom_spoke(aes(x = 18.46520, y = -34.31050, angle = 180, radius = -2), linetype = "dotted", colour = "ivory") +
  geom_spoke(aes(x = 18.46520, y = -34.31050, angle = 180, radius = 2), linetype = "dotted", colour = "ivory") +
  geom_spoke(aes(x = 27.48889, y = -33.28611, angle = 40, radius = -2), linetype = "dotted", colour = "ivory") +
  geom_spoke(aes(x = 27.48889, y = -33.28611, angle = 40, radius = 2), linetype = "dotted", colour = "ivory") +
  annotate("text", label = "West\nCoast", x = 19.5, y = -31.2, size = 3.0, angle = 0, colour = "ivory") +
  annotate("text", label = "South\nCoast", x = 23, y = -33.0, size = 3.0, angle = 0, colour = "ivory") +
  annotate("text", label = "East\nCoast", x = 28, y = -31, size = 3.0, angle = 0, colour = "ivory") +
  # South Atlantic Anticyclone
  annotate("text", label = "SOUTH\nATLANTIC\nANTICYCLONE", x = 13.5, y = -33.5, size = 3.0, angle = 0, colour = "ivory") +
  geom_curve(data = cyc_atlantic, aes(x = x, y = y, xend = xend, yend = yend), curvature = 0.2, colour = "ivory",
             arrow = arrow(angle = 40, type = "open", length = unit(0.25,"cm"))) +
  # South Indian Anticyclone
  annotate("text", label = "SOUTH\nINDIAN\nANTICYCLONE", x = 36.5, y = -34.0, size = 3.0, angle = 0, colour = "ivory") +
  geom_curve(data = cyc_indian, aes(x = x, y = y, xend = xend, yend = yend), curvature = 0.2, colour = "ivory",
             arrow = arrow(angle = 40, type = "open", length = unit(0.25,"cm"))) +
  # Westerlies
  annotate("text", label = "WESTERLIES", x = 22.5, y = -37.0, size = 3.0, angle = 0, colour = "ivory") +
  geom_curve(data = westerlies, aes(x = x, y = y, xend = xend, yend = yend), colour = "ivory",
             arrow = arrow(angle = 40, type = "open", length = unit(0.25,"cm")), curvature = -0.01) +
  # Improve on the x and y axis labels
  scale_x_continuous(breaks = seq(15, 35, 5),
                     labels = scales::unit_format("°E", sep = "")) +
  scale_y_continuous(breaks = seq(-35, -30, 5),
                     labels = c("35°S", "30°S")) +
  labs(x = NULL, y = NULL) +
  # Scale bar
  scaleBar(lon = 22.0, lat = -29.5, distanceLon = 200, distanceLat = 50, distanceLegend = 90, dist.unit = "km",
           arrow.length = 100, arrow.distance = 130, arrow.North.size = 3, 
           legend.colour = "ivory", arrow.colour = "ivory", N.colour = "ivory") +
  # Slightly shrink the plotting area
  coord_cartesian(xlim = c(10.5, 39.5), ylim = c(-39.5, -25.5), expand = F) +
  # Use viridis colour scheme
  scale_fill_viridis(name = "Temp.\n(°C)", option = "A", breaks = c(24, 20, 16)) +
  # Adjust the theme
  theme_bw() +
  theme(panel.border = element_rect(fill = NA, colour = "black", size = 1),
        axis.text = element_text(size = 12, colour = "black"),
        axis.ticks = element_line(colour = "black"))
# fig_1_bottom

# Convert the figures to grobs
fig_1_top_grob <- ggplotGrob(fig_1_top)
fb_inset_grob <- ggplotGrob(fb_inset)
fig_1_bottom_grob <- ggplotGrob(fig_1_bottom)

# Stick them together
fig_1 <- ggplot() +
  # First set the x and y axis values so we know what the ranges are
  # in order to make it easier to place our facets
  coord_equal(xlim = c(1, 10), ylim = c(1, 10), expand = F) +
  # Then we place our facetsover one another using the coordinates we created
  annotation_custom(fig_1_top_grob,
                    xmin = 1, xmax = 10, ymin = 5.5, ymax = 10) +
  annotation_custom(fb_inset_grob,
                    xmin = 3.5, xmax = 5.5, ymin = 7.2, ymax = 8.8) +
  annotation_custom(fig_1_bottom_grob,
                    xmin = 1, xmax = 10, ymin = 1, ymax = 5.5)
# save
ggsave(plot = fig_1, filename = "graph/fig_1.pdf", height = 8, width = 8)

Oliver, E. C., Lago, V., Hobday, A. J., Holbrook, N. J., Ling, S. D., and Mundy, C. N. (2018). Marine heatwaves off eastern tasmania: Trends, interannual variability, and predictability. Progress in oceanography 161, 116–130.

Schlegel, R. W., Oliver, E. C., Perkins-Kirkpatrick, S., Kruger, A., and Smit, A. J. (2017). Predominant atmospheric and oceanic patterns during coastal marine heatwaves. Frontiers in Marine Science 4, 323.

Session information

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.18.so

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
 [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2    data.table_1.11.6 yasomi_0.3       
 [4] proxy_0.4-22      e1071_1.7-0       lubridate_1.7.4  
 [7] ncdf4_1.16        forcats_0.3.0     stringr_1.3.1    
[10] dplyr_0.7.6       purrr_0.2.5       readr_1.1.1      
[13] tidyr_0.8.1       tibble_1.4.2      ggplot2_3.0.0    
[16] tidyverse_1.2.1   jsonlite_1.6     

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4  haven_1.1.2       lattice_0.20-35  
 [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.2.0       
 [7] rlang_0.2.2       R.oo_1.22.0       pillar_1.3.0     
[10] glue_1.3.0        withr_2.1.2       R.utils_2.7.0    
[13] doMC_1.3.5        modelr_0.1.2      readxl_1.1.0     
[16] foreach_1.4.4     bindr_0.1.1       plyr_1.8.4       
[19] munsell_0.5.0     gtable_0.2.0      workflowr_1.1.1  
[22] cellranger_1.1.0  rvest_0.3.2       R.methodsS3_1.7.1
[25] codetools_0.2-15  evaluate_0.11     knitr_1.20       
[28] parallel_3.6.0    class_7.3-14      broom_0.5.0      
[31] Rcpp_0.12.18      backports_1.1.2   scales_1.0.0     
[34] hms_0.4.2         digest_0.6.16     stringi_1.2.4    
[37] grid_3.6.0        rprojroot_1.3-2   cli_1.0.0        
[40] tools_3.6.0       maps_3.3.0        magrittr_1.5     
[43] lazyeval_0.2.1    crayon_1.3.4      whisker_0.3-2    
[46] pkgconfig_2.0.2   xml2_1.2.0        iterators_1.0.10 
[49] assertthat_0.2.0  rmarkdown_1.10    httr_1.3.1       
[52] rstudioapi_0.7    R6_2.2.2          nlme_3.1-137     
[55] git2r_0.23.0      compiler_3.6.0   

This reproducible R Markdown analysis was created with workflowr 1.1.1