README
This document contains the R workflow for processing and analysing the raw data, and generating figures for the manuscript titled “Global exposure risk of frogs to increasing environmental dryness”.

GitHub Repository The Rmarkdown file can be downloaded from the Code drop down menu (top right).

Abbreviations

  • AI: Aridity Index
  • \(e_a\): actual vapour pressure
  • \(e_s\): saturation vapour pressure
  • EWL: Evaporative water loss
  • IUCN: International Union for Conservation of Nature
  • PDSI: Palmer Drought Severity Index
  • PRISMA: Preferred Reporting Items for Systematic Reviews and Meta-Analyses
  • RH: Relative humidity
  • \(r_i\): Relative skin resistance or resistance to water loss
  • SA: surface-area
  • SPP: Shared Socioeconomic Pathway
  • VPD: Vapour pressure deficit
  • WoS: Web of Science
  • WU: water uptake

Spatial risk

Descriptors

The following series of code is to estimate the risk of anuran amphibians to increasing aridity (Aridity Index; AI) and extreme drought (Palmer Drought Severity Index; PDSI). The AI was categorised to five categories (Humid, Dry sub-humid, Semi-arid, Arid, and Hyper-arid) and the PDSI was categorised to seven categories (Extremely moist, Very moist, Moderate moist, Normal, Moderate drought, Severe drought, Extreme drought) based on descriptions from Budyko (1961) and Palmer (1965), respectively (Table S1).

Generalised ecological types or ‘ecotype’ of each anuran species were classified based on descriptions from Moen and Wiens (2017) focusing on adult behaviour and microhabitat preferences outside the breeding season, given that many anurans breed in water but are not adapted to live in water all year (Table S2).

Table S1

Table S1 Aridity index (Budyko, 1961) and the Palmer Drought Severity Index (Palmer, 1965) categorised.

Aridity index (AI) AI value Palmer Drought Severity Index (PDSI) PDSI Value
Humid ≥ 0.65 Extremely moist ≥ 4
Dry sub-humid 0.50 – ≤ 0.65 Very moist 3 – ≤ 4
Semi-arid 0.20 – ≤ 0.50 Moderate moist 2 – ≤ 3
Arid 0.05 – ≤ 0.20 Normal -2 – ≤ 2
Hyper-arid < 0.05 Moderate drought -3 – ≤ -2
Severe drought -4 – ≤ -3
Extreme drought < -4

Table S2

Table S2 Microhabitat preference or ecotype (Moen and Wiens, 2017).

Ecotype Description
Aquatic Almost always in water.
Arboreal Ability to climb vegetation and spends a substantial amount of time above ground (e.g. forests, on top of vegetations in wetlands).
Semi-aquatic Stay near water bodies, usually permanent water (e.g. ground-level wetlands). Often observed in water and nearby vegetation.
Ground-dwelling Terrestrial only, not always near water bodies, however always in a moist environment (e.g. forest leaf litter, moist soils, scrublands, grassland). Sometimes above ground and sometimes semi-fossorial. Can be classified as broad.
Fossorial The non-breeding season is spent underground in burrows. Capable of burrowing and aestivating.
Stream-dwelling Restricted to near fast-flowing streams usually on rocks or vegetation near streams.

Species richness

Anuran (frogs and toads) distribution shape files were obtained from the International Union for Conservation of Nature’s Red List of Threatened Species (IUCN Red List; extracted on 11/01/2021). Species richness or species assemblages was defined as the sum of species in each grid cell (0.5°), based on the geographic range, and was calculated using the rasterizeIUCN and calcSR function from the from the rasterSp package. Ecotype-specific species richness was also calculated.

filedir <- "rasterSp/"

# Clean raw data
frog_dat <- read.csv(file.path(data_path, "anuran_species_list.csv")) %>%
    dplyr::select(genus:strategy, SVL_cm:binomial_tree_phylo) %>%
    dplyr::mutate(lnSVL = log(SVL_cm), lnMass = log(mass_g), lnArea = log(shape_area),
        IUCN = factor(IUCN, levels = c("Least Concern", "Near Threatened", "Vulnerable",
            "Endangered", "Critically Endangered", "Extinct")), order_name = factor(order_name),
        family_name = factor(family_name), ecotype = factor(ecotype)) %>%
    dplyr::filter(order_name == "ANURA" & ecotype != "" & IUCN != "" & IUCN != "Extinct") %>%
    droplevels()

# Convert shape files into rasters and save to file once (downloaded shape
# files 08/01/2020) only extracted distribution of native range anuran_rast <-
# rasterizeIUCN(dsn = paste0(filedir, 'ANURA/ANURA.shp'), resolution = 0.5,
# seasonal = c(1, 2), origin = 1, presence = c(1,2), save = TRUE, path =
# paste0(rasterSp_path))

# Calculate anuran richness and by ecotype The calcSR function uses a stepwise
# procedure to calculate the sum of species for each grid cell.
anuran_sr <- rasterSp::calcSR(species_names = frog_dat$binomial_IUCN, path = paste0(rasterSp_path))
aquatic_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Aquatic",
    "binomial_IUCN"], path = paste0(rasterSp_path))
arboreal_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Arboreal",
    "binomial_IUCN"], path = paste0(rasterSp_path))
fossorial_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Fossorial",
    "binomial_IUCN"], path = paste0(rasterSp_path))
ground_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Ground-dwelling",
    "binomial_IUCN"], path = paste0(rasterSp_path))
semi_aq_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Semi-aquatic",
    "binomial_IUCN"], path = paste0(rasterSp_path))
stream_sr <- rasterSp::calcSR(species_names = frog_dat[frog_dat$ecotype == "Stream-dwelling",
    "binomial_IUCN"], path = paste0(rasterSp_path))

# Convert raster to matrix then to data frame
anuran_sr_df <- raster::as.data.frame(raster::rasterToPoints(anuran_sr)) %>%
    dplyr::rename(species_n = layer)
aquatic_sr_df <- raster::as.data.frame(raster::rasterToPoints(aquatic_sr)) %>%
    dplyr::rename(aquatic_n = layer)
arboreal_sr_df <- raster::as.data.frame(raster::rasterToPoints(arboreal_sr)) %>%
    dplyr::rename(arboreal_n = layer)
fossorial_sr_df <- raster::as.data.frame(raster::rasterToPoints(fossorial_sr)) %>%
    dplyr::rename(fossorial_n = layer)
ground_sr_df <- raster::as.data.frame(raster::rasterToPoints(ground_sr)) %>%
    dplyr::rename(ground_n = layer)
semi_aq_sr_df <- raster::as.data.frame(raster::rasterToPoints(semi_aq_sr)) %>%
    dplyr::rename(semi_aq_n = layer)
stream_sr_df <- raster::as.data.frame(raster::rasterToPoints(stream_sr)) %>%
    dplyr::rename(stream_n = layer)

There were 5636 anuran species from the IUCN Red List available for the AI and PDSI analysis.

Calculate AI and PDSI

High resolution (~4 km2) global dataset on precipitation (mm/month) and potential evapotranspiration (mm/month) were obtained from Abatzoglou et al. (2018) under 1) the current climate (1970–2000), 2) an intermediate greenhouse gas emission scenario of +2°C (Shared Socioeconomic Pathway 2–4.5; SPP2–4.5), and 3) a high greenhouse gas emission or “business-as-usual” scenario of +4°C (SSP5–8.5) by 2080–2099 (Fig. S1a–f).

We obtained a self-calibrated PDSI with Penman–Monteith potential evapotranspiration representing the current climate (1970–2000) and an intermediate and high emission scenario by 2080–2099 (SSP2–4.5 and SSP5-8.5) from Zhao and Dai (2022). The SSP2–4.5 and SSP5–8.5 scenarios were based on the average of 25 CMIP6 models of precipitation, evapotranspiration, soil moisture, and runoff (Eyring et al., 2016), where the mean annual surface temperature is expected to increase by 2.7°C (2.1–3.5°C range) and 4.4°C (3.3–5.7°C range), respectively by 2080–2100 (IPCC, 2021).

# Import the downloaded files
ppt_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate19812010_ppt.nc")),
    anuran_sr)  # Precipitation year 1981-2010 - Pr (mm)
pet_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate19812010_pet.nc")),
    anuran_sr)  # Evapotranspiration 1981-2010 - ET0 (mm)

# +2C and +4C future scenarios
ppt_2C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate2C_ppt.nc")),
    anuran_sr)
pet_2C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate2C_pet.nc")),
    anuran_sr)
ppt_4C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate4C_ppt.nc")),
    anuran_sr)
pet_4C_rast <- raster::projectRaster(raster::stack(x = file.path(spatial_path, "TerraClimate4C_pet.nc")),
    anuran_sr)

# Obtain mean values from 1981-2010
ppt_mean_rast <- raster::calc(ppt_rast, fun = mean, na.rm = TRUE)
pet_mean_rast <- raster::calc(pet_rast, fun = mean, na.rm = TRUE)
ppt_2C_mean_rast <- raster::calc(ppt_2C_rast, fun = mean, na.rm = TRUE)
pet_2C_mean_rast <- raster::calc(pet_2C_rast, fun = mean, na.rm = TRUE)
ppt_4C_mean_rast <- raster::calc(ppt_4C_rast, fun = mean, na.rm = TRUE)
pet_4C_mean_rast <- raster::calc(pet_4C_rast, fun = mean, na.rm = TRUE)

# Obtain mean values from 10-year monthly PDSI
PDSI_ssp245_rast <- raster::stack(x = file.path(spatial_path, "pdsisc.monthly.1900-2100.r2.5x2.5.EnsAvg25Models.TP2.ipe-2.ssp245.nc"))
PDSI_ssp585_rast <- raster::stack(x = file.path(spatial_path, "pdsisc.monthly.1900-2100.r2.5x2.5.EnsAvg25Models.TP2.ipe-2.ssp585.nc"))

PDSI_cur_rast <- raster::subset(PDSI_ssp245_rast, grep("X197.*|X198.*|X199.*", names(PDSI_ssp245_rast),
    value = T))
PDSI_2C_rast <- raster::subset(PDSI_ssp245_rast, grep("X207.*|X208.*|X209.*", names(PDSI_ssp245_rast),
    value = T))
PDSI_4C_rast <- raster::subset(PDSI_ssp585_rast, grep("X207.*|X208.*|X209.*", names(PDSI_ssp585_rast),
    value = T))

PDSI_cur_mean_rast <- raster::calc(PDSI_cur_rast, fun = mean, na.rm = TRUE)
PDSI_2C_mean_rast <- raster::calc(PDSI_2C_rast, fun = mean, na.rm = TRUE)
PDSI_4C_mean_rast <- raster::calc(PDSI_4C_rast, fun = mean, na.rm = TRUE)

# Calculate aridity index [Precipitation (ppt) / Evapotranspiration (pet)]
ai_rast <- raster::overlay(x = ppt_mean_rast, y = pet_mean_rast, fun = function(x,
    y) {
    return(x/y)
})
ai_2C_rast <- raster::overlay(x = ppt_2C_mean_rast, y = pet_2C_mean_rast, fun = function(x,
    y) {
    return(x/y)
})
ai_4C_rast <- raster::overlay(x = ppt_4C_mean_rast, y = pet_4C_mean_rast, fun = function(x,
    y) {
    return(x/y)
})

# Calculate change in PDSI Reproject intensity raster to match anuran_sr
PDSI_cur_mean_rast <- raster::projectRaster(PDSI_cur_mean_rast, anuran_sr)
PDSI_2C_mean_rast <- raster::projectRaster(PDSI_2C_mean_rast, anuran_sr)
PDSI_4C_mean_rast <- raster::projectRaster(PDSI_4C_mean_rast, anuran_sr)

PDSI_2C_diff_rast <- raster::overlay(x = PDSI_cur_mean_rast, y = PDSI_2C_mean_rast,
    fun = function(x, y) {
        return(y - x)
    })
PDSI_4C_diff_rast <- raster::overlay(x = PDSI_cur_mean_rast, y = PDSI_4C_mean_rast,
    fun = function(x, y) {
        return(y - x)
    })

# Convert raster to matrix then to data frame
ai_df <- raster::as.data.frame(raster::rasterToPoints(ai_rast))
ai_2C_df <- raster::as.data.frame(raster::rasterToPoints(ai_2C_rast))
ai_4C_df <- raster::as.data.frame(raster::rasterToPoints(ai_4C_rast))

PDSI_2C_diff_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_2C_diff_rast))
PDSI_4C_diff_df <- raster::as.data.frame(raster::rasterToPoints(PDSI_4C_diff_rast))

# Convert aridity index into category
ai_df <- ai_df %>%
    # Recode
dplyr::mutate(category = case_when(is.infinite(layer) ~ "Humid", layer >= 0.65 ~
    "Humid", layer >= 0.5 & layer < 0.65 ~ "Dry sub-humid", layer >= 0.2 & layer <
    0.5 ~ "Semi-arid", layer >= 0.05 & layer < 0.2 ~ "Arid", layer < 0.05 ~ "Hyper-arid")) %>%
    # Convert to ordered factor
dplyr::mutate(category = factor(category, levels = c("Hyper-arid", "Arid", "Semi-arid",
    "Dry sub-humid", "Humid"), ordered = TRUE))

ai_2C_df <- ai_2C_df %>%
    # Recode
dplyr::mutate(category = case_when(is.infinite(layer) ~ "Humid", layer >= 0.65 ~
    "Humid", layer >= 0.5 & layer < 0.65 ~ "Dry sub-humid", layer >= 0.2 & layer <
    0.5 ~ "Semi-arid", layer >= 0.05 & layer < 0.2 ~ "Arid", layer < 0.05 ~ "Hyper-arid")) %>%
    # Convert to ordered factor
dplyr::mutate(category = factor(category, levels = c("Hyper-arid", "Arid", "Semi-arid",
    "Dry sub-humid", "Humid"), ordered = TRUE))

ai_4C_df <- ai_4C_df %>%
    # Recode
dplyr::mutate(category = case_when(is.infinite(layer) ~ "Humid", layer >= 0.65 ~
    "Humid", layer >= 0.5 & layer < 0.65 ~ "Dry sub-humid", layer >= 0.2 & layer <
    0.5 ~ "Semi-arid", layer >= 0.05 & layer < 0.2 ~ "Arid", layer < 0.05 ~ "Hyper-arid")) %>%
    # Convert to ordered factor
dplyr::mutate(category = factor(category, levels = c("Hyper-arid", "Arid", "Semi-arid",
    "Dry sub-humid", "Humid"), ordered = TRUE))

# Convert PDSI to category
PDSI_2C_diff_df <- PDSI_2C_diff_df %>%
    # Recode
dplyr::mutate(change_2C = case_when(layer >= 4 ~ "4", layer >= 3 & layer < 4 ~ "3",
    layer >= 2 & layer < 3 ~ "2", layer >= 1 & layer < 2 ~ "1", layer >= -1 & layer <
        1 ~ "0", layer >= -2 & layer < -1 ~ "-1", layer >= -3 & layer < -2 ~ "-2",
    layer >= -4 & layer < -3 ~ "-3", layer < -4 ~ "-4")) %>%
    # Convert to ordered factor
dplyr::mutate(change_2C = factor(change_2C, levels = c("-4", "-3", "-2", "-1", "0",
    "1", "2", "3", "4"), ordered = TRUE))

PDSI_4C_diff_df <- PDSI_4C_diff_df %>%
    # Recode
dplyr::mutate(change_4C = case_when(layer >= 4 ~ "4", layer >= 3 & layer < 4 ~ "3",
    layer >= 2 & layer < 3 ~ "2", layer >= 1 & layer < 2 ~ "1", layer >= -1 & layer <
        1 ~ "0", layer >= -2 & layer < -1 ~ "-1", layer >= -3 & layer < -2 ~ "-2",
    layer >= -4 & layer < -3 ~ "-3", layer < -4 ~ "-4")) %>%
    # Convert to ordered factor
dplyr::mutate(change_4C = factor(change_4C, levels = c("-4", "-3", "-2", "-1", "0",
    "1", "2", "3", "4"), ordered = TRUE))
# Crop and mask
world <- rgdal::readOGR(file.path(spatial_path, "ne_50m_land/ne_50m_land.shp"))
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/Drought project/Spatial data/ne_50m_land/ne_50m_land.shp", layer: "ne_50m_land"
## with 1420 features
## It has 3 fields
## Integer64 fields read as strings:  scalerank
world_spdf <- sp::SpatialPolygonsDataFrame(world, world@data)

# Precipitation
ppt_curr_plot <- raster::as.data.frame(raster::rasterToPoints(ppt_mean_rast)) %>%
    ggplot() + geom_raster(aes(y = y, x = x, fill = layer)) + geom_polygon(data = world_spdf,
    aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
    scale_fill_continuous_sequential(palette = "YlGnBu", lim = c(0, 650), name = "Precipitation (mm)") +
    ggtitle("Current (1981-2010)") + scale_y_continuous(limits = c(-60, 90), expand = c(0,
    0)) + scale_x_continuous(limits = c(-180, 180), expand = c(0, 0)) + guides(fill = guide_colourbar(barwidth = 10,
    barheight = 0.3, label.position = "bottom")) + theme_void() + ylab(NULL) + xlab(NULL) +
    theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(),
        plot.title = element_text(size = 10), legend.title = element_text(size = 8),
        legend.position = "bottom") + coord_fixed(ratio = 1)

ppt_2C_plot <- raster::as.data.frame(raster::rasterToPoints(ppt_2C_mean_rast)) %>%
    ggplot() + geom_raster(aes(y = y, x = x, fill = layer)) + geom_polygon(data = world_spdf,
    aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
    scale_fill_continuous_sequential(palette = "YlGnBu", lim = c(0, 650)) + ggtitle("+2°C (2080-2100)") +
    scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180,
    180), expand = c(0, 0)) + theme_void() + ylab(NULL) + xlab(NULL) + theme(axis.title = element_blank(),
    axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(size = 10)) +
    coord_fixed(ratio = 1)

ppt_4C_plot <- raster::as.data.frame(raster::rasterToPoints(ppt_4C_mean_rast)) %>%
    ggplot() + geom_raster(aes(y = y, x = x, fill = layer)) + geom_polygon(data = world_spdf,
    aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
    scale_fill_continuous_sequential(palette = "YlGnBu", lim = c(0, 650)) + ggtitle("+4°C (2080-2100)") +
    scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180,
    180), expand = c(0, 0)) + theme_void() + ylab(NULL) + xlab(NULL) + theme(axis.title = element_blank(),
    axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(size = 10)) +
    coord_fixed(ratio = 1)

ppt_legend <- cowplot::get_legend(ppt_curr_plot)

ppt_prow <- cowplot::plot_grid(ppt_curr_plot + theme(legend.position = "none"), ppt_2C_plot +
    theme(legend.position = "none"), ppt_4C_plot + theme(legend.position = "none"),
    align = "v", ncol = 3, labels = c("a", "b", "c"))

ppt_plots <- cowplot::plot_grid(ppt_prow, ppt_legend, ncol = 1, rel_heights = c(1,
    0.1))

# Evapotranspiration
pet_curr_plot <- raster::as.data.frame(raster::rasterToPoints(pet_mean_rast)) %>%
    ggplot() + geom_raster(aes(y = y, x = x, fill = layer)) + geom_polygon(data = world_spdf,
    aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
    scale_fill_continuous_sequential(palette = "Rocket", lim = c(0, 300), name = "Evapotranspiration (mm)") +
    scale_y_continuous(limits = c(-60, 90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180,
    180), expand = c(0, 0)) + guides(fill = guide_colourbar(barwidth = 10, barheight = 0.3,
    label.position = "bottom")) + theme_void() + ylab(NULL) + xlab(NULL) + theme(axis.title = element_blank(),
    axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(size = 10),
    legend.title = element_text(size = 8), legend.position = "bottom") + coord_fixed(ratio = 1)

pet_2C_plot <- raster::as.data.frame(raster::rasterToPoints(pet_2C_mean_rast)) %>%
    ggplot() + geom_raster(aes(y = y, x = x, fill = layer)) + geom_polygon(data = world_spdf,
    aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
    scale_fill_continuous_sequential(palette = "Rocket", lim = c(0, 300)) + scale_y_continuous(limits = c(-60,
    90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180, 180), expand = c(0,
    0)) + theme_void() + ylab(NULL) + xlab(NULL) + theme(axis.title = element_blank(),
    axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(size = 10)) +
    coord_fixed(ratio = 1)

pet_4C_plot <- raster::as.data.frame(raster::rasterToPoints(pet_4C_mean_rast)) %>%
    ggplot() + geom_raster(aes(y = y, x = x, fill = layer)) + geom_polygon(data = world_spdf,
    aes(x = long, y = lat, group = group), colour = "#64686b", fill = NA, size = 0.1) +
    scale_fill_continuous_sequential(palette = "Rocket", lim = c(0, 300)) + scale_y_continuous(limits = c(-60,
    90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180, 180), expand = c(0,
    0)) + theme_void() + ylab(NULL) + xlab(NULL) + theme(axis.title = element_blank(),
    axis.text = element_blank(), axis.ticks = element_blank(), plot.title = element_text(size = 10)) +
    coord_fixed(ratio = 1)

pet_legend <- cowplot::get_legend(pet_curr_plot)

pet_prow <- cowplot::plot_grid(pet_curr_plot + theme(legend.position = "none"), pet_2C_plot +
    theme(legend.position = "none"), pet_4C_plot + theme(legend.position = "none"),
    align = "v", ncol = 3, labels = c("d", "e", "f"))

pet_plots <- cowplot::plot_grid(pet_prow, pet_legend, ncol = 1, rel_heights = c(1,
    0.1))

# Aridity Index
arid_col <- c("#8E063B", "#CB6D53", "#E99A2C", "#F5D579", "white")
ai_curr_plot <- ggplot() + geom_raster(data = ai_df, aes(y = y, x = x, fill = category)) +
    geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b",
        fill = NA, size = 0.1) + scale_fill_manual(values = arid_col, guide = guide_legend(reverse = TRUE),
    name = "Aridity Index") + ylab(NULL) + xlab(NULL) + scale_y_continuous(limits = c(-60,
    90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180, 180), expand = c(0,
    0)) + theme_void() + theme(axis.title = element_blank(), axis.text = element_blank(),
    axis.ticks = element_blank(), plot.title = element_text(size = 10), legend.title = element_text(size = 8),
    legend.key.size = unit(0.2, "cm"), legend.position = "bottom") + coord_fixed(ratio = 1)

ai_2C_plot <- ggplot() + geom_raster(data = ai_2C_df, aes(y = y, x = x, fill = category)) +
    geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b",
        fill = NA, size = 0.1) + scale_fill_manual(values = arid_col, guide = guide_legend(reverse = TRUE),
    name = "Aridity Index") + ylab(NULL) + xlab(NULL) + scale_y_continuous(limits = c(-60,
    90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180, 180), expand = c(0,
    0)) + theme_void() + theme(axis.title = element_blank(), axis.text = element_blank(),
    axis.ticks = element_blank(), plot.title = element_text(size = 10)) + coord_fixed(ratio = 1)

ai_4C_plot <- ggplot() + geom_raster(data = ai_4C_df, aes(y = y, x = x, fill = category)) +
    geom_polygon(data = world_spdf, aes(x = long, y = lat, group = group), colour = "#64686b",
        fill = NA, size = 0.1) + scale_fill_manual(values = arid_col, guide = guide_legend(reverse = TRUE),
    name = "Aridity Index") + ylab(NULL) + xlab(NULL) + scale_y_continuous(limits = c(-60,
    90), expand = c(0, 0)) + scale_x_continuous(limits = c(-180, 180), expand = c(0,
    0)) + theme_void() + theme(axis.title = element_blank(), axis.text = element_blank(),
    axis.ticks = element_blank(), plot.title = element_text(size = 10)) + coord_fixed(ratio = 1)

ai_legend <- cowplot::get_legend(ai_curr_plot)

ai_prow <- cowplot::plot_grid(ai_curr_plot + theme(legend.position = "none"), ai_2C_plot +
    theme(legend.position = "none"), ai_4C_plot + theme(legend.position = "none"),
    align = "v", ncol = 3, labels = c("g", "h", "i"))

ai_plots <- cowplot::plot_grid(ai_prow, ai_legend, ncol = 1, rel_heights = c(1, 0.1))

cowplot::plot_grid(ppt_plots, pet_plots, ai_plots, ncol = 1)