---
title: "Climate_Projections"
date: "2026-02-15"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Libraries

Load 

```{r}
library(tidyverse)
library(reshape2)
library(RColorBrewer)
library(cowplot)
library(ncdf4)
library(sf)
library(rnaturalearth)
library(oce)
library(sf)
library(terra)
```



# Red Sea Temperatures

## Import data

Import and join data files.

Data is CORA sourced from CORIOLIS (https://data-selection.odatis-ocean.fr/coriolis) and was downloaded on the 16 Feb 2026.

```{r}
files <- list.files(
  path = "Depth_Data",
  pattern = "\\.csv$",   # change if needed
  full.names = TRUE
)

## join, make all columns be characters for now to ensure robustness of the join to different column types
depth_profiles <- files %>%
  map_dfr(
    ~ readr::read_csv(
        .x,
        col_types = cols(.default = col_character())
      ) %>%
      mutate(profile_id = basename(.x))
  )

## set column type for key variables - lat, lon, temperature, depth, pressure, depth, temp
depth_profiles <- depth_profiles %>%
  rename(lat = `LATITUDE (degree_north)`,
         lon = `LONGITUDE (degree_east)`,
         pres = `PRES (decibar)`,
         depth = `DEPH (meter)`,
         temp = `TEMP (degree_Celsius)`
         ) %>%
  mutate(lat == as.numeric(lat),
         lon == as.numeric(lon),
         pres == as.numeric(pres),
         depth == as.numeric(depth),
         temp == as.numeric(temp)
         )
  
```

## clean data

Remove all locations outside of the Red Sea.

```{r}
depth_profiles <- depth_profiles %>%
  filter(lon >= 32, lon <= 43.5,
         lat >= 12.54, lat <= 30)
```

estimate depth from pressure for those data with pressure only (note this is no 100% accurate because salinity does matter...)

```{r}
depth_profiles <- depth_profiles %>%
  mutate(depth = ifelse(is.na(depth) & !is.na(pres), 
                          swDepth(pressure = pres, latitude = lat, eos = "gsw"), NA
                          ))
```

filter to only data with temperature profiles and only with QC flags that are for good data (1, 2, 5, and 8)

```{r}
depth_profiles <- depth_profiles %>%
  drop_na(temp) %>%
  filter(depth > 1) %>%
  ## filter out points which are flagged as bad
  filter(POSITION_QC %in% c(1, 2, 5, 8)) %>%
  filter(TEMP_QC %in% c(1, 2, 5, 8, NA)) %>%
  filter(PRES_QC %in% c(1, 2, 5, 8, NA)) %>%
  filter(DEPH_QC %in% c(1, 2, 5, 8, NA)) %>%
  filter(if_any(c(POSITION_QC, TEMP_QC, PRES_QC, DEPH_QC), ~!is.na(.))) %>%
  group_by(profile_id) %>%
  filter(max(depth, na.rm = TRUE) > 200 &
           n() >= 10) %>%
  ungroup() %>%
  ## keep only what we need
  dplyr::select(lat, lon, temp, depth, profile_id) %>%
  ## arrange for plotting
  arrange(lat)

depth_profiles <- depth_profiles %>%
  mutate(across(c(lat, lon, temp, depth), as.numeric))

## get depth summaries
rs_temp_summary <- depth_profiles %>%
  mutate(
    depth_bin = floor(depth / 10) * 10
    ) %>%
  group_by(depth_bin) %>%
  summarise(
    temp_med = median(temp, na.rm = TRUE),
    n_profiles = n(),
    .groups = "drop")
```

plot

```{r}
pTP <-
  ggplot(data = depth_profiles, aes(y = temp, x = as.numeric(depth))) +
  geom_line(aes(group = profile_id), colour = "grey75", alpha = .1) +
  geom_smooth(colour = "black", se = FALSE) +
  scale_x_reverse() +
  labs(y = "Temperature (\u00B0C)" , x = "Depth (m)") +
  coord_flip() +
  theme_classic() +
  theme(text = element_text(size = 6), 
        legend.position = "none")
pTP
```



# Climate Projections

## Import Data

Extract data for target years and depths from the NetCDF files

Senario = CMIP6.ScenarioMIP.IPSL.IPSL-CM6A-LR.ssp585.r1i1p1f1.Odec.thetao.gn

```{r}
## files
files <- c(
  file.path("climate_projections/thetao_Odec_IPSL-CM6A-LR_ssp585_r1i1p1f1_gn_2020-2090.nc"),
  file.path("climate_projections/thetao_Odec_IPSL-CM6A-LR_ssp585_r1i1p1f1_gn_2105-2295.nc")
)

## targets
target_depths <- c(0, 200, 1000)
target_years  <- c(2020, 2050, 2100, 2200)

## helper: nearest index
nearest_index <- function(vec, value){
  which.min(abs(vec - value))
}

## prep storage
all_data <- list()

## precompute time ranges for each file
file_years <- lapply(files, function(f){
  nc <- nc_open(f)
  time_bounds <- ncvar_get(nc, "time_bounds")   # 2 x N
  time_mid <- colMeans(time_bounds)             # <-- correct: mean per column
  time_units <- ncatt_get(nc, "time", "units")$value
  ref_year <- as.numeric(substr(time_units, 12, 15))
  years <- ref_year + time_mid / 365
  nc_close(nc)
  years
})

## loop over target years
for(y in target_years){
  
  # pick the file containing the nearest year
  file_idx <- which.min(sapply(file_years, function(yrs) min(abs(yrs - y))))
  f <- files[file_idx]
  
  nc <- nc_open(f)
  
  # extract variables
  temp  <- ncvar_get(nc, "thetao")               # x,y,depth,time
  depth_bounds <- ncvar_get(nc, "olevel_bounds") 
  depth <- colMeans(depth_bounds)               # <-- correct: mean per depth level
  lon <- ncvar_get(nc, "nav_lon")
  lat <- ncvar_get(nc, "nav_lat")
  
  # nearest depth indices
  depth_idx <- sapply(target_depths, nearest_index, vec=depth)
  
  # compute years in this file
  time_bounds <- ncvar_get(nc, "time_bounds")
  time_mid <- colMeans(time_bounds)
  time_units <- ncatt_get(nc, "time", "units")$value
  ref_year <- as.numeric(substr(time_units, 12, 15))
  years_in_file <- ref_year + time_mid / 365
  
  # nearest time index
  time_idx <- nearest_index(years_in_file, y)
  
  # extract slices for each depth
  for(d in seq_along(target_depths)){
    
    t_slice <- temp[,,depth_idx[d],time_idx]
    
    df <- data.frame(
      lon   = as.vector(lon),
      lat   = as.vector(lat),
      depth = target_depths[d],
      year  = y,
      temperature = as.vector(t_slice)
    )
    
    df <- df[!is.na(df$temperature), ]
    
    all_data[[length(all_data)+1]] <- df
  }
  
  nc_close(nc)
}

## bind all data into a single data.frame
final_df <- bind_rows(all_data)
```

Now interpolate the data (files are not on regular lat/lon grids, they are curvelinear)

```{r}
interp_results <- list()

for(y in target_years){
  for(d in target_depths){
    
    sub <- final_df %>%
      filter(year == y,
             depth == d)
    
    sub <- sub[!is.na(sub$temperature), ]
    
    cat("Interpolating year:", y,
        "depth:", d, "m\n")
    
    interp_data <- with(sub,
      akima::interp(
        x = lon,
        y = lat,
        z = temperature,
        nx = 150,     # resolution (increase if desired)
        ny = 150,
        duplicate = "mean"
      )
    )
    
    interp_df <- expand.grid(
      lon = interp_data$x,
      lat = interp_data$y
    )
    
    interp_df$temperature <- as.vector(interp_data$z)
    interp_df$year  <- y
    interp_df$depth <- d
    
    interp_results[[length(interp_results)+1]] <- interp_df
  }
}

interp_all_df <- bind_rows(interp_results)
```

clean the environment

```{r}
rm(list = setdiff(ls(), c("depth_profiles", "pTP", "rs_temp_summary", "files", "target_depths", "target_years", "nearest_index", "final_df", "interp_all_df")))
```


## summarise average by depth

summarise average temperature by depth by year

```{r}
temp_summary <- final_df %>%
  group_by(year, depth) %>%
  summarise(
    temperature_mean = mean(temperature, na.rm = TRUE),
    temperature_sd   = sd(temperature, na.rm = TRUE),
    n = sum(!is.na(temperature)),
    .groups = "drop"
  ) %>%
  mutate(
    ci_lower = temperature_mean - 1.96 * temperature_sd / sqrt(n),
    ci_upper = temperature_mean + 1.96 * temperature_sd / sqrt(n)
  ) %>%
  select(year, depth, temperature_mean, ci_lower, ci_upper)
```


## summarise change from 2020

Build data frame where Red Sea temperatures are fixed at 2020 averages (global predictions do not work well for the region)
Add projected warming.

calculate projects delta T from current Red Sea temperatures at 1000m

```{r}
# Compute CMIP6 ΔT relative to 2020 for each location/depth
cmip_delta <- interp_all_df %>%
  filter(lon >= 32, lon <= 43.5,
         lat >= 12.54, lat <= 30) %>%
  group_by(lon, lat, depth) %>%
  mutate(delta_temp = temperature - temperature[year == 2020]) %>%
  ungroup()

# Apply the adjusted temperature using vectorized case_when
interp_all_adjusted <- interp_all_df %>%
  left_join(
    cmip_delta %>% select(lon, lat, depth, year, delta_temp),
    by = c("lon","lat","depth","year")
  ) %>%
  mutate(
    temperature = case_when(
      lon >= 32 & lon <= 43.5 &
      lat >= 12.54 & lat <= 30 &
      depth == 0 ~ 27.88 + delta_temp,
      lon >= 32 & lon <= 43.5 &
      lat >= 12.54 & lat <= 30 &
      depth == 200 ~ 22.10 + delta_temp,
      lon >= 32 & lon <= 43.5 &
      lat >= 12.54 & lat <= 30 &
      depth == 1000 ~ 21.52 + delta_temp,
      TRUE ~ temperature
    )
  ) %>%
  select(-delta_temp)

# calculate temperatures relative to current RS temperature at 1000m
interp_all_rs_adjusted <- interp_all_adjusted %>%
  mutate(temperature = temperature - 21.5)
```


## Plot

Plot the data.
Get shapefile for land to overlay on the maps

```{r}
land <- ne_countries(scale = "medium", returnclass = "sf")
```

Plot maps of projected temperature

```{r}
## 2020m
p1 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2020 & depth == 0), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24), breaks = c(-24, 0, 24), 
                       guide = guide_colourbar(direction = "horizontal", barheight  = 0.5, 
                                               title.position = "left", label.position = "top")) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "top", legend.title = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p2 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2020 & depth == 200), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p3 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2020 & depth == 1000), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))


## 2050
p4 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2050 & depth == 0), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p5 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2050 & depth == 200), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p6 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2050 & depth == 1000), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))


## 2100
p7 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2100 & depth == 0), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p8 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2100 & depth == 200), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p9 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2100 & depth == 1000), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))


## 2200
p10 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2200 & depth == 0), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p11 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2200 & depth == 200), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

p12 <- ggplot() +
  geom_raster(data = interp_all_rs_adjusted %>%
                filter(year == 2200 & depth == 1000), 
              aes(x=lon, y=lat, fill=temperature)) +
  geom_sf(data = land, fill = "grey75", colour = NA) +
  scale_fill_distiller(palette = "RdBu", limits = c(-24, 24)) +
  scale_x_continuous(limits = c(-179, 179)) +
  scale_y_continuous(limits = c(-75, 89), breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf() +
  theme_minimal() +
  theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
        legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

## extract legends
legend <- get_plot_component(p1 + theme(legend.text = element_text(size = 5)), "guide-box", return_all = TRUE)

## combine
pCP <- ggdraw() +
  draw_plot(plot_grid(p1 + theme(legend.position = "none"),
                      p2, p3, 
                      p4, p5, p6, 
                      p7, p8, p9, 
                      p10, p11, p12,
                      ncol = 3),
            x = .02, y = .04, width = .98, height = .9) +
  draw_text(c("0m", "200m", "1000m"), 
            x = c(.02+(.98/6), .02+(.98/6)*3, .02+(.98/6)*5), 
            y = c(.01, .01, .01), 
            size = 8, vjust = 0) +
  draw_text(c("2020", "2050", "2100", "2200"), 
            x = c(.01, .01, .01, .01), 
            y = c(.83, .605, .38, .155), 
            size = 8, angle = 90) +
  draw_text(text = "Delta*T~from~Red~Sea~at~1000*m~(degree*C)", parse = TRUE, x = .25, y = .96, size = 8) +
  draw_plot(legend[[4]], x = .02 + (.98/3), y = .96, height = .02, width = .98/3) 
pCP

rm(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, legend)
```

Plot summary of temp at depth

```{r}
pGTD <- ggplot(final_df, aes(x = as.factor(year), y = temperature, fill = as.factor(depth))) +
  geom_boxplot(outliers = FALSE) +
  geom_hline(yintercept = 27.88, linetype = 2) +
  geom_hline(yintercept = 21.5, linetype = 2) +
  scale_fill_manual(values = c("#B8F1ED", "#2BB3A8", "#0B4F57")) +
  annotate("text", x = 4.5, y = 27.88+.5, label = "Red Sea 0m (2020)", size = 5/.pt, hjust = 0, vjust = 0, angle = 90) +
  annotate("text", x = 4.5, y = 21.5-.5, label = "Red Sea 200m-3040m (2020)", size = 5/.pt, hjust = 1, vjust = 0, angle = 90) +
  #scale_x_discrete(expand = expansion(add = c(0.5, .5))) +
  labs(x = "Year", y = "Temperature under SSP5-8.5 (\u00B0C)", fill = "Depth (m)") +
  theme_classic() +
  theme(text = element_text(size = 6),
        legend.position = "inside",
    legend.position.inside = c(0.02, 0.98),
    legend.justification = c("left", "top"),

    legend.direction = "horizontal")
pGTD
```



# Plot Red Sea Map

Get Gebco data and land in higher resolution

```{r}
## Load GEBCO, filter to only water, flip to positive scale
gebcoRS <- as.data.frame(rast("GEBCO_Data/GEBCO24_RS.tif"), xy = TRUE) %>%
  filter(GEBCO24_RS <= 0) %>%
  mutate(depth = GEBCO24_RS * -1) %>%
  select(-GEBCO24_RS)

## convert to zonations
gebcoRS <- gebcoRS %>%
  mutate(
    depth_zone = case_when(
      depth < 0 ~ "Land",
      depth <= 200 ~ "0–200 m",
      depth <= 1000 ~ "200–1000 m",
      TRUE ~ ">1000 m"
    ),
    depth_zone = factor(depth_zone,
                        levels = c("0–200 m", "200–1000 m", ">1000 m"))
    )

## land
land2 <- ne_countries(scale = "large", returnclass = "sf")
```

Map of Red Sea

```{r}
pRS <- ggplot() +
  geom_tile(data = gebcoRS, aes(x = x, y = y, fill = depth_zone)) +
  geom_sf(data = land2, fill = "grey75", colour = "grey75") +
  scale_fill_manual(
    values = c("0–200 m" = "#B8F1ED", "200–1000 m" = "#2BB3A8", ">1000 m" = "#0B4F57"),
    name = "Maximum\nDepth") +
  labs(x = "Longitude", y = "Latitude") +
  coord_sf(expand = FALSE,
           xlim = c(32.2, 43.5),
           ylim = c(12.5, 30)) +
  theme_classic() +
  theme(text = element_text(size = 6),
        axis.line = element_blank(),
        legend.position = "inside",
        legend.position.inside = c(0.02, 0.02),
        legend.justification = c("left", "bottom"),
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.border = element_blank())
pRS
```



# Plot combined results

```{r}
pt <- 
  ggdraw() +
  draw_plot(pRS,
            x = 0, y = 11/18, width = 5/15, height = 7/18) +
  draw_plot(pTP,
            x = 5/15, y = 11/18, width = 3/15, height = 7/18) +
  draw_plot(pGTD,
            x = 8/15, y = 11/18, width = 7/15, height = 7/18) +
  draw_plot(pCP,
            x = 0, y = 0, width = 1, height = 11/18) +
  draw_plot_label(c("a", "b", "c", "d"), 
             x = c(.01, 5/15+.01, 8/15+.01, .01), y = c(.99, .99, .99, 11/18-.01), size = 8)

ggsave("RS_Climate.png", pt, width = 18, height = 21.5, units = "cm", dpi = 600, bg = "white")
```