Ocean Acidification: pH

Figure. Change in pH from hindcast. The map (top) shows the change in mean surface pH in the vicinity of Palau over the period 1993- 2022. The grey line is the Palau EEZ. The line plot (bottom) shows the change in mean surface pH averaged over the area within the top plot. The solid black line represents the trend, which is statistically significant (p < 0.05). The colored dots represent the 10 years with the lowest pH on record. Link to data

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import os
import os.path as op
import sys

import pandas as pd
import numpy as np
import xarray as xr
import geopandas as gpd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from myst_nb import glue 

sys.path.append("../../../../indicators_setup")

from ind_setup.plotting_int import plot_timeseries_interactive, plot_oni_index_th
from ind_setup.plotting import plot_base_map, plot_map_subplots, add_oni_cat, plot_bar_probs, fontsize


sys.path.append("../../../functions")
from data_downloaders import download_oni_index
from ocean import process_trend_with_nan

Setup

Define area of interest

#Area of interest
lon_range  = [129.4088, 137.0541]
lat_range = [1.5214, 11.6587]

EEZ shapefile

path_figs = "../../../matrix_cc/figures"

shp_f = op.join(os.getcwd(), '..', '..','..', 'data/Palau_EEZ/pw_eez_pol_april2022.shp')
shp_eez = gpd.read_file(shp_f)

Load Data

data_xr = xr.open_dataset(op.join(os.getcwd(), '..', '..','..', 'data/data_phyc_o2_ph.nc'))
dataset_id = 'ph'
label = 'pH'

Show code cell source Hide code cell source

# data_xr = xr.open_dataset(op.join(os.getcwd(), '..', '..','..', 'data/data_phyc_o2_ph_2022.nc')).isel(depth = 0).drop('depth')
# data_2025 = xr.open_dataset('/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/data/cmems_mod_glo_bgc_my_0.25deg_P1M-m_1768393001661_2024_2025.nc')
# data_xr = data_xr.merge(data_2025.sel(longitude = data_xr.longitude, latitude = data_xr.latitude).isel(depth = 0).drop('depth'))
# data_xr.to_netcdf(op.join(os.getcwd(), '..', '..','..', 'data/data_phyc_o2_ph.nc'))
# data_xr 

Analysis

Plotting

Average

fig, ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
im = ax.pcolor(data_xr.longitude, data_xr.latitude, data_xr.mean(dim='time')[dataset_id], transform=ccrs.PlateCarree(), 
                cmap = 'magma_r', 
                vmin = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 1), 
                vmax = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 99))
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
plt.colorbar(im, ax=ax, label= label)
glue("average_map", fig, display=False)
plt.savefig(op.join(path_figs, 'F14_pH_mean_map.png'), dpi=300, bbox_inches='tight')

Change

trend_m, _, _, _, _ = process_trend_with_nan(data_xr[dataset_id])

fig, ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
im = ax.pcolor(data_xr.longitude, data_xr.latitude, 
               trend_m,
                transform=ccrs.PlateCarree(), 
                cmap = 'RdBu_r', 
                vmin = -.07,
                vmax = .07,
                )
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
plt.colorbar(im, ax=ax, label=  dataset_id)
plt.savefig(op.join(path_figs, 'F14_pH_mean_map_trend.png'), dpi=300, bbox_inches='tight')

Seasonal average

data_month = data_xr.groupby('time.season').mean().sel(season = ['DJF', 'MAM', 'JJA', 'SON'])
im = plot_map_subplots(data_month, dataset_id, shp_eez = shp_eez, cmap = 'magma_r', 
                  vmin = np.nanpercentile(data_month.min(dim = 'season')[dataset_id], 1), 
                  vmax = np.nanpercentile(data_month.max(dim = 'season')[dataset_id], 99),
                  figsize = (15,11), sub_plot = [1, 4], cbar_pad = 0.05,
                  cbar = 1)

Seasonal anomaly

data_month = data_xr.groupby('time.season').mean().sel(season = ['DJF', 'MAM', 'JJA', 'SON']) - data_xr.mean(dim='time')
im = plot_map_subplots(data_month, dataset_id, shp_eez = shp_eez, 
                  cmap = 'RdBu_r', vmin=-.01, vmax=.01,
                  figsize = (15,11), sub_plot = [1, 4], cbar_pad = 0.05,
                  cbar = 1)

Annual average

data_y = data_xr.resample(time='1YE').mean()
im = plot_map_subplots(data_y, dataset_id, shp_eez = shp_eez, cmap = 'magma_r', 
                  vmin = np.nanpercentile(data_xr.min(dim = 'time')[dataset_id], 1), 
                  vmax = np.nanpercentile(data_xr.max(dim = 'time')[dataset_id], 99),
                  cbar = 1)

Annual anomaly

data_an = data_y - data_xr.mean(dim='time')
fig = plot_map_subplots(data_an, dataset_id, shp_eez = shp_eez, cmap='RdBu_r', vmin=-.1, vmax=.1, cbar = 1)

Average over area

data_month = data_xr.groupby('time.season').mean()
im = plot_map_subplots(data_month, dataset_id, shp_eez = shp_eez, cmap = 'magma_r', 
                  vmin = np.nanpercentile(data_month.min(dim = 'season')[dataset_id], 1), 
                  vmax = np.nanpercentile(data_month.max(dim = 'season')[dataset_id], 99),
                  figsize = (15,11), sub_plot = [1, 4], cbar_pad = 0.05,
                  cbar = 1)

dict_plot = [{'data' : data_xr.mean(dim = ['longitude', 'latitude']).to_dataframe(), 
              'var' : dataset_id, 'ax' : 1, 'label' : 'pH - MEAN AREA'},]

fig, trend = plot_timeseries_interactive(dict_plot, trendline=True, scatter_dict = None, return_trend=True, figsize = (25, 12));
fig.write_html(op.join(path_figs, 'F14_pH_mean_trend.html'), include_plotlyjs="cdn")

Seasonal average

fig, ax = plt.subplots(figsize=(12,2))
data_xr.mean(dim = ['longitude', 'latitude']).groupby('time.month').mean()[dataset_id].plot(ax = ax, marker = 'o', color = 'k')

[<matplotlib.lines.Line2D at 0x189395a60>]

Timeseries at a given point

loc = [7.37, 134.7]
dict_plot = [{'data' : data_xr.sel(longitude=loc[1], latitude=loc[0], method='nearest').to_dataframe(), 
              'var' : dataset_id, 'ax' : 1, 'label' : f'{label} at [{loc[0]}, {loc[1]}]'},]

fig, ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
ax.plot(loc[1], loc[0], '*', markersize = 12, color = 'royalblue', transform=ccrs.PlateCarree(), label = 'Location Analysis')
ax.legend()

<matplotlib.legend.Legend at 0x185392840>

fig = plot_timeseries_interactive(dict_plot, trendline=True, scatter_dict = None, figsize = (25, 12));

ONI index analysis

p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
df1 = download_oni_index(p_data)

lims = [-.5, .5]
plot_oni_index_th(df1, lims = lims)

Group by ONI category

df1 = add_oni_cat(df1, lims = lims)
df1['ONI'] = df1['oni_cat']
data_xr['ONI'] = (('time'), df1.iloc[np.intersect1d(data_xr.time, df1.index, return_indices=True)[2]].ONI.values)
data_xr['ONI_cat'] = (('time'), np.where(data_xr.ONI < lims[0], -1, np.where(data_xr.ONI > lims[1], 1, 0)))
data_oni = data_xr.groupby('ONI_cat').mean()

Average

fig = plot_map_subplots(data_oni, dataset_id, shp_eez = shp_eez, cmap = 'magma_r', 
                  vmin = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 1)-.005, 
                  vmax = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 99) + .005,
                  sub_plot= [1, 3], figsize = (20, 9), cbar = True, cbar_pad = 0.1,
                  titles = ['La Niña', 'Neutral', 'El Niño'],)

plt.savefig(op.join(path_figs, 'F14_pH_ENSO.png'), dpi=300, bbox_inches='tight')

Anomaly

data_an = data_oni - data_xr.mean(dim='time')
fig = plot_map_subplots(data_an, dataset_id, shp_eez = shp_eez, cmap='RdBu_r', vmin=-.005, vmax=.005,
                  sub_plot= [1, 3], figsize = (20, 9), cbar = True, cbar_pad = 0.1,
                  titles = ['La Niña', 'Neutral', 'El Niño'],)

Table

from ind_setup.tables import style_matrix, table_ocean
style_matrix(table_ocean(data_xr, trend[0], data_oni, dataset_id))

Key Metrics Summary

Metric	Value
Monthly Average	8.054
Monthly Maximum 01/09/1994	8.079
Monthly Minimum 01/05/2019	8.017
Maximum Annual Average	8.074
Minimum Annual Average	8.026
Rate of change [pH/year]	-0.001
Change between 1993 and 2025 [pH]	-0.032
Average La Niña ph	8.054
Average El Niño ph	8.052
Average Neutral ph	8.054