pH_full_figures.R

# This scripts create a standard set of ecoregion-level figures that include
# time series, maps, and warming summary plots for pH data downloaded through
# the prepare_data_exp script


# ---------- Configuration -----------


#Some initial packages
library(here)
library(terra)
library(tidyverse)
library(maps)
library(ggplot2)
library(tidyterra)
library(dplyr)

plot_dir <- "YOURFILEPATH"

# read in the sst data stored in the folder designated in DATA_DIR the
# prepare_data_exp script for each scenario
ph_data_ssp126<- readRDS("YOURFILEPATH")
ph_data_ssp370<- readRDS("YOURFILEPATH")
ph_data_ssp245<- readRDS("YOURFILEPATH")
ph_data_ssp585<- readRDS("YOURFILEPATH")

# Rename the pH column in each dataset before joining, so they don't collide
ph_data_ssp126 <- ph_data_ssp126 %>% rename(ph_ssp126 = ph)
ph_data_ssp370 <- ph_data_ssp370 %>% rename(ph_ssp370 = ph)
ph_data_ssp245 <- ph_data_ssp245 %>% rename(ph_ssp245 = ph)
ph_data_ssp585 <- ph_data_ssp585 %>% rename(ph_ssp585 = ph)


# Combine all 4 scenarios into one dataset
ph_data_spatial <- ph_data_ssp126 %>%
  left_join(ph_data_ssp370, by = c("lon", "lat", "year")) %>%
  left_join(ph_data_ssp245, by = c("lon", "lat", "year")) %>%
  left_join(ph_data_ssp585, by = c("lon", "lat", "year"))


# --- Preparing for plotting by reading in the ecoregion shape file and the world map ---

world_land <- map_data("world")

## ICES ecoregion outline (to overlay on the ph maps)
# load the ICES ecoregions, keep only the polygon, and draw its
# boundary on top of the raster so you can see the ecoregion sitting inside a
# larger ph map. The enlarged bounding box in scripts/prepare_data_exp.R
# is what provides that surrounding area.
library(sf)
ices_shp <- "YOURFILEPATH"
ices     <- st_read(ices_shp, quiet = TRUE)
norw     <- ices[ices$Ecoregion == "YOUR ECOREGION NAME", ]

# turn the polygon boundary into a plain data frame (lon, lat, group) so it can
# be drawn with geom_path alongside geom_raster + coord_quickmap (the ICES layer
# is EPSG:4326 / -180..180, matching the raster's display longitudes)
norw_coords  <- as.data.frame(st_coordinates(norw))
norw_outline <- data.frame(
  lon   = norw_coords$X,
  lat   = norw_coords$Y,
  group = interaction(norw_coords$L1, norw_coords$L2, norw_coords$L3)
)

# ---------- Time series CLIPPED to the shapefile (2015-2100) -------
# Instead of the bounding-box mean saved by prepare_data, build the annual ph
# time series using only the grid cells that fall inside the ICES ecoregion
# polygon. We do this with terra::extract():
#   1. turn each scenario's (lon, lat, year, ph) grid into a SpatRaster whose
#      layers are the years,
#   2. extract(raster, polygon, fun = mean, weights = TRUE) returns the polygon
#      mean for every layer/year. weights = TRUE area-weights cells that are only
#      partly inside the polygon, so the edge is respected instead of a hard box.
# NB: tidyr also defines extract(), so we call terra::extract() explicitly.

# terra needs a SpatVector polygon (norw above is an sf object for ggplot)
norw_vect <- vect(ices_shp)
norw_vect <- norw_vect[norw_vect$Ecoregion == "YOUR ECOREGION NAME", ]

# the 4 scenario ph columns in ph_data_spatial
ts_scenario_cols <- c("ph_ssp126", "ph_ssp370", "ph_ssp245", "ph_ssp585")

clipped_ts <- data.frame()

for (scen in ts_scenario_cols) {
  message(paste0("Clipping ", scen, " to the ecoregion polygon"))

  # one row per (lon, lat), one column per year -> becomes the raster layers
  scen_wide <- ph_data_spatial %>%
    select(lon, lat, year, ph = all_of(scen)) %>%
    pivot_wider(names_from = year, values_from = ph)

  r <- rast(scen_wide, type = "xyz", crs = "EPSG:4326")

  # area-weighted polygon mean, one value per layer (year)
  m <- terra::extract(r, norw_vect, fun = mean, na.rm = TRUE,
                      weights = TRUE, ID = FALSE)

  # m is a 1-row data frame with one column per year -> reshape to long
  scen_long <- data.frame(
    year     = as.integer(names(m)),
    mean_ph = as.numeric(m[1, ]),
    scenario = sub("ph_", "", scen)
  )

  clipped_ts <- bind_rows(clipped_ts, scen_long)
}

clipped_ts <- clipped_ts %>% filter(year >= 2015 & year <= 2100)


# absolute clipped mean ph time series
clipped_ts %>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  labs(title = "YOUR ECOREGION NAME Mean ph clipped to ICES ecoregion (2015-2100)",
       x = "Year", y = "ph (°C)", color = "Scenario") +
  theme_minimal()
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_2015_2100.png"),
       width = 10, height = 7)

clipped_ts%>%
  filter(year>2015)%>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  facet_grid(rows=vars(scenario),scales="free_y")+
  labs(title = "Mean Yearly pH Timeseries Separated (2015-2100)",
       x = "Year", y = "Temperature (C)")
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_2015_2100_separated.png"),
       width = 10, height = 7)

##Mean Yearly pH Timeseries
clipped_ts%>%
  filter(year>2015)%>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  labs(title = "Mean Yearly pH Timeseries (2015-2100)",
       x = "Year", y = "Temperature (C)")
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_2015_2100.png"), width = 10, height = 7)

##Baseline Mean Yearly pH Timeseries Separated
clipped_ts%>%
  filter(year >= 2015 & year <= 2025)%>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  facet_grid(rows=vars(scenario),scales="free_y") +
  labs(title = "Baseline (2015-2025) Mean Yearly pH Timeseries Separated",
       x = "Year", y = "Temperature (C)")
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_baseline_separated_2015_2025.png"), width = 10, height = 7)

##Baseline Mean Yearly pH Timeseries
clipped_ts%>%
  filter(year >= 2015 & year <= 2025)%>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  labs(title = "Baseline (2015-2025) Mean Yearly pH Timeseries",
       x = "Year", y = "Temperature (C)")
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_baseline_2015_2025.png"), width = 10, height = 7)

##Future Predicted Mean Yearly pH Timeseries Separated
clipped_ts%>%
  filter(year>2025)%>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  facet_grid(rows=vars(scenario),scales="free_y")+
  labs(title = "Future (2025-2100) Predicted Mean Yearly pH Timeseries Separated",
       x = "Year", y = "Temperature (C)")
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_future_separated_2025_2100.png"), width = 10, height = 7)

##Future Predicted Mean Yearly pH Timeseries
clipped_ts%>%
  filter(year>2025)%>%
  ggplot(aes(x = year, y = mean_ph, color = scenario)) +
  geom_line() +
  labs(title = "Future (2025-2100) Predicted Mean Yearly pH Timeseries",
       x = "Year", y = "Temperature (C)")
ggsave(file.path(plot_dir, "timeseries_YOUR_ECOREGION_clipped_future_2025_2100.png"), width = 10, height = 7)



# ---------- making the baseline of 2015-2025 -------------

# the 4 scenario columns to loop through
scenario_cols <- c("ph_ssp126", "ph_ssp370", "ph_ssp245", "ph_ssp585")

# start with the unique lon/lat grid; we'll attach one baseline column per scenario
baseline_spatial <- ph_data_spatial %>%
  distinct(lon, lat)

for (scen in scenario_cols) {
  message(paste0("Computing baseline for ", scen))

  baseline_col <- paste0("baseline_", scen)

  scen_baseline <- ph_data_spatial %>%
    filter(year >= 2015 & year <= 2025) %>%
    group_by(lon, lat) %>%
    summarise(!!baseline_col := mean(.data[[scen]], na.rm = TRUE),
              .groups = "drop")

  baseline_spatial <- baseline_spatial %>%
    left_join(scen_baseline, by = c("lon", "lat"))
}


# -------- calculating differences between baselines and predicted futures -----

#create a dataframe that just has past 2025 data
ph_data_spatial_future <- ph_data_spatial%>%filter(year>2025)

ph_spatial_dif_map <- ph_data_spatial%>%filter(year %in% c(2050, 2100))

# attach each point's baseline (matched by lon/lat) to every future year row,
# then compute the difference from baseline for each scenario
ph_data_spatial_future <- ph_data_spatial_future %>%
  left_join(baseline_spatial, by = c("lon", "lat"))

for (scen in scenario_cols) {
  message(paste0("Computing difference from baseline for ", scen))

  diff_col     <- paste0("diff_", scen)       # e.g. diff_ph_ssp126
  baseline_col <- paste0("baseline_", scen)   # e.g. baseline_ph_ssp126

  ph_data_spatial_future <- ph_data_spatial_future %>%
    mutate(!!diff_col := .data[[scen]] - .data[[baseline_col]])
}

# -------- mapping ---------

#creating a for loop to map all of the baselines

# the baseline columns we made above
baseline_cols <- paste0("baseline_", scenario_cols)

# stack each baseline column into one long dataframe with a scenario label,
# so we can facet by scenario in a single plot
baseline_long <- data.frame()

for (bcol in baseline_cols) {
  message(paste0("Adding ", bcol, " to the long dataframe"))

  scen_df <- baseline_spatial %>%
    select(lon, lat, ph = all_of(bcol)) %>%   # rename the current column to a common name "ph"
    mutate(scenario = bcol)                     # tag which scenario this row came from

  baseline_long <- bind_rows(baseline_long, scen_df)
}

# pull the map extent straight from the data's lon/lat points
lon_range <- range(baseline_long$lon, na.rm = TRUE)
lat_range <- range(baseline_long$lat, na.rm = TRUE)

# one faceted map: same color scale across all scenarios, one panel each
# --- ORIGINAL (tight crop, no ecoregion outline) ---
# baseline_long %>%
#   ggplot() +
#   geom_raster(aes(x = lon, y = lat, fill = ph)) +
#   geom_polygon(data = world_land, aes(x = long, y = lat, group = group),
#                fill = "gray90", color = "gray60") +
#   scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
#                        name = "ph (°C)") +
#   coord_quickmap(xlim = lon_range, ylim = lat_range, expand = FALSE) +
#   facet_wrap(~ scenario) +
#   theme_minimal() +
#   labs(title = "Baseline (2015-2025) pH by Scenario",
#        x = "Longitude", y = "Latitude")
# ggsave(file.path(plot_dir, "map_baseline_2015_2025.png"), width = 10, height = 7)

# --- Mapping---
baseline_long %>%
  ggplot() +
  geom_raster(aes(x = lon, y = lat, fill = ph)) +
  geom_polygon(data = world_land, aes(x = long, y = lat, group = group),
               fill = "gray90", color = "gray60") +
  geom_path(data = norw_outline, aes(x = lon, y = lat, group = group),
            color = "black", linewidth = 0.5) +
  scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
                       name = "ph (°C)") +
  coord_quickmap(xlim = lon_range, ylim = lat_range, expand = FALSE) +
  facet_wrap(~ scenario) +
  theme_minimal() +
  labs(title = "Baseline (2015-2025) pH by Scenario",
       x = "Longitude", y = "Latitude")
ggsave(file.path(plot_dir, "map_baseline_2015_2025.png"), width = 10, height = 7)

#mapping the future scenarios

# build a long dataframe with three panels per scenario:
#   1) the baseline ph, 2) difference from baseline at 2050, 3) difference at 2100
map_panels <- data.frame()

for (scen in scenario_cols) {
  message(paste0("Building map panels for ", scen))

  scen_label   <- sub("ph_", "", scen)       # "ssp126"
  baseline_col <- paste0("baseline_", scen)   # baseline_ph_ssp126
  diff_col     <- paste0("diff_", scen)       # diff_ph_ssp126

  # baseline (one value per lon/lat)
  base_df <- baseline_spatial %>%
    select(lon, lat, value = all_of(baseline_col)) %>%
    mutate(scenario = scen_label, panel = "baseline")

  # difference from baseline in 2050
  d2050 <- ph_data_spatial_future %>%
    filter(year == 2050) %>%
    select(lon, lat, value = all_of(diff_col)) %>%
    mutate(scenario = scen_label, panel = "diff to 2050")

  # difference from baseline in 2100
  d2100 <- ph_data_spatial_future %>%
    filter(year == 2100) %>%
    select(lon, lat, value = all_of(diff_col)) %>%
    mutate(scenario = scen_label, panel = "diff to 2100")

  map_panels <- bind_rows(map_panels, base_df, d2050, d2100)
}

# keep the panels in a sensible left-to-right order
map_panels$panel <- factor(map_panels$panel,
                           levels = c("baseline", "diff to 2050", "diff to 2100"))

# map extent from the data's lon/lat points
lon_range <- range(map_panels$lon, na.rm = TRUE)
lat_range <- range(map_panels$lat, na.rm = TRUE)

# grid of maps: one row per scenario, one column per panel
# --- ORIGINAL (tight crop, no ecoregion outline) ---
# map_panels %>%
#   ggplot() +
#   geom_raster(aes(x = lon, y = lat, fill = value)) +
#   geom_polygon(data = world_land, aes(x = long, y = lat, group = group),
#                fill = "gray90", color = "gray60") +
#   scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
#                        name = "ph / diff (°C)") +
#   coord_quickmap(xlim = lon_range, ylim = lat_range, expand = FALSE) +
#   facet_grid(scenario ~ panel) +
#   theme_minimal() +
#   labs(title = "Baseline ph and Difference from Baseline (2050, 2100) by Scenario",
#        x = "Longitude", y = "Latitude")
# ggsave(file.path(plot_dir, "map_baseline_differences_2050_2100.png"), width = 10, height = 7)

# --- NEW (larger area + ecoregion outline, SINGLE shared fill scale) ---
# map_panels %>%
#   ggplot() +
#   geom_raster(aes(x = lon, y = lat, fill = value)) +
#   geom_polygon(data = world_land, aes(x = long, y = lat, group = group),
#                fill = "gray90", color = "gray60") +
#   geom_path(data = norw_outline, aes(x = lon, y = lat, group = group),
#             color = "black", linewidth = 0.5) +
#   scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
#                        name = "pH / diff") +
#   coord_quickmap(xlim = lon_range, ylim = lat_range, expand = FALSE) +
#   facet_grid(scenario ~ panel) +
#   theme_minimal() +
#   labs(title = "Baseline pH and Difference from Baseline (2050, 2100) by Scenario",
#        x = "Longitude", y = "Latitude")
# ggsave(file.path(plot_dir, "map_baseline_differences_2050_2100.png"), width = 10, height = 7)

# --- NEW (SEPARATE fill scales: baseline pH vs difference in pH) ---
# The baseline pH (~7.5-8.1) and the difference from baseline (~ -0.4-0) sit on
# very different ranges, so we give each its own colour bar. We split into two
# sub-plots (baseline / difference) and join them with cowplot; rel_widths 1:2
# keeps every panel the same size so the rows align, and the narrow 9x8 canvas
# closes the gap between the baseline colour bar and the difference maps.
library(cowplot)

diff_map_layers <- list(
  geom_polygon(data = world_land, aes(x = long, y = lat, group = group),
               fill = "gray90", color = "gray60"),
  geom_path(data = norw_outline, aes(x = lon, y = lat, group = group),
            color = "black", linewidth = 0.5),
  coord_quickmap(xlim = lon_range, ylim = lat_range, expand = FALSE),
  theme_minimal(),
  labs(x = "Longitude", y = "Latitude")
)

# ONE baseline map (not one per scenario): the 2015-2025 baseline is essentially
# identical across scenarios (they only diverge in the future), so collapse the
# 4 scenario baselines into a single averaged baseline map for the left column.
baseline_single <- baseline_spatial %>%
  mutate(value = rowMeans(select(., all_of(baseline_cols)), na.rm = TRUE)) %>%
  select(lon, lat, value)

# left sub-plot: a SINGLE baseline pH map (no scenario faceting)
p_base_diff <- baseline_single %>%
  ggplot() +
  geom_raster(aes(x = lon, y = lat, fill = value)) +
  diff_map_layers +
  scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
                       name = "pH") +
  theme(legend.box.spacing = unit(2, "pt")) +   # colour bar hugs the map
  labs(title = "Baseline pH (2015-2025)")

# right sub-plot: difference from baseline, its own pH scale
p_diff <- map_panels %>%
  filter(panel != "baseline") %>%
  mutate(panel = droplevels(panel)) %>%
  ggplot() +
  geom_raster(aes(x = lon, y = lat, fill = value)) +
  diff_map_layers +
  scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
                       name = "diff (pH)") +
  facet_grid(scenario ~ panel) +
  labs(title = "Difference from baseline")

diff_map <- plot_grid(p_base_diff, p_diff, nrow = 1, rel_widths = c(1, 2),
                      align = "h", axis = "tb")
# narrower canvas so the height-limited panels fill the width, closing the gap
# between the baseline colour bar and the difference maps
ggsave(file.path(plot_dir, "map_baseline_differences_2050_2100.png"),
       diff_map, width = 9, height = 8, bg = "white")


# -------- third map: PERCENT change from baseline ---------
# Same idea as the difference map, but expressed as a percentage of each cell's
# 2015-2025 baseline instead of an absolute pH difference:
#   percent change = (future - baseline) / baseline * 100
# Separate fill scales again (baseline pH vs % change).

for (scen in scenario_cols) {
  message(paste0("Computing percent change from baseline for ", scen))

  pct_col      <- paste0("pct_", scen)        # e.g. pct_ph_ssp126
  diff_col     <- paste0("diff_", scen)       # diff_ph_ssp126 (future - baseline)
  baseline_col <- paste0("baseline_", scen)   # baseline_ph_ssp126

  ph_data_spatial_future <- ph_data_spatial_future %>%
    mutate(!!pct_col := .data[[diff_col]] / .data[[baseline_col]] * 100)
}

# build a long dataframe with three panels per scenario, mirroring the diff map:
#   1) the baseline pH, 2) % change at 2050, 3) % change at 2100
pct_panels <- data.frame()

for (scen in scenario_cols) {
  message(paste0("Building percent-change map panels for ", scen))

  scen_label   <- sub("ph_", "", scen)        # "ssp126"
  baseline_col <- paste0("baseline_", scen)   # baseline_ph_ssp126
  pct_col      <- paste0("pct_", scen)        # pct_ph_ssp126

  # baseline pH (one value per lon/lat) -- NB this panel is in pH units, not %
  base_df <- baseline_spatial %>%
    select(lon, lat, value = all_of(baseline_col)) %>%
    mutate(scenario = scen_label, panel = "baseline (pH)")

  # percent change from baseline in 2050
  p2050 <- ph_data_spatial_future %>%
    filter(year == 2050) %>%
    select(lon, lat, value = all_of(pct_col)) %>%
    mutate(scenario = scen_label, panel = "% change to 2050")

  # percent change from baseline in 2100
  p2100 <- ph_data_spatial_future %>%
    filter(year == 2100) %>%
    select(lon, lat, value = all_of(pct_col)) %>%
    mutate(scenario = scen_label, panel = "% change to 2100")

  pct_panels <- bind_rows(pct_panels, base_df, p2050, p2100)
}

# keep the panels in a sensible left-to-right order
pct_panels$panel <- factor(pct_panels$panel,
                           levels = c("baseline (pH)", "% change to 2050", "% change to 2100"))

# map extent from the data's lon/lat points
lon_range <- range(pct_panels$lon, na.rm = TRUE)
lat_range <- range(pct_panels$lat, na.rm = TRUE)

# common map layers, reused by both sub-plots (see diff map above for the layout
# reasoning: rel_widths 1:2 + a narrow 9x8 canvas keep the rows aligned and the
# baseline colour bar close to the percent-change maps)
pct_map_layers <- list(
  geom_polygon(data = world_land, aes(x = long, y = lat, group = group),
               fill = "gray90", color = "gray60"),
  geom_path(data = norw_outline, aes(x = lon, y = lat, group = group),
            color = "black", linewidth = 0.5),
  coord_quickmap(xlim = lon_range, ylim = lat_range, expand = FALSE),
  theme_minimal(),
  labs(x = "Longitude", y = "Latitude")
)

# left sub-plot: a SINGLE baseline pH map (no scenario faceting) --
# reuse the averaged baseline built for the difference map above
p_base_pct <- baseline_single %>%
  ggplot() +
  geom_raster(aes(x = lon, y = lat, fill = value)) +
  pct_map_layers +
  scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
                       name = "pH") +
  theme(legend.box.spacing = unit(2, "pt")) +   # colour bar hugs the map
  labs(title = "Baseline pH (2015-2025)")

# right sub-plot: percent change (its own % scale), scenario strip kept on right
p_change <- pct_panels %>%
  filter(panel != "baseline (pH)") %>%
  mutate(panel = droplevels(panel)) %>%
  ggplot() +
  geom_raster(aes(x = lon, y = lat, fill = value)) +
  pct_map_layers +
  scale_fill_gradientn(colors = c("blue", "cyan", "yellow", "orange", "red"),
                       name = "% change") +
  facet_grid(scenario ~ panel) +
  labs(title = "Percent change from baseline")

pct_map <- plot_grid(p_base_pct, p_change, nrow = 1, rel_widths = c(1, 2),
                     align = "h", axis = "tb")
ggsave(file.path(plot_dir, "map_baseline_percent_change_2050_2100.png"),
       pct_map, width = 9, height = 8, bg = "white")


#------------------ pH CHANGE SUMMARY (clipped to the shapefile) ------------------
# Equivalent of the SST warming summary, but for pH ("acidification"). The saved
# summary csv only holds the bounding-box spread and the clipped raster_*.rds
# files keep just the ensemble MEAN per cell, so ensemble p10/p90 can't be
# recovered from them. Instead we go back to the raw CMIP6 NetCDFs (which still
# carry the [lon, lat, time, member] array) and recompute using ONLY the cells
# inside the ICES polygon.
library(ncdf4)
library(lubridate)

# folder holding the raw scenario NetCDFs for this region/variable
nc_dir <- "YOURFILEPATH"

# baseline + future periods. Baseline is the 2015-2025 near-present, which lives
# in the SSP scenarios (the historical block ends in 2014).
BASELINE_PERIOD <- c(2015, 2025)
change_periods <- list(
  "Near-term (2021-2040)" = c(2021, 2040),
  "Mid-term (2041-2060)"  = c(2041, 2060),
  "Long-term (2081-2100)" = c(2081, 2100)
)

# annual, polygon-clipped ensemble stats (mean / p10 / p90) from one NetCDF.
# The pH NetCDFs store the variable as "ph" (SST used "sst").
region_ts_from_nc <- function(nc_path, scen_key, poly) {
  nc <- nc_open(nc_path); on.exit(nc_close(nc))

  lon <- ncvar_get(nc, "lon"); lat <- ncvar_get(nc, "lat")

  # read only the polygon's bounding box (+1 deg margin) instead of the globe
  e <- ext(poly)
  lon_idx <- which(lon >= e$xmin - 1 & lon <= e$xmax + 1)
  lat_idx <- which(lat >= e$ymin - 1 & lat <= e$ymax + 1)
  ph <- ncvar_get(nc, "ph",
                  start = c(lon_idx[1], lat_idx[1], 1, 1),
                  count = c(length(lon_idx), length(lat_idx), -1, -1))
  if (length(dim(ph)) == 3) dim(ph) <- c(dim(ph), 1L)   # guard: single member

  fill <- ncatt_get(nc, "ph", "_FillValue")$value
  if (!is.null(fill) && !is.na(fill)) ph[ph == fill] <- NA

  # as.numeric() strips the array dims ncvar_get keeps (otherwise the weight
  # vector is a 1-D array and matrix * array is "non-conformable")
  lon_c <- as.numeric(lon[lon_idx]); lat_c <- as.numeric(lat[lat_idx])
  n_lon <- length(lon_c); n_lat <- length(lat_c)
  n_time <- dim(ph)[3];  n_mem <- dim(ph)[4]

  # which of the clipped cells fall INSIDE the ecoregion polygon
  grid   <- expand.grid(lon = lon_c, lat = lat_c)   # lon varies fastest
  pts    <- vect(as.matrix(grid), type = "points", crs = crs(poly))
  inside <- !is.na(terra::extract(poly, pts)[, 2])

  # area weight = cos(lat) for in-polygon cells, 0 elsewhere
  wt <- as.numeric(cos(grid$lat * pi / 180) * inside)

  # area-weighted regional mean per (time, member)
  ph_mat <- matrix(ph, nrow = n_lon * n_lat, ncol = n_time * n_mem)
  num <- colSums(ph_mat * wt, na.rm = TRUE)
  den <- colSums((!is.na(ph_mat)) * wt)
  reg <- matrix(num / den, nrow = n_time, ncol = n_mem)   # [time, member]
  reg[!is.finite(reg)] <- NA_real_

  # ensemble stats across members, per time step
  ens_q   <- function(p) apply(reg, 1, quantile, probs = p, na.rm = TRUE)
  ts_mean <- rowMeans(reg, na.rm = TRUE)
  ts_p10  <- ens_q(0.10)
  ts_p90  <- ens_q(0.90)

  # annual aggregation
  tu     <- ncatt_get(nc, "time", "units")$value
  origin <- as.Date(regmatches(tu, regexpr("\\d{4}-\\d{1,2}-\\d{1,2}", tu)))
  dates  <- as.Date(ncvar_get(nc, "time"), origin = origin)
  yrs    <- year(dates)
  yr_u   <- sort(unique(yrs))
  ann    <- function(v) tapply(v, yrs, mean, na.rm = TRUE)[as.character(yr_u)]

  data.frame(
    year     = yr_u,
    scenario = scen_key,
    mean     = as.numeric(ann(ts_mean)),
    p10      = as.numeric(ann(ts_p10)),
    p90      = as.numeric(ann(ts_p90))
  )
}

# build the annual clipped time series for every SSP scenario (2015-2100)
scen_nc <- c(ssp126 = "tos_cmip6_ssp126.nc",
             ssp245 = "tos_cmip6_ssp245.nc",
             ssp370 = "tos_cmip6_ssp370.nc",
             ssp585 = "tos_cmip6_ssp585.nc")

clipped_annual <- bind_rows(lapply(names(scen_nc), function(k) {
  message("Clipping pH change series for ", k)
  region_ts_from_nc(file.path(nc_dir, scen_nc[k]), k, norw_vect)
}))

# baseline = each scenario's own clipped mean over 2015-2025 (one ref per scenario)
baseline_by_scen <- clipped_annual %>%
  filter(year >= BASELINE_PERIOD[1], year <= BASELINE_PERIOD[2]) %>%
  group_by(scenario) %>%
  summarise(ref = mean(mean, na.rm = TRUE), .groups = "drop")

# delta (change) from the 2015-2025 baseline for each future period & scenario
# (pH falls with acidification, so these deltas are negative)
ph_change_summary <- bind_rows(lapply(names(change_periods), function(pn) {
  p <- change_periods[[pn]]
  clipped_annual %>%
    filter(year >= p[1], year <= p[2]) %>%
    group_by(scenario) %>%
    summarise(
      mean_p = mean(mean, na.rm = TRUE),
      p10_p  = mean(p10,  na.rm = TRUE),
      p90_p  = mean(p90,  na.rm = TRUE),
      .groups = "drop"
    ) %>%
    left_join(baseline_by_scen, by = "scenario") %>%
    mutate(
      delta_mean = mean_p - ref,
      delta_p10  = p10_p  - ref,
      delta_p90  = p90_p  - ref,
      period     = pn
    ) %>%
    select(scenario, delta_mean, delta_p10, delta_p90, period)
})) %>%
  mutate(scenario_label = recode(scenario,
                                 ssp126 = "SSP1-2.6", ssp245 = "SSP2-4.5",
                                 ssp370 = "SSP3-7.0", ssp585 = "SSP5-8.5"))

# 12 point-ranges (3 periods x 4 scenarios), clipped to the ecoregion:
#   delta_mean -> centre point, delta_p10 -> bottom, delta_p90 -> top
ph_change_summary %>%
  mutate(
    period = factor(period, levels = c("Near-term (2021-2040)",
                                       "Mid-term (2041-2060)",
                                       "Long-term (2081-2100)")),
    scenario_label = factor(scenario_label,
                            levels = c("SSP1-2.6", "SSP2-4.5", "SSP3-7.0", "SSP5-8.5"))
  ) %>%
  ggplot(aes(x = scenario_label, color = scenario_label)) +
  geom_pointrange(
    aes(y    = delta_mean,   # centre point
        ymin = delta_p10,    # bottom of line
        ymax = delta_p90)    # top of line
  ) +
  facet_wrap(~ period) +
  labs(
    title    = " YOUR_ECOREGION projected pH change by period and scenario",
    subtitle = "Clipped to ICES ecoregion. Point = ensemble mean; line spans p10-p90 (Δ from 2015-2025 baseline)",
    x = "Scenario", y = "Change in pH", color = "Scenario"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file.path(plot_dir, "ph_change_summary_pointrange_YOUR_ECOREGION_clipped.png"),
       width = 10, height = 7)