Assessment of WOMBATlite global metics¶

In [1]:

parameters

# These first two cells must be in all notebooks!
# It allows us to run all the notebooks at once, this cell has a tag "parameters" which allows us to pass in 
# arguments externally using papermill (see mkfigs.sh for details)

# Set esm_file to the datastore for the main experiment of interest
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json"

# papermill settings. *No need to modify these if running interactively.* 
papermill = False                      # `cwd` and `nbname` will be populated by papermill.
cwd = None                             # current working directory 
nbname = None                          # notebook name

# These first two cells must be in all notebooks! # It allows us to run all the notebooks at once, this cell has a tag "parameters" which allows us to pass in # arguments externally using papermill (see mkfigs.sh for details) # Set esm_file to the datastore for the main experiment of interest esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json" # papermill settings. *No need to modify these if running interactively.* papermill = False # `cwd` and `nbname` will be populated by papermill. cwd = None # current working directory nbname = None # notebook name

In [2]:

injected-parameters

# Parameters
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/"
nbname = "wombatlite_global.ipynb"

# Parameters esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json" papermill = True cwd = "/g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/" nbname = "wombatlite_global.ipynb"

In [3]:

if not papermill: 
    import nci_ipynb, os  # requires conda/analysis3-26.03 or later
    cwd = nci_ipynb.dir()
    nbname = nci_ipynb.name()
    os.chdir(cwd)
import mkfigs_bootstrap  # noqa: adds external/access-model-mkfigs/src to sys.path (stop-gap)
from mkfigs import MkmdWriter
mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)

if not papermill: import nci_ipynb, os # requires conda/analysis3-26.03 or later cwd = nci_ipynb.dir() nbname = nci_ipynb.name() os.chdir(cwd) import mkfigs_bootstrap # noqa: adds external/access-model-mkfigs/src to sys.path (stop-gap) from mkfigs import MkmdWriter mkmd = MkmdWriter(esm_file, nbname, str(cwd), pm=papermill)

In [4]:

import string

import numpy as np
import xarray as xr
from xarray.indexes.nd_point_index import TreeAdapter
import pandas as pd
import h5py as h5
import intake
import xesmf
from scipy.io import loadmat
from scipy.stats import binned_statistic_2d, binned_statistic_dd
import sklearn

#currently installed in analysis3, needs at least 26.07
import skill_metrics as sm

from distributed import Client

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.lines as mlines
from matplotlib.gridspec import GridSpec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cmocean.cm as cmo
from cmocean.tools import lighten

import string import numpy as np import xarray as xr from xarray.indexes.nd_point_index import TreeAdapter import pandas as pd import h5py as h5 import intake import xesmf from scipy.io import loadmat from scipy.stats import binned_statistic_2d, binned_statistic_dd import sklearn #currently installed in analysis3, needs at least 26.07 import skill_metrics as sm from distributed import Client import matplotlib.pyplot as plt import matplotlib.colors as mcolors import matplotlib.lines as mlines from matplotlib.gridspec import GridSpec import cartopy.crs as ccrs import cartopy.feature as cfeature import cmocean.cm as cmo from cmocean.tools import lighten

In [5]:

import warnings
warnings.filterwarnings('ignore')

import warnings warnings.filterwarnings('ignore')

In [6]:

client = Client(threads_per_worker=1)
print(client.dashboard_link)

client = Client(threads_per_worker=1) print(client.dashboard_link)

http://127.0.0.1:8787/status

Open datastore for model data¶

In [7]:

model_datastore = intake.open_esm_datastore(
    esm_file,
    columns_with_iterables=[
            "variable",
            "variable_long_name",
            "variable_standard_name",
            "variable_cell_methods",
            "variable_units"
    ]
)

model_datastore = intake.open_esm_datastore( esm_file, columns_with_iterables=[ "variable", "variable_long_name", "variable_standard_name", "variable_cell_methods", "variable_units" ] )

In [8]:

IAF = esm_file.find('iaf') > 0

# Get the final complete year
no3 = model_datastore.search(variable="no3", frequency="1mon").to_dask(
    xarray_open_kwargs={"use_cftime": True},
)
last = no3.time[-1].astype('datetime64[s]').item()
final_year = last.year if last.month == 12 else last.year - 1

IAF = esm_file.find('iaf') > 0 # Get the final complete year no3 = model_datastore.search(variable="no3", frequency="1mon").to_dask( xarray_open_kwargs={"use_cftime": True}, ) last = no3.time[-1].astype('datetime64[s]').item() final_year = last.year if last.month == 12 else last.year - 1

Prepare gridded reference datasets¶

1 deg regular grid to interpolate to

In [9]:

grid_reg = xesmf.util.grid_global(1, 1, cf=True)

grid_reg = xesmf.util.grid_global(1, 1, cf=True)

In [10]:

def regrid_to_common(da, method="bilinear"):
    """Regrid to a common regular grid"""
    regridder = xesmf.Regridder(da, grid_reg, method, periodic=True)
    return regridder(da)

def regrid_to_common(da, method="bilinear"): """Regrid to a common regular grid""" regridder = xesmf.Regridder(da, grid_reg, method, periodic=True) return regridder(da)

Nitrate¶

From WOA 2023

In [11]:

def rename_time_to_month(da):
    """Reset the time coordinate to integer month"""
    assert len(da.time) == 12
    return da.rename({"time":"month"}).assign_coords({"month": np.arange(1,13)})

def rename_time_to_month(da): """Reset the time coordinate to integer month""" assert len(da.time) == 12 return da.rename({"time":"month"}).assign_coords({"month": np.arange(1,13)})

In [12]:

WOA_PATH = "/g/data/av17/access-nri/OM3/woa23"
files = [
    f"{WOA_PATH}/woa23_all_n{str(m+1).zfill(2)}_01.nc" for m in range(12)
]

# Get surface no3
surno3_woa = xr.open_mfdataset(
    files,
    decode_times=False,
    preprocess=lambda ds: ds["n_an"],
    compat="override",
    coords="minimal",
)["n_an"].isel(depth=0, drop=True)
surno3_woa = rename_time_to_month(surno3_woa)

# Regrid to common grid
surno3_woa_reg = regrid_to_common(surno3_woa, method="nearest_s2d").compute()

WOA_PATH = "/g/data/av17/access-nri/OM3/woa23" files = [ f"{WOA_PATH}/woa23_all_n{str(m+1).zfill(2)}_01.nc" for m in range(12) ] # Get surface no3 surno3_woa = xr.open_mfdataset( files, decode_times=False, preprocess=lambda ds: ds["n_an"], compat="override", coords="minimal", )["n_an"].isel(depth=0, drop=True) surno3_woa = rename_time_to_month(surno3_woa) # Regrid to common grid surno3_woa_reg = regrid_to_common(surno3_woa, method="nearest_s2d").compute()

Dissolved Iron¶

From Huang et al. 2022: Data-driven modeling of dissolved iron in the global ocean

In [13]:

HUANG_PATH = "/g/data/av17/access-nri/OM3/huang-2022-dFe"

# Get surface dFe
surdfe_Huang = xr.open_dataset(
    f"{HUANG_PATH}/Monthly_dFe_V2.nc"
)["dFe_RF"].isel(Depth=0, drop=True)
surdfe_Huang = surdfe_Huang.isel(Month=slice(12)).rename(
    {"Month":"month"}
)

# Regrid to common grid
surdfe_Huang_reg = regrid_to_common(surdfe_Huang, method="nearest_s2d").compute()

HUANG_PATH = "/g/data/av17/access-nri/OM3/huang-2022-dFe" # Get surface dFe surdfe_Huang = xr.open_dataset( f"{HUANG_PATH}/Monthly_dFe_V2.nc" )["dFe_RF"].isel(Depth=0, drop=True) surdfe_Huang = surdfe_Huang.isel(Month=slice(12)).rename( {"Month":"month"} ) # Regrid to common grid surdfe_Huang_reg = regrid_to_common(surdfe_Huang, method="nearest_s2d").compute()

Chlorophyll¶

From Sauzede et al. 2016: A neural network-based method for merging ocean color and Argo data to extend surface bio-optical properties to depth: Retrieval of the particulate backscattering coefficient

In [14]:

SAUZEDE_PATH = "/g/data/av17/access-nri/OM3/MULTIOBS_GLO_BIO_BGC/MULTIOBS_GLO_BIO_BGC_3D_REP_015_010/cmems_obs-mob_glo_bgc-chl-poc_my_0.25deg-climatology_P1M-m_202411"

# Get 3D chl
chl_cop = xr.open_mfdataset(
    f"{SAUZEDE_PATH}/cmems_obs-mob_glo_bgc-chl-poc_my_0.25deg-climatology_P1M-m_*_P202411.nc",
    decode_times=False
)["chl"]
chl_cop = rename_time_to_month(chl_cop)

SAUZEDE_PATH = "/g/data/av17/access-nri/OM3/MULTIOBS_GLO_BIO_BGC/MULTIOBS_GLO_BIO_BGC_3D_REP_015_010/cmems_obs-mob_glo_bgc-chl-poc_my_0.25deg-climatology_P1M-m_202411" # Get 3D chl chl_cop = xr.open_mfdataset( f"{SAUZEDE_PATH}/cmems_obs-mob_glo_bgc-chl-poc_my_0.25deg-climatology_P1M-m_*_P202411.nc", decode_times=False )["chl"] chl_cop = rename_time_to_month(chl_cop)

In [15]:

# Get surface chl
surchl_cop = chl_cop.isel(depth=0, drop=True)

# Regrid to common grid
surchl_cop_reg = regrid_to_common(surchl_cop, method="nearest_s2d").compute()

# Get surface chl surchl_cop = chl_cop.isel(depth=0, drop=True) # Regrid to common grid surchl_cop_reg = regrid_to_common(surchl_cop, method="nearest_s2d").compute()

In [16]:

def depth_of_max(da, z=None, dim="depth"):
    """
    Return the depth of the maximum, determined from parabolic interpolation
    of the three points nearest the maximum (in index space)

    If z is not provided, the "dim" coordinate values are used.
    """
    if z is None:
        z = da[dim]
    Nk = da.sizes[dim]
    
    finite = da.notnull().any(dim)
    
    k_max = da.fillna(-np.inf).argmax(dim=dim).compute()
    k_above = (k_max - 1).clip(0, Nk - 1)
    k_below = (k_max + 1).clip(0, Nk - 1)
    
    da_max = da.isel({dim: k_max}).drop_vars(dim)
    da_above = da.isel({dim: k_above}).drop_vars(dim)
    da_below = da.isel({dim: k_below}).drop_vars(dim)
    
    z_max = z.isel({dim: k_max}).drop_vars(dim)
    z_above = z.isel({dim: k_above}).drop_vars(dim)
    z_below = z.isel({dim: k_below}).drop_vars(dim)
    
    d_above = z_above - z_max
    d_below = z_below - z_max
    
    # Parabolic interpolation
    num = (da_above - da_max) * d_below**2 - (da_below - da_max) * d_above**2
    den = 2.0 * ((da_above - da_max) * d_below - (da_below - da_max) * d_above)
    
    interior = (k_max >= 1) & (k_max <= Nk - 2)
    z_peak = z_max + num / den.where(den != 0)
    z_peak = z_peak.where(interior & z_peak.notnull(), z_max)
    return z_peak.where(finite) 

def depth_of_max(da, z=None, dim="depth"): """ Return the depth of the maximum, determined from parabolic interpolation of the three points nearest the maximum (in index space) If z is not provided, the "dim" coordinate values are used. """ if z is None: z = da[dim] Nk = da.sizes[dim] finite = da.notnull().any(dim) k_max = da.fillna(-np.inf).argmax(dim=dim).compute() k_above = (k_max - 1).clip(0, Nk - 1) k_below = (k_max + 1).clip(0, Nk - 1) da_max = da.isel({dim: k_max}).drop_vars(dim) da_above = da.isel({dim: k_above}).drop_vars(dim) da_below = da.isel({dim: k_below}).drop_vars(dim) z_max = z.isel({dim: k_max}).drop_vars(dim) z_above = z.isel({dim: k_above}).drop_vars(dim) z_below = z.isel({dim: k_below}).drop_vars(dim) d_above = z_above - z_max d_below = z_below - z_max # Parabolic interpolation num = (da_above - da_max) * d_below**2 - (da_below - da_max) * d_above**2 den = 2.0 * ((da_above - da_max) * d_below - (da_below - da_max) * d_above) interior = (k_max >= 1) & (k_max <= Nk - 2) z_peak = z_max + num / den.where(den != 0) z_peak = z_peak.where(interior & z_peak.notnull(), z_max) return z_peak.where(finite)

In [17]:

# Get depth of max chl
chmaxz_cop = depth_of_max(chl_cop, dim="depth")

# Regrid to common grid
chmaxz_cop_reg = regrid_to_common(chmaxz_cop, method="nearest_s2d").compute()

# Get depth of max chl chmaxz_cop = depth_of_max(chl_cop, dim="depth") # Regrid to common grid chmaxz_cop_reg = regrid_to_common(chmaxz_cop, method="nearest_s2d").compute()

CO2 flux¶

From Chau et al 2022: A seamless ensemble-based reconstruction of surface ocean pCO2 and air–sea CO2 fluxes over the global coastal and open oceans

TODO: Need to copy

In [18]:

CHAU_PATH = "/g/data/av17/access-nri/OM3/chau-2022-carbon"

if IAF:
    co2_year = final_year
else:
    co2_year = 1990

# Get co2 surface flux
co2_cop = xr.open_mfdataset(
    f"{CHAU_PATH}/cmems_obs-mob_glo_bgc-car_my_irr-i_{co2_year}*.nc",
    decode_times=False
)["fgco2"]

# Regrid to common grid
co2_cop_reg = regrid_to_common(co2_cop, method="nearest_s2d").compute()

CHAU_PATH = "/g/data/av17/access-nri/OM3/chau-2022-carbon" if IAF: co2_year = final_year else: co2_year = 1990 # Get co2 surface flux co2_cop = xr.open_mfdataset( f"{CHAU_PATH}/cmems_obs-mob_glo_bgc-car_my_irr-i_{co2_year}*.nc", decode_times=False )["fgco2"] # Regrid to common grid co2_cop_reg = regrid_to_common(co2_cop, method="nearest_s2d").compute()

NPP (MODIS-based)¶

TODO: Need to source and copy

In [19]:

MODIS_PATH = "/g/data/ik11/observations/OSU_ocean_productivity/"

times = xr.date_range("2005-01-05", "2006-01-01", freq="8D", use_cftime=True)

with h5.File(f"{MODIS_PATH}/MODIS_Carbon_NPP.mat", 'r') as file:
    assert all(times.dayofyear == np.array(file['time_steps']).flatten())
    npp_modis = xr.DataArray(
        np.array(file['NPP']).transpose((2,1,0)),
        coords={
            "time": times,
            "lat": np.array(file['Lat'])[0,:],
            "lon": np.array(file['Lon'])[:,0],
        }
    )

MODIS_PATH = "/g/data/ik11/observations/OSU_ocean_productivity/" times = xr.date_range("2005-01-05", "2006-01-01", freq="8D", use_cftime=True) with h5.File(f"{MODIS_PATH}/MODIS_Carbon_NPP.mat", 'r') as file: assert all(times.dayofyear == np.array(file['time_steps']).flatten()) npp_modis = xr.DataArray( np.array(file['NPP']).transpose((2,1,0)), coords={ "time": times, "lat": np.array(file['Lat'])[0,:], "lon": np.array(file['Lon'])[:,0], } )

In [20]:

# Resample to monthly
npp_modis = rename_time_to_month(npp_modis.resample(time='MS').mean())

# Regrid to common grid
npp_modis_reg = regrid_to_common(npp_modis, method="nearest_s2d")

# Resample to monthly npp_modis = rename_time_to_month(npp_modis.resample(time='MS').mean()) # Regrid to common grid npp_modis_reg = regrid_to_common(npp_modis, method="nearest_s2d")

Prepare non-gridded reference datasets¶

Primary limiting nutrient¶

From Browning and Moore 2023: Global analysis of ocean phytoplankton nutrient limitation reveals high prevalence of co-limitation

In [21]:

BROWNING_PATH = "/g/data/av17/access-nri/OM3/browning-2023-bioassay"

lim_brown = pd.read_excel(
    f"{BROWNING_PATH}/Browning_Moore_Bioassay_Dataset_v1.0.xlsx",
    skiprows=1,
    usecols=(2, 3, 24, 25, 26),
    names=['lon', 'lat', 'primary', 'secondary', 'tertiary'],
)

# Encode Primary numerically (dropping Mn)
mapping = {
    'P': 1.0,
    'N_P': 1.5,
    'N': 2.0,
    'N_Fe': 2.5,
    'N_Fe_sep': 2.5,
    'Fe': 3.0,
    'Fe_Mn': 3.0,
    'Fe_Co': 3.0,
    'Fe_Zn_sep': 3.0,

}
lim_brown['primary'] = lim_brown['primary'].map(mapping)
lim_brown = lim_brown.dropna(subset=['primary'])

BROWNING_PATH = "/g/data/av17/access-nri/OM3/browning-2023-bioassay" lim_brown = pd.read_excel( f"{BROWNING_PATH}/Browning_Moore_Bioassay_Dataset_v1.0.xlsx", skiprows=1, usecols=(2, 3, 24, 25, 26), names=['lon', 'lat', 'primary', 'secondary', 'tertiary'], ) # Encode Primary numerically (dropping Mn) mapping = { 'P': 1.0, 'N_P': 1.5, 'N': 2.0, 'N_Fe': 2.5, 'N_Fe_sep': 2.5, 'Fe': 3.0, 'Fe_Mn': 3.0, 'Fe_Co': 3.0, 'Fe_Zn_sep': 3.0, } lim_brown['primary'] = lim_brown['primary'].map(mapping) lim_brown = lim_brown.dropna(subset=['primary'])

Bin the nutrient limitation onto grid for comparison

In [22]:

def bin_to_gridded(coords, values, edges, dims):
    """
    Grid scattered points onto an N-D array by per-cell mean.
    """
    sample = np.column_stack([np.asarray(c) for c in coords])
    coord_ok = np.isfinite(sample).all(axis=1)

    single = not isinstance(values, dict)
    items = {'_': values} if single else values

    out = {}
    for name, series in items.items():
        v = np.asarray(series, dtype=float)
        ok = coord_ok & np.isfinite(v)
        grid, _, _ = binned_statistic_dd(sample[ok], v[ok], statistic='mean', bins=edges)
        out[name] = xr.DataArray(grid, dims=dims)

    return out['_'] if single else out

def bin_to_gridded(coords, values, edges, dims): """ Grid scattered points onto an N-D array by per-cell mean. """ sample = np.column_stack([np.asarray(c) for c in coords]) coord_ok = np.isfinite(sample).all(axis=1) single = not isinstance(values, dict) items = {'_': values} if single else values out = {} for name, series in items.items(): v = np.asarray(series, dtype=float) ok = coord_ok & np.isfinite(v) grid, _, _ = binned_statistic_dd(sample[ok], v[ok], statistic='mean', bins=edges) out[name] = xr.DataArray(grid, dims=dims) return out['_'] if single else out

In [23]:

def bin_to_gridded(coords, values, edges, dims):
    """
    Grid scattered points onto an N-D array by per-cell mean.
    """
    sample = np.column_stack([np.asarray(c) for c in coords])
    coord_ok = np.isfinite(sample).all(axis=1)
    out = {}
    for name, series in values.items():
        v = np.asarray(series, dtype=float)
        ok = coord_ok & np.isfinite(v)
        grid, _, _ = binned_statistic_dd(sample[ok], v[ok], statistic='mean', bins=edges)
        out[name] = xr.DataArray(grid, dims=dims)
    return xr.Dataset(out)

def bin_to_gridded(coords, values, edges, dims): """ Grid scattered points onto an N-D array by per-cell mean. """ sample = np.column_stack([np.asarray(c) for c in coords]) coord_ok = np.isfinite(sample).all(axis=1) out = {} for name, series in values.items(): v = np.asarray(series, dtype=float) ok = coord_ok & np.isfinite(v) grid, _, _ = binned_statistic_dd(sample[ok], v[ok], statistic='mean', bins=edges) out[name] = xr.DataArray(grid, dims=dims) return xr.Dataset(out)

In [24]:

lon_c = grid_reg.lon
lat_c = grid_reg.lat
lon_edges = np.append(lon_c - 0.5, lon_c[-1] + 0.5)
lat_edges = np.append(lat_c - 0.5, lat_c[-1] + 0.5)

lim_brown_reg = bin_to_gridded(
    coords=[lim_brown['lat'], lim_brown['lon']],
    values={"primary": lim_brown['primary']},
    edges=[lat_edges, lon_edges],
    dims=('lat', 'lon'),
)["primary"]

print("Number of limiting nutrient measurements before binning?", lim_brown['primary'].count())
print("Number of limiting nutrient measurements after binning?", lim_brown_reg.count().values)

lon_c = grid_reg.lon lat_c = grid_reg.lat lon_edges = np.append(lon_c - 0.5, lon_c[-1] + 0.5) lat_edges = np.append(lat_c - 0.5, lat_c[-1] + 0.5) lim_brown_reg = bin_to_gridded( coords=[lim_brown['lat'], lim_brown['lon']], values={"primary": lim_brown['primary']}, edges=[lat_edges, lon_edges], dims=('lat', 'lon'), )["primary"] print("Number of limiting nutrient measurements before binning?", lim_brown['primary'].count()) print("Number of limiting nutrient measurements after binning?", lim_brown_reg.count().values)

Number of limiting nutrient measurements before binning? 942
Number of limiting nutrient measurements after binning? 137

Particle fluxes (POC, PIC)¶

From Mouw et al. 2016: Global ocean particulate organic carbon flux merged with satellite parameters

In [25]:

MOUW_PATH = "/g/data/av17/access-nri/OM3/mouw-2016-particle-flux"

partflux_raw = pd.read_table(
    f"{MOUW_PATH}/datasets/GO_flux.tab",
    skiprows=88,
    usecols=(1, 3, 4, 8, 9, 11, 15, 16, 17, 18, 19, 20, 40),
    names=[
        'Location ID', 'Latitude', 'Longitude', 'Depth (m)',
        'Datetime deployed','Duration (days)', 
        'C flux [mg/m**2/day]', 'C flux std dev',
        'POC flux [mg/m**2/day]', 'POC flux std dev', 
        'PIC flux [mg/m**2/day]', 'PIC flux std dev', 
        'Reference'
    ]
)

# Bin the months into 13 different integer categories representing
# all months (1-12) and a > 330 day average (0)
dur = partflux_raw['Duration (days)']
dep = pd.to_datetime(partflux_raw['Datetime deployed'], format='mixed', errors='coerce')
bad = dep.isna() | dur.isna()

print(f"Dropping {bad.sum()} rows with no deploy date or duration")
partflux_raw = partflux_raw[~bad].copy()
dur = dur[~bad]
dep = dep[~bad]

mid = dep + pd.to_timedelta(dur / 2, unit='D')
partflux_raw['Month'] = mid.dt.month.where(dur < 330, 0).astype('int8')

# Average replicate measurements per (month, location) — multiple months per
# location allowed
g = partflux_raw.groupby(['Month', 'Location ID'])
partflux = pd.DataFrame({
    'lon': g['Longitude'].mean(),
    'lat': g['Latitude'].mean(),
    'depth': g['Depth (m)'].mean(),
    'month': g['Month'].mean(),
    'pocflux_obs': g['POC flux [mg/m**2/day]'].mean(),
    'picflux_obs': g['PIC flux [mg/m**2/day]'].mean(),
})

MOUW_PATH = "/g/data/av17/access-nri/OM3/mouw-2016-particle-flux" partflux_raw = pd.read_table( f"{MOUW_PATH}/datasets/GO_flux.tab", skiprows=88, usecols=(1, 3, 4, 8, 9, 11, 15, 16, 17, 18, 19, 20, 40), names=[ 'Location ID', 'Latitude', 'Longitude', 'Depth (m)', 'Datetime deployed','Duration (days)', 'C flux [mg/m**2/day]', 'C flux std dev', 'POC flux [mg/m**2/day]', 'POC flux std dev', 'PIC flux [mg/m**2/day]', 'PIC flux std dev', 'Reference' ] ) # Bin the months into 13 different integer categories representing # all months (1-12) and a > 330 day average (0) dur = partflux_raw['Duration (days)'] dep = pd.to_datetime(partflux_raw['Datetime deployed'], format='mixed', errors='coerce') bad = dep.isna() | dur.isna() print(f"Dropping {bad.sum()} rows with no deploy date or duration") partflux_raw = partflux_raw[~bad].copy() dur = dur[~bad] dep = dep[~bad] mid = dep + pd.to_timedelta(dur / 2, unit='D') partflux_raw['Month'] = mid.dt.month.where(dur < 330, 0).astype('int8') # Average replicate measurements per (month, location) — multiple months per # location allowed g = partflux_raw.groupby(['Month', 'Location ID']) partflux = pd.DataFrame({ 'lon': g['Longitude'].mean(), 'lat': g['Latitude'].mean(), 'depth': g['Depth (m)'].mean(), 'month': g['Month'].mean(), 'pocflux_obs': g['POC flux [mg/m**2/day]'].mean(), 'picflux_obs': g['PIC flux [mg/m**2/day]'].mean(), })

Dropping 8 rows with no deploy date or duration

Bin the fluxes onto grid for comparison (not actually used at the moment)

# Get depth interfaces for binning datastore_filt = model_datastore.search(variable=["z_i"]) zi = datastore_filt.search( filename=datastore_filt.unique().filename[0] ).to_dask()["z_i"] month_edges = np.arange(-0.5, 13.5, 1.0) depth_edges = zi.values partflux_reg = bin_to_gridded( coords=[partflux["month"], partflux["depth"], partflux["lat"], partflux["lon"]], values={'pocflux': partflux["pocflux_obs"], 'picflux': partflux["picflux_obs"]}, edges=[month_edges, depth_edges, lat_edges, lon_edges], dims=('month', 'depth', 'lat', 'lon'), ) print("Number of POC flux measurements before binning?", partflux["pocflux_obs"].count()) print("Number of POC flux measurements after binning?", partflux_reg["pocflux"].count().values) print("Number of PIC flux measurements before binning?", partflux["picflux_obs"].count()) print("Number of PIC flux measurements after binning?", partflux_reg["picflux"].count().values)

Zooplankton grazing pressure¶

In [26]:

ZOO_GRAZE_PATH = "/g/data/ik11/observations/Zoo_Graze_OBS"

mat_data = loadmat(
    f"{ZOO_GRAZE_PATH}/GrazingPressure_Longhurst.mat"
)

provinces = [str(item[0]) for item in mat_data['Lonhurst_Prov'][0]]
total = 0.5 * (np.ravel(mat_data['GP_micro']) + np.ravel(mat_data['GP_meso']))

zoograz_long = xr.DataArray(
    total,
    coords={"province": provinces},
)

ZOO_GRAZE_PATH = "/g/data/ik11/observations/Zoo_Graze_OBS" mat_data = loadmat( f"{ZOO_GRAZE_PATH}/GrazingPressure_Longhurst.mat" ) provinces = [str(item[0]) for item in mat_data['Lonhurst_Prov'][0]] total = 0.5 * (np.ravel(mat_data['GP_micro']) + np.ravel(mat_data['GP_meso'])) zoograz_long = xr.DataArray( total, coords={"province": provinces}, )

The grazing pressures are for the following provinces

In [27]:

mat_data = loadmat(
    f"{ZOO_GRAZE_PATH}/Longhurst.mat"
)

lh_regions = mat_data['Longhurst_COARSE']
lh_lon = mat_data['LonGrid']
lh_lat = mat_data['LatGrid']
lh_provs = provinces

# Longhurst provinces
lh_regions = xr.DataArray(
    lh_regions,
    dims=("x", "y"),
    coords={"lat": (("x", "y"), lh_lat), "lon": (("x", "y"), lh_lon)},
)
p1 = (lh_regions-0.5).plot(figsize=(10,5), levels=range(12), cmap=cmo.thermal, add_colorbar=False)
cbar = plt.colorbar(p1, ticks=range(12))
_ = cbar.ax.set_yticks(np.arange(0.5,11), lh_provs)

mat_data = loadmat( f"{ZOO_GRAZE_PATH}/Longhurst.mat" ) lh_regions = mat_data['Longhurst_COARSE'] lh_lon = mat_data['LonGrid'] lh_lat = mat_data['LatGrid'] lh_provs = provinces # Longhurst provinces lh_regions = xr.DataArray( lh_regions, dims=("x", "y"), coords={"lat": (("x", "y"), lh_lat), "lon": (("x", "y"), lh_lon)}, ) p1 = (lh_regions-0.5).plot(figsize=(10,5), levels=range(12), cmap=cmo.thermal, add_colorbar=False) cbar = plt.colorbar(p1, ticks=range(12)) _ = cbar.ax.set_yticks(np.arange(0.5,11), lh_provs)

Make figure of observational targets¶

In [28]:

def add_cyclic_point(da, lon_dim="lon"):
    """Add cyclic longitude point for plotting"""
    from cartopy.util import add_cyclic_point
    
    lon_idx = da.dims.index(lon_dim)
    wrap_data, wrap_lon = add_cyclic_point(da.values, coord=da.coords[lon_dim], axis=lon_idx)

    return xr.DataArray(
        wrap_data, 
        dims=da.dims, 
        coords={**da.coords, lon_dim: wrap_lon}
    )

def add_cyclic_point(da, lon_dim="lon"): """Add cyclic longitude point for plotting""" from cartopy.util import add_cyclic_point lon_idx = da.dims.index(lon_dim) wrap_data, wrap_lon = add_cyclic_point(da.values, coord=da.coords[lon_dim], axis=lon_idx) return xr.DataArray( wrap_data, dims=da.dims, coords={**da.coords, lon_dim: wrap_lon} )

In [29]:

fig = plt.figure(figsize=(14,20))
gs = GridSpec(20, 8)
proj = ccrs.PlateCarree(central_longitude=205)
# Note Robinson can plot incorrectly 
# See https://github.com/SciTools/cartopy/issues/2457

fslab = 15; fstic = 13
axs = []; levs = []; ps = []

# 1. NO3 climatology
# =============================================
axs.append(plt.subplot(gs[0:4,0:4], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0,1,0.1), np.arange(1,10,1), np.arange(10,30.1,2))
    )
)
colmap = lighten(cmo.turbid, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(
    surno3_woa_reg.mean(dim='month')
)
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='max',
    )
)

# Add boxes around areas of interest
axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2)
axs[-1].plot((-180.0,180.0),(20,20), transform=ccrs.PlateCarree(), color='k', linewidth=2)
axs[-1].annotate(
    '',  # Text to display
    xy=(190.0,0.0),   # The arrowhead location
    xytext=(190.0,-20.0),  # The starting point of the arrow
    size=2,
    arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5),
    transform=ccrs.PlateCarree()
)
axs[-1].annotate(
    '',  # Text to display
    xy=(210.0,0.0),   # The arrowhead location
    xytext=(210.0,20.0),  # The starting point of the arrow
    size=2,
    arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5),
    transform=ccrs.PlateCarree()
)

# 2. dFe climatology
# =============================================
axs.append(plt.subplot(gs[0:4,4:8], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0.1,0.2,0.01), np.arange(0.2,1,0.05))
    )
)
colmap = lighten(cmo.turbid, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(
    surdfe_Huang_reg.mean(dim='month')
)
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='both',
    )
)

# 3. Surface chlorophyll
# =============================================
axs.append(plt.subplot(gs[4:8,0:4], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0,0.1,0.01), np.arange(0.1,0.5,0.05), np.arange(0.5,1.51,0.1))
    )
)
colmap = lighten(cmo.algae, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(
    surchl_cop_reg.mean(dim='month')
)
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='max',
    )
)

# 4. Depth of max Chl
# =============================================
axs.append(plt.subplot(gs[4:8,4:8], projection=proj))
levs.append(np.arange(0,151,5))
colmap = lighten(cmo.deep, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(chmaxz_cop_reg.mean(dim='month'))
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='max',
    )
)

# 5. Depth-integrated NPP (MODIS-based CbPM algorithm)
# =============================================
axs.append(plt.subplot(gs[8:12,0:4], projection=proj))
levs.append(np.arange(100,1001,50))
colmap = lighten(cmo.amp, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(npp_modis_reg.mean(dim='month'))
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='both',
    )
)

# 6. Downward flux of CO2
# =============================================
axs.append(plt.subplot(gs[8:12,4:8], projection=proj))
levs.append(np.arange(-4,4.1,0.4))
colmap = lighten(cmo.balance, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(-1 * co2_cop_reg.mean(dim='time'))
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='both',
    )
)

# Add box around area of interest
axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2)
axs[-1].annotate(
    '',  # Text to display
    xy=(210.0,-35.0),   # The arrowhead location
    xytext=(210.0,-20.0),  # The starting point of the arrow
    size=2,
    arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5),
    transform=ccrs.PlateCarree()
)

# 7. Downward flux of POC through ocean interior
# =============================================
axs.append(plt.subplot(gs[12:16,0:4], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0,10,1), np.arange(10,100,10), np.arange(100,1001,100))
    )
)
colmap = lighten(cmo.amp, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

partflux_cat = pd.cut(
    partflux["depth"].groupby('Location ID').mean(),
    bins=[0, 150, 500, 1000, 10000],
    labels=['Shallow', 'Upper mesopelagic', 'Lower mesopelagic', 'Deep']
)
partflux_siz = 10 * (partflux_cat.cat.codes + 1)
ps.append(
    axs[-1].scatter(
        partflux["lon"].groupby('Location ID').mean(),
        partflux["lat"].groupby('Location ID').mean(),
        c=partflux["pocflux_obs"].groupby('Location ID').mean(),
        s=partflux_siz,
        transform=ccrs.PlateCarree(),
        alpha=0.75, edgecolor='none', cmap=colmap, norm=norm
    )
)
points = [
    mlines.Line2D(
        [], [],
        color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[0]*0.2, label='Shallow'
    ),
    mlines.Line2D(
        [], [],
        color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[1]*0.2, label='Upper meso'
    ),
    mlines.Line2D(
        [], [],
        color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[2]*0.2, label='Lower meso'
    ),
    mlines.Line2D(
        [], [],
        color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[3]*0.2, label='Deep'
    )
]
axs[-1].legend(handles=points)

# 8. Downward flux of PIC through ocean interior
# =============================================
axs.append(plt.subplot(gs[12:16,4:8], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0,10,1), np.arange(10,50,5), np.arange(50,501,50))
    )
)
colmap = lighten(cmo.amp, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

ps.append(
    axs[-1].scatter(
        partflux["lon"].groupby('Location ID').mean(),
        partflux["lat"].groupby('Location ID').mean(),
        c=partflux["picflux_obs"].groupby('Location ID').mean(),
        s=partflux_siz,
        transform=ccrs.PlateCarree(),
        alpha=0.75, edgecolor='none', cmap=colmap, norm=norm
    )
)
points = [
    mlines.Line2D(
        [], [], color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[0]*0.2, label='Shallow'
    ),
    mlines.Line2D(
        [], [], color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[1]*0.2, label='Upper meso'
    ),
    mlines.Line2D(
        [], [], color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[2]*0.2, label='Lower meso'
    ),
    mlines.Line2D(
        [], [], color='w', marker='o', markeredgecolor='k', linestyle='None',
        markersize=partflux_siz.sort_values().unique()[3]*0.2, label='Deep'
    )
]
axs[-1].legend(handles=points)

# 9. Primary limiting nutrient
# =============================================
axs.append(plt.subplot(gs[16:20,0:4], projection=proj))
levs.append(None)
colmap = {2.0:'dodgerblue', 3.0:'firebrick', 2.5:'goldenrod', 1.5:'darkorchid'}
colors = [colmap[lim] for lim in lim_brown['primary']] 
ps.append(
    axs[-1].scatter(
        lim_brown['lon'],
        lim_brown['lat'],
        transform=ccrs.PlateCarree(),
        c=colors, alpha=0.75, edgecolor='k'
    )
)
points = [
    mlines.Line2D(
        [], [],
        color=colmap[2.0], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='N'),
    mlines.Line2D(
        [], [],
        color=colmap[3.0], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='Fe'),
    mlines.Line2D(
        [], [],
        color=colmap[2.5], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='N-Fe'),
    mlines.Line2D(
        [], [],
        color=colmap[1.5], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='N-P')
]
axs[-1].legend(
    handles=points,
    loc='lower center', bbox_to_anchor=(0.5,-0.2), ncols=4, frameon=False
)

# 10. Grazing pressure
# =============================================
axs.append(plt.subplot(gs[17:19,4:8]))
levs.append(None)
ps.append(
    axs[-1].scatter(zoograz_long.province[3:], zoograz_long[3:], marker='o', color='k')
)
axs[-1].tick_params(axis='x', labelrotation=90)
axs[-1].tick_params(
    left=False, labelleft=False, right=True, labelright=True, labelsize=fstic
)
axs[-1].set_ylim(0.0,0.3)

# All plots
# =============================================
title_x = 0.5; title_y = 1.05
label_x = 0.05; label_y = 1.06
titles = [
    'Surface NO$_3$ ($\mu$M)',
    'Surface dFe (nM)',
    'Surface Chlorophyll (mg m$^{-3}$)',
    'Depth of Chlorophyll max (m)',
    'Depth-integrated NPP (mg C m$^{-2}$ day$^{-1}$)', 
    'Air-sea CO$_2$ exchange (mol year$^{-1}$)',
    'Particle flux (mg C m$^{-2}$ day$^{-1}$)',
    'CaCO$_3$ flux (mg CaCO$_3$ m$^{-2}$ day$^{-1}$)',
    'Primary limiting nutrient',
    'Zooplankton grazing pressure (day$^{-1}$)'
]
for idx in range(len(axs)):
    plt.text(
        title_x, title_y, titles[idx],
        ha='center', va='center', fontweight='bold', transform=axs[idx].transAxes
    )
    plt.text(
        label_x, label_y, string.ascii_lowercase[idx],
        fontweight='bold', fontsize=fslab+2, ha='center', va='center', transform=axs[idx].transAxes
    )
    if hasattr(axs[idx], 'coastlines'):
        axs[idx].add_feature(cfeature.LAND, color='grey', zorder=3)
        axs[idx].add_feature(cfeature.COASTLINE, edgecolor='k', zorder=4)
    if levs[idx] is not None:
        plt.colorbar(
            ps[idx], ax=axs[idx], orientation='horizontal', ticks=levs[idx][::5], fraction=0.035, pad=0.05
        )
    
plt.text(
    0.5, 0.96,
    'Observations and observation-based products for model assessment',
    va='center', ha='center', fontsize=15, transform=fig.transFigure
)
plt.subplots_adjust(left=0.1, bottom=0.05, top=0.95, right=0.9, hspace=0.35, wspace=0.25)

fig = plt.figure(figsize=(14,20)) gs = GridSpec(20, 8) proj = ccrs.PlateCarree(central_longitude=205) # Note Robinson can plot incorrectly # See https://github.com/SciTools/cartopy/issues/2457 fslab = 15; fstic = 13 axs = []; levs = []; ps = [] # 1. NO3 climatology # ============================================= axs.append(plt.subplot(gs[0:4,0:4], projection=proj)) levs.append( np.concatenate( (np.arange(0,1,0.1), np.arange(1,10,1), np.arange(10,30.1,2)) ) ) colmap = lighten(cmo.turbid, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point( surno3_woa_reg.mean(dim='month') ) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='max', ) ) # Add boxes around areas of interest axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2) axs[-1].plot((-180.0,180.0),(20,20), transform=ccrs.PlateCarree(), color='k', linewidth=2) axs[-1].annotate( '', # Text to display xy=(190.0,0.0), # The arrowhead location xytext=(190.0,-20.0), # The starting point of the arrow size=2, arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5), transform=ccrs.PlateCarree() ) axs[-1].annotate( '', # Text to display xy=(210.0,0.0), # The arrowhead location xytext=(210.0,20.0), # The starting point of the arrow size=2, arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5), transform=ccrs.PlateCarree() ) # 2. dFe climatology # ============================================= axs.append(plt.subplot(gs[0:4,4:8], projection=proj)) levs.append( np.concatenate( (np.arange(0.1,0.2,0.01), np.arange(0.2,1,0.05)) ) ) colmap = lighten(cmo.turbid, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point( surdfe_Huang_reg.mean(dim='month') ) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='both', ) ) # 3. Surface chlorophyll # ============================================= axs.append(plt.subplot(gs[4:8,0:4], projection=proj)) levs.append( np.concatenate( (np.arange(0,0.1,0.01), np.arange(0.1,0.5,0.05), np.arange(0.5,1.51,0.1)) ) ) colmap = lighten(cmo.algae, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point( surchl_cop_reg.mean(dim='month') ) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='max', ) ) # 4. Depth of max Chl # ============================================= axs.append(plt.subplot(gs[4:8,4:8], projection=proj)) levs.append(np.arange(0,151,5)) colmap = lighten(cmo.deep, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point(chmaxz_cop_reg.mean(dim='month')) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='max', ) ) # 5. Depth-integrated NPP (MODIS-based CbPM algorithm) # ============================================= axs.append(plt.subplot(gs[8:12,0:4], projection=proj)) levs.append(np.arange(100,1001,50)) colmap = lighten(cmo.amp, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point(npp_modis_reg.mean(dim='month')) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='both', ) ) # 6. Downward flux of CO2 # ============================================= axs.append(plt.subplot(gs[8:12,4:8], projection=proj)) levs.append(np.arange(-4,4.1,0.4)) colmap = lighten(cmo.balance, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point(-1 * co2_cop_reg.mean(dim='time')) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='both', ) ) # Add box around area of interest axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2) axs[-1].annotate( '', # Text to display xy=(210.0,-35.0), # The arrowhead location xytext=(210.0,-20.0), # The starting point of the arrow size=2, arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5), transform=ccrs.PlateCarree() ) # 7. Downward flux of POC through ocean interior # ============================================= axs.append(plt.subplot(gs[12:16,0:4], projection=proj)) levs.append( np.concatenate( (np.arange(0,10,1), np.arange(10,100,10), np.arange(100,1001,100)) ) ) colmap = lighten(cmo.amp, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) partflux_cat = pd.cut( partflux["depth"].groupby('Location ID').mean(), bins=[0, 150, 500, 1000, 10000], labels=['Shallow', 'Upper mesopelagic', 'Lower mesopelagic', 'Deep'] ) partflux_siz = 10 * (partflux_cat.cat.codes + 1) ps.append( axs[-1].scatter( partflux["lon"].groupby('Location ID').mean(), partflux["lat"].groupby('Location ID').mean(), c=partflux["pocflux_obs"].groupby('Location ID').mean(), s=partflux_siz, transform=ccrs.PlateCarree(), alpha=0.75, edgecolor='none', cmap=colmap, norm=norm ) ) points = [ mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[0]*0.2, label='Shallow' ), mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[1]*0.2, label='Upper meso' ), mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[2]*0.2, label='Lower meso' ), mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[3]*0.2, label='Deep' ) ] axs[-1].legend(handles=points) # 8. Downward flux of PIC through ocean interior # ============================================= axs.append(plt.subplot(gs[12:16,4:8], projection=proj)) levs.append( np.concatenate( (np.arange(0,10,1), np.arange(10,50,5), np.arange(50,501,50)) ) ) colmap = lighten(cmo.amp, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) ps.append( axs[-1].scatter( partflux["lon"].groupby('Location ID').mean(), partflux["lat"].groupby('Location ID').mean(), c=partflux["picflux_obs"].groupby('Location ID').mean(), s=partflux_siz, transform=ccrs.PlateCarree(), alpha=0.75, edgecolor='none', cmap=colmap, norm=norm ) ) points = [ mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[0]*0.2, label='Shallow' ), mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[1]*0.2, label='Upper meso' ), mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[2]*0.2, label='Lower meso' ), mlines.Line2D( [], [], color='w', marker='o', markeredgecolor='k', linestyle='None', markersize=partflux_siz.sort_values().unique()[3]*0.2, label='Deep' ) ] axs[-1].legend(handles=points) # 9. Primary limiting nutrient # ============================================= axs.append(plt.subplot(gs[16:20,0:4], projection=proj)) levs.append(None) colmap = {2.0:'dodgerblue', 3.0:'firebrick', 2.5:'goldenrod', 1.5:'darkorchid'} colors = [colmap[lim] for lim in lim_brown['primary']] ps.append( axs[-1].scatter( lim_brown['lon'], lim_brown['lat'], transform=ccrs.PlateCarree(), c=colors, alpha=0.75, edgecolor='k' ) ) points = [ mlines.Line2D( [], [], color=colmap[2.0], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='N'), mlines.Line2D( [], [], color=colmap[3.0], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='Fe'), mlines.Line2D( [], [], color=colmap[2.5], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='N-Fe'), mlines.Line2D( [], [], color=colmap[1.5], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='N-P') ] axs[-1].legend( handles=points, loc='lower center', bbox_to_anchor=(0.5,-0.2), ncols=4, frameon=False ) # 10. Grazing pressure # ============================================= axs.append(plt.subplot(gs[17:19,4:8])) levs.append(None) ps.append( axs[-1].scatter(zoograz_long.province[3:], zoograz_long[3:], marker='o', color='k') ) axs[-1].tick_params(axis='x', labelrotation=90) axs[-1].tick_params( left=False, labelleft=False, right=True, labelright=True, labelsize=fstic ) axs[-1].set_ylim(0.0,0.3) # All plots # ============================================= title_x = 0.5; title_y = 1.05 label_x = 0.05; label_y = 1.06 titles = [ 'Surface NO$_3$ ($\mu$M)', 'Surface dFe (nM)', 'Surface Chlorophyll (mg m$^{-3}$)', 'Depth of Chlorophyll max (m)', 'Depth-integrated NPP (mg C m$^{-2}$ day$^{-1}$)', 'Air-sea CO$_2$ exchange (mol year$^{-1}$)', 'Particle flux (mg C m$^{-2}$ day$^{-1}$)', 'CaCO$_3$ flux (mg CaCO$_3$ m$^{-2}$ day$^{-1}$)', 'Primary limiting nutrient', 'Zooplankton grazing pressure (day$^{-1}$)' ] for idx in range(len(axs)): plt.text( title_x, title_y, titles[idx], ha='center', va='center', fontweight='bold', transform=axs[idx].transAxes ) plt.text( label_x, label_y, string.ascii_lowercase[idx], fontweight='bold', fontsize=fslab+2, ha='center', va='center', transform=axs[idx].transAxes ) if hasattr(axs[idx], 'coastlines'): axs[idx].add_feature(cfeature.LAND, color='grey', zorder=3) axs[idx].add_feature(cfeature.COASTLINE, edgecolor='k', zorder=4) if levs[idx] is not None: plt.colorbar( ps[idx], ax=axs[idx], orientation='horizontal', ticks=levs[idx][::5], fraction=0.035, pad=0.05 ) plt.text( 0.5, 0.96, 'Observations and observation-based products for model assessment', va='center', ha='center', fontsize=15, transform=fig.transFigure ) plt.subplots_adjust(left=0.1, bottom=0.05, top=0.95, right=0.9, hspace=0.35, wspace=0.25)

In [30]:

mkmd.savefig(fig,
    "Reference metrics for WOMBATlite global assessment.",
    "Observations and observation-based metrics for WOMBATlite global assessment. "
    "[GitHub issue: Evaluation: WOMBATlite global metrics](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/47)",
    dpi=150
)

mkmd.savefig(fig, "Reference metrics for WOMBATlite global assessment.", "Observations and observation-based metrics for WOMBATlite global assessment. " "[GitHub issue: Evaluation: WOMBATlite global metrics](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/47)", dpi=150 )

Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global_01.png
Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global.md

Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global.md successfully.

Prepare model output to compare with the observations¶

In [31]:

def clim_mean(da, clim=True):
    """Calculate the climatological mean over the assessment period"""
    clim_period = slice(-60, None) # Last 5 years of simulation
    return da.isel(time=clim_period).groupby("time.month").mean("time")
    
def interfaces_to_centres(da_zi):
    """
    Average interface quantities to get centre quantities
    """
    zl = (da_zi.zi.values[:-1] + da_zi.zi.values[1:]) / 2
    da_zi = da_zi.rename(zi="zl")
    lowers = da_zi.isel(zl=slice(1,None)).assign_coords({"zl": zl})
    uppers = da_zi.isel(zl=slice(None,-1)).assign_coords({"zl": zl})
    da_zl = (lowers + uppers) / 2
    da_zl.zl.attrs["axis"] = "Z"
    return da_zl

def clim_mean(da, clim=True): """Calculate the climatological mean over the assessment period""" clim_period = slice(-60, None) # Last 5 years of simulation return da.isel(time=clim_period).groupby("time.month").mean("time") def interfaces_to_centres(da_zi): """ Average interface quantities to get centre quantities """ zl = (da_zi.zi.values[:-1] + da_zi.zi.values[1:]) / 2 da_zi = da_zi.rename(zi="zl") lowers = da_zi.isel(zl=slice(1,None)).assign_coords({"zl": zl}) uppers = da_zi.isel(zl=slice(None,-1)).assign_coords({"zl": zl}) da_zl = (lowers + uppers) / 2 da_zl.zl.attrs["axis"] = "Z" return da_zl

Get model curvilinear lat/lons

In [32]:

coords = model_datastore.search(
    variable=["geolon", "geolat"], filename=".*geometry.*"
).to_dask().rename({
    "geolon": "lon",
    "geolat": "lat",
    "lonh": "xh",
    "lath": "yh"
}).compute()

coords = model_datastore.search( variable=["geolon", "geolat"], filename=".*geometry.*" ).to_dask().rename({ "geolon": "lon", "geolat": "lat", "lonh": "xh", "lath": "yh" }).compute()

Get model layer depths

In [33]:

ei = model_datastore.search(variable="e", file_id=".zi.*").to_dask()
el = interfaces_to_centres(ei)
e_om3 = xr.merge([
    ei.rename({"e": "ei"}),
    el.rename({"e": "el"})
]).assign_coords(coords)
e_om3_clim = clim_mean(e_om3)
e_om3_final_year = e_om3.sel(time=str(final_year))

ei = model_datastore.search(variable="e", file_id=".zi.*").to_dask() el = interfaces_to_centres(ei) e_om3 = xr.merge([ ei.rename({"e": "ei"}), el.rename({"e": "el"}) ]).assign_coords(coords) e_om3_clim = clim_mean(e_om3) e_om3_final_year = e_om3.sel(time=str(final_year))

Get model monthly BGC variables

In [34]:

variables = [
    "no3",
    "fe",
    "pchl",
    "npp3d",
    "dic_stf",
    "phy_lnit",
    "phy_lfer",
    "zoograzphy",
    "zoograzdet",
    "phy",
    "det",
    "det_vmove",
    "caco3",
    "caco3_vmove",
]

ds_dict = model_datastore.search(
    variable=variables, frequency="1mon"
).to_dataset_dict(progressbar=False)
bgc_om3 = xr.merge(ds_dict.values()).assign_coords(coords)
bgc_om3_clim = clim_mean(bgc_om3)
bgc_om3_final_year = bgc_om3.sel(time=str(final_year))

variables = [ "no3", "fe", "pchl", "npp3d", "dic_stf", "phy_lnit", "phy_lfer", "zoograzphy", "zoograzdet", "phy", "det", "det_vmove", "caco3", "caco3_vmove", ] ds_dict = model_datastore.search( variable=variables, frequency="1mon" ).to_dataset_dict(progressbar=False) bgc_om3 = xr.merge(ds_dict.values()).assign_coords(coords) bgc_om3_clim = clim_mean(bgc_om3) bgc_om3_final_year = bgc_om3.sel(time=str(final_year))

Nitrate¶

In [35]:

surno3_om3 = 1e6 * 1.035 * bgc_om3_clim["no3"].isel(zl=0, drop=True)

surno3_om3_reg = regrid_to_common(surno3_om3).compute()

surno3_om3 = 1e6 * 1.035 * bgc_om3_clim["no3"].isel(zl=0, drop=True) surno3_om3_reg = regrid_to_common(surno3_om3).compute()

Dissolved Iron¶

In [36]:

surdfe_om3 = 1e9 * 1.035 * bgc_om3_clim["fe"].isel(zl=0, drop=True)

surdfe_om3_reg = regrid_to_common(surdfe_om3).compute()

surdfe_om3 = 1e9 * 1.035 * bgc_om3_clim["fe"].isel(zl=0, drop=True) surdfe_om3_reg = regrid_to_common(surdfe_om3).compute()

Chlorophyll¶

In [37]:

# Compute the depth of the climatological mean, consistent with obs
chl_om3 = 1e6 * 1.035 * 12 * bgc_om3_clim["pchl"]
el_om3 = e_om3_clim["el"]

surchl_om3 = chl_om3.isel(zl=0, drop=True)
surchl_om3_reg = regrid_to_common(surchl_om3).compute()

chmaxz_om3 = depth_of_max(chl_om3, z=-1*el_om3, dim="zl")
chmaxz_om3_reg = regrid_to_common(chmaxz_om3).compute()

# Compute the depth of the climatological mean, consistent with obs chl_om3 = 1e6 * 1.035 * 12 * bgc_om3_clim["pchl"] el_om3 = e_om3_clim["el"] surchl_om3 = chl_om3.isel(zl=0, drop=True) surchl_om3_reg = regrid_to_common(surchl_om3).compute() chmaxz_om3 = depth_of_max(chl_om3, z=-1*el_om3, dim="zl") chmaxz_om3_reg = regrid_to_common(chmaxz_om3).compute()

CO2 flux¶

In [38]:

co2_om3 = 86400 * 365 * bgc_om3_final_year["dic_stf"]

co2_om3_reg = regrid_to_common(co2_om3).compute()

co2_om3 = 86400 * 365 * bgc_om3_final_year["dic_stf"] co2_om3_reg = regrid_to_common(co2_om3).compute()

NPP¶

In [39]:

# Integrate over depth
thk_om3 = -1 * e_om3["ei"].diff("zi").rename({"zi": "zl"}).assign_coords({"zl": bgc_om3.zl})
npp_om3 = 86400 * 1e6 * 1.035 * 12 * clim_mean((bgc_om3["npp3d"] * thk_om3).sum("zl"))

npp_om3_reg = regrid_to_common(npp_om3).compute()

# Integrate over depth thk_om3 = -1 * e_om3["ei"].diff("zi").rename({"zi": "zl"}).assign_coords({"zl": bgc_om3.zl}) npp_om3 = 86400 * 1e6 * 1.035 * 12 * clim_mean((bgc_om3["npp3d"] * thk_om3).sum("zl")) npp_om3_reg = regrid_to_common(npp_om3).compute()

Primary limiting nutrient¶

In [40]:

lim_om3 = xr.where(
    bgc_om3_final_year["phy_lnit"].isel(zl=0, drop=True) > 
    bgc_om3_final_year["phy_lfer"].isel(zl=0, drop=True),
    3.0, 2.0
)
lim_om3 = lim_om3.where(lim_om3 != 0).mean("time")

lim_om3_reg = regrid_to_common(lim_om3).compute()

lim_om3 = xr.where( bgc_om3_final_year["phy_lnit"].isel(zl=0, drop=True) > bgc_om3_final_year["phy_lfer"].isel(zl=0, drop=True), 3.0, 2.0 ) lim_om3 = lim_om3.where(lim_om3 != 0).mean("time") lim_om3_reg = regrid_to_common(lim_om3).compute()

Zooplankton grazing pressure¶

TODO: This just averages over the top 20 layers, which is very adhoc..

In [41]:

lh_regions_reg = regrid_to_common(lh_regions).rename("province").compute()

lh_regions_reg = regrid_to_common(lh_regions).rename("province").compute()

In [42]:

zoograz_om3 =  86400.0 * (
    (bgc_om3_final_year["zoograzphy"] + bgc_om3_final_year["zoograzdet"]) /
    bgc_om3_final_year["phy"]
).sel(zl=slice(0, 20)).mean(dim="zl")

zoograz_om3_reg = regrid_to_common(zoograz_om3).compute()

zoograz_om3 = 86400.0 * ( (bgc_om3_final_year["zoograzphy"] + bgc_om3_final_year["zoograzdet"]) / bgc_om3_final_year["phy"] ).sel(zl=slice(0, 20)).mean(dim="zl") zoograz_om3_reg = regrid_to_common(zoograz_om3).compute()

In [43]:

zoograz_long_om3 = zoograz_om3_reg.groupby(lh_regions_reg).median(...)
zoograz_long_om3 = zoograz_long_om3.assign_coords({"province": zoograz_long.province})

zoograz_long_om3 = zoograz_om3_reg.groupby(lh_regions_reg).median(...) zoograz_long_om3 = zoograz_long_om3.assign_coords({"province": zoograz_long.province})

Particle fluxes at obs locations¶

In [44]:

class SklearnGeoBallTreeAdapter(TreeAdapter):
    """
    Custom adapter for xarray.indexes.NDPointIndex using Ball tree with the haversine metric
    Taken from xoak, which is not in conda/analysis3
    """

    def __init__(self, points, options):
        from sklearn.neighbors import BallTree

        opts = dict(options)
        opts.update({"metric": "haversine"})

        self._balltree = BallTree(np.deg2rad(points), **opts)

    def query(self, points: np.ndarray):
        return self._balltree.query(np.deg2rad(points))

    def equals(self, other):
        return np.array_equal(self._balltree.data, other._balltree.data)
        
def sample_at_obs(ds, obs, check_plots=False):
    """
    Nearest-neighbour look up in ds at a provided set of points

    ds must include the layer depths, -1*el
    obs is
    Parameters
    ----------
    ds: xarray.Dataset
        Dataset containing variables to sample. Must also include the layer depths, -1*e
    obs: pandas.dataframe
        Dataframe of sample locations "lat", "lon", "depth", "month"
    check_plots: boolean
        If True, make plots comparing the model and obs points
    """
    # Set lons -180 - 180
    ds = ds.assign_coords(lon=((ds.lon + 180) % 360) - 180)

    # Set xarray index for geo-BallTree lookup
    ds = ds.set_xindex(
        ["lat", "lon"],
        xr.indexes.NDPointIndex,
        tree_adapter_cls=SklearnGeoBallTreeAdapter,
    )

    # Sample at obs points in horizontal and by month
    obs_mon = xr.DataArray(obs["month"].values, dims="obs")
    obs_lat = xr.DataArray(obs["lat"].values, dims="obs")
    obs_lon = xr.DataArray(obs["lon"].values, dims="obs")
    ds_h = ds.sel(lat=obs_lat, lon=obs_lon, method="nearest")
    ds_h = ds_h.sel(month=obs_mon)

    # Sample in depth
    obs_dep = xr.DataArray(obs["depth"].values, dims="obs")
    ind_depth = abs(ds_h["el"] - obs_dep).fillna(np.inf).argmin("zl").compute()
    ds_out = ds_h.isel(zl=ind_depth)

    if check_plots:
        dlat = ds_out.lat.values - obs_lat.values
        dlon = ds_out.lon.values - obs_lon.values
        ddep = ds_out["el"].values - obs_dep.values
        
        fig, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True)
        for ax, (name, d) in zip(axes, [
            ("Δlat (°)", dlat),
            ("Δlon (°)", dlon),
            ("Δdepth (m)", ddep)
        ]):
            ax.axhline(0, c="k", lw=0.5)
            ax.plot(d, ".", ms=2)
            ax.set_ylabel(name)
        axes[2].set_xlabel("obs index")
        fig.tight_layout()

    return ds_out.drop_vars("el")

class SklearnGeoBallTreeAdapter(TreeAdapter): """ Custom adapter for xarray.indexes.NDPointIndex using Ball tree with the haversine metric Taken from xoak, which is not in conda/analysis3 """ def __init__(self, points, options): from sklearn.neighbors import BallTree opts = dict(options) opts.update({"metric": "haversine"}) self._balltree = BallTree(np.deg2rad(points), **opts) def query(self, points: np.ndarray): return self._balltree.query(np.deg2rad(points)) def equals(self, other): return np.array_equal(self._balltree.data, other._balltree.data) def sample_at_obs(ds, obs, check_plots=False): """ Nearest-neighbour look up in ds at a provided set of points ds must include the layer depths, -1*el obs is Parameters ---------- ds: xarray.Dataset Dataset containing variables to sample. Must also include the layer depths, -1*e obs: pandas.dataframe Dataframe of sample locations "lat", "lon", "depth", "month" check_plots: boolean If True, make plots comparing the model and obs points """ # Set lons -180 - 180 ds = ds.assign_coords(lon=((ds.lon + 180) % 360) - 180) # Set xarray index for geo-BallTree lookup ds = ds.set_xindex( ["lat", "lon"], xr.indexes.NDPointIndex, tree_adapter_cls=SklearnGeoBallTreeAdapter, ) # Sample at obs points in horizontal and by month obs_mon = xr.DataArray(obs["month"].values, dims="obs") obs_lat = xr.DataArray(obs["lat"].values, dims="obs") obs_lon = xr.DataArray(obs["lon"].values, dims="obs") ds_h = ds.sel(lat=obs_lat, lon=obs_lon, method="nearest") ds_h = ds_h.sel(month=obs_mon) # Sample in depth obs_dep = xr.DataArray(obs["depth"].values, dims="obs") ind_depth = abs(ds_h["el"] - obs_dep).fillna(np.inf).argmin("zl").compute() ds_out = ds_h.isel(zl=ind_depth) if check_plots: dlat = ds_out.lat.values - obs_lat.values dlon = ds_out.lon.values - obs_lon.values ddep = ds_out["el"].values - obs_dep.values fig, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True) for ax, (name, d) in zip(axes, [ ("Δlat (°)", dlat), ("Δlon (°)", dlon), ("Δdepth (m)", ddep) ]): ax.axhline(0, c="k", lw=0.5) ax.plot(d, ".", ms=2) ax.set_ylabel(name) axes[2].set_xlabel("obs index") fig.tight_layout() return ds_out.drop_vars("el")

In [45]:

partflux_om3 = xr.merge([
    (bgc_om3_final_year["det"] * bgc_om3_final_year["det_vmove"]).rename("pocflux"),
    (bgc_om3_final_year["caco3"] * bgc_om3_final_year["caco3_vmove"]).rename("picflux"),
    -1 * e_om3_final_year["el"]
])
    
# Prefix annual mean and relabel time coordinate to match obs partflux months
partflux_om3 = xr.concat((
    partflux_om3.mean("time").assign_coords({"month": 0}),
    partflux_om3.rename({"time":"month"}).assign_coords({"month": np.arange(1,13)})
),
    dim="month"
)

# Sample at observation points
partflux_om3 = sample_at_obs(partflux_om3, partflux, check_plots=True)

# Add to dataframe
partflux["pocflux_model"] = 86400 * 1e6 * 1.035 * 12 * partflux_om3["pocflux"].values
partflux["picflux_model"] = 86400 * 1e6 * 1.035 * 12 * partflux_om3["picflux"].values

partflux_om3 = xr.merge([ (bgc_om3_final_year["det"] * bgc_om3_final_year["det_vmove"]).rename("pocflux"), (bgc_om3_final_year["caco3"] * bgc_om3_final_year["caco3_vmove"]).rename("picflux"), -1 * e_om3_final_year["el"] ]) # Prefix annual mean and relabel time coordinate to match obs partflux months partflux_om3 = xr.concat(( partflux_om3.mean("time").assign_coords({"month": 0}), partflux_om3.rename({"time":"month"}).assign_coords({"month": np.arange(1,13)}) ), dim="month" ) # Sample at observation points partflux_om3 = sample_at_obs(partflux_om3, partflux, check_plots=True) # Add to dataframe partflux["pocflux_model"] = 86400 * 1e6 * 1.035 * 12 * partflux_om3["pocflux"].values partflux["picflux_model"] = 86400 * 1e6 * 1.035 * 12 * partflux_om3["picflux"].values

Make figure of model data¶

In [46]:

fig = plt.figure(figsize=(14,20))
gs = GridSpec(20, 8)
proj = ccrs.PlateCarree(central_longitude=205)
# Note Robinson can plot incorrectly 
# See https://github.com/SciTools/cartopy/issues/2457

fslab = 15; fstic = 13
axs = []; levs = []; ps = []

# 1. NO3 climatology
# =============================================
axs.append(plt.subplot(gs[0:4,0:4], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0,1,0.1), np.arange(1,10,1), np.arange(10,30.1,2))
    )
)
colmap = lighten(cmo.turbid, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(
    surno3_om3_reg.mean(dim='month')
)
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='max'
    )
)

# Add boxes around areas of interest
axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2)
axs[-1].plot((-180.0,180.0),(20,20), transform=ccrs.PlateCarree(), color='k', linewidth=2)
axs[-1].annotate(
    '',  # Text to display
    xy=(190.0,0.0),   # The arrowhead location
    xytext=(190.0,-20.0),  # The starting point of the arrow
    size=2,
    arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5),
    transform=ccrs.PlateCarree()
)
axs[-1].annotate(
    '',  # Text to display
    xy=(210.0,0.0),   # The arrowhead location
    xytext=(210.0,20.0),  # The starting point of the arrow
    size=2,
    arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5),
    transform=ccrs.PlateCarree()
)

# 2. dFe climatology
# =============================================
axs.append(plt.subplot(gs[0:4,4:8], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0.1,0.2,0.01), np.arange(0.2,1,0.05))
    )
)
colmap = lighten(cmo.turbid, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(
    surdfe_om3_reg.mean(dim='month')
)
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='both'
    )
)

# 3. Surface chlorophyll
# =============================================
axs.append(plt.subplot(gs[4:8,0:4], projection=proj))
levs.append(
    np.concatenate(
        (np.arange(0,0.1,0.01), np.arange(0.1,0.5,0.05), np.arange(0.5,1.51,0.1))
    )
)
colmap = lighten(cmo.algae, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(
    surchl_om3_reg.mean(dim='month')
)
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='max'
    )
)

# 4. Depth of max Chl
# =============================================
axs.append(plt.subplot(gs[4:8,4:8], projection=proj))
levs.append(np.arange(0,151,5))
colmap = lighten(cmo.deep, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(chmaxz_om3_reg.mean(dim='month'))
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='max'
    )
)

# 5. Depth-integrated NPP (MODIS-based CbPM algorithm)
# =============================================
axs.append(plt.subplot(gs[8:12,0:4], projection=proj))
levs.append(np.arange(100,1001,50))
colmap = lighten(cmo.amp, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(npp_om3_reg.mean(dim='month'))
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='both'
    )
)

# 6. Downward flux of CO2
# =============================================
axs.append(plt.subplot(gs[8:12,4:8], projection=proj))
levs.append(np.arange(-4,4.1,0.4))
colmap = lighten(cmo.balance, 0.8)
norm = mcolors.BoundaryNorm(levs[-1], ncolors=256)

to_plot = add_cyclic_point(-1 * co2_om3_reg.mean(dim='time'))
ps.append(
    to_plot.plot.contourf(
        ax=axs[-1],
        transform=ccrs.PlateCarree(),
        add_colorbar=False,
        cmap=colmap, norm=norm, levels=levs[-1], extend='both'
    )
)

# Add box around area of interest
axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2)
axs[-1].annotate(
    '',  # Text to display
    xy=(210.0,-35.0),   # The arrowhead location
    xytext=(210.0,-20.0),  # The starting point of the arrow
    size=2,
    arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5),
    transform=ccrs.PlateCarree()
)

# 7. Downward flux of POC through ocean interior
# =============================================
axs.append(plt.subplot(gs[13:15,0:4]))
levs.append(None)
colmap = cmo.deep
ps.append(
    axs[-1].scatter(
        partflux['pocflux_obs'], partflux['pocflux_model'], c=partflux['depth'],
        marker='o', cmap=colmap, alpha=0.25, s=20, edgecolor='grey'
    )
)
axs[-1].set_xscale('log')
axs[-1].set_yscale('log')
axs[-1].set_xlim(1e-2, 1e3)
axs[-1].set_ylabel('model', fontsize=fstic)
axs[-1].set_xlabel('observations', fontsize=fstic)
axs[-1].plot((1e-5, 1e3), (1e-5, 1e3), 'k--')

# 8. Downward flux of PIC through ocean interior
# =============================================
axs.append(plt.subplot(gs[13:15,4:8]))
levs.append(None)
colmap = cmo.deep
ps.append(
    axs[-1].scatter(
        partflux['picflux_obs'], partflux['picflux_model'], c=partflux['depth'],
        marker='o', cmap=colmap, alpha=0.25, s=20, edgecolor='grey'
    )
)
axs[-1].set_xscale('log')
axs[-1].set_yscale('log')
axs[-1].set_xlim(1e-2, 1e3)
axs[-1].set_xlabel('observations', fontsize=fstic)
axs[-1].plot((1e-5, 1e3), (1e-5, 1e3), 'k--')
axs[-1].tick_params(right=True, left=False, labelleft=False, labelright=True)

# 9. Primary limiting nutrient
# =============================================
axs.append(plt.subplot(gs[16:20,0:4], projection=proj))
levs.append(None)
colmap = {2.0:'dodgerblue', 3.0:'firebrick', 2.5:'goldenrod', 1.5:'darkorchid'}
colors = [colmap[lim] for lim in lim_brown['primary']]
to_plot = add_cyclic_point(lim_om3_reg)
to_plot.plot.contourf(
    ax=axs[-1],
    transform=ccrs.PlateCarree(),
    add_colorbar=False,
    cmap="RdBu_r", levels=np.array([0,1e-6,2.1,2.1+1e-6,3.5]), extend='neither',
)
ps.append(
    axs[-1].scatter(
        lim_brown['lon'],
        lim_brown['lat'],
        transform=ccrs.PlateCarree(),
        c=colors, alpha=0.75, edgecolor='k'
    )
)
points = [
    mlines.Line2D(
        [], [],
        color=colmap[2.0], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='N'),
    mlines.Line2D(
        [], [],
        color=colmap[3.0], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='Fe'),
    mlines.Line2D(
        [], [],
        color=colmap[2.5], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='N-Fe'),
    mlines.Line2D(
        [], [],
        color=colmap[1.5], marker='o', markeredgecolor='k', linestyle='None',
        markersize=10, label='N-P')
]
axs[-1].legend(
    handles=points,
    loc='lower center', bbox_to_anchor=(0.5,-0.2), ncols=4, frameon=False
)

# # 10. Grazing pressure
# # =============================================
axs.append(plt.subplot(gs[17:19,4:8]))
levs.append(None)
axs[-1].scatter(
    zoograz_long.province[3:], zoograz_long[3:],
    marker='o', color='k', label="Rohr et al. (2024)"
)

ps.append(
    axs[-1].scatter(
        zoograz_long_om3.province[3:], zoograz_long_om3[3:],
        marker='o', color='firebrick', label="ACCESS-OM3-WOMBATlite"
    )
)
axs[-1].tick_params(axis='x', labelrotation=90)
axs[-1].tick_params(
    left=False, labelleft=False, right=True, labelright=True, labelsize=fstic
)
axs[-1].set_ylim(-0.025,0.3)

# All plots
# =============================================
title_x = 0.5; title_y = 1.05
label_x = 0.05; label_y = 1.06
titles = [
    'Surface NO$_3$ ($\mu$M)',
    'Surface dFe (nM)',
    'Surface Chlorophyll (mg m$^{-3}$)',
    'Depth of Chlorophyll max (m)',
    'Depth-integrated NPP (mg C m$^{-2}$ day$^{-1}$)', 
    'Air-sea CO$_2$ exchange (mol year$^{-1}$)',
    'Particle flux (mg C m$^{-2}$ day$^{-1}$)',
    'CaCO$_3$ flux (mg CaCO$_3$ m$^{-2}$ day$^{-1}$)',
    'Primary limiting nutrient',
    'Zooplankton grazing pressure (day$^{-1}$)'
]
for idx in range(len(axs)):
    plt.text(
        title_x, title_y, titles[idx],
        ha='center', va='center', fontweight='bold', transform=axs[idx].transAxes
    )
    plt.text(
        label_x, label_y, string.ascii_lowercase[idx],
        fontweight='bold', fontsize=fslab+2, ha='center', va='center', transform=axs[idx].transAxes
    )
    if hasattr(axs[idx], 'coastlines'):
        axs[idx].add_feature(cfeature.LAND, color='grey', zorder=3)
        axs[idx].add_feature(cfeature.COASTLINE, edgecolor='k', zorder=4)
    if levs[idx] is not None:
        plt.colorbar(
            ps[idx], ax=axs[idx], orientation='horizontal', ticks=levs[idx][::5], fraction=0.035, pad=0.05
        )
    
plt.text(
    0.5, 0.96,
    'ACCESS-OM3-WOMBATlite model assessment',
    va='center', ha='center', fontsize=15, transform=fig.transFigure
)
plt.subplots_adjust(left=0.1, bottom=0.05, top=0.95, right=0.9, hspace=0.35, wspace=0.25)

# Colorbar for particle fluxes
# =============================================
cbax8a = fig.add_axes([0.4, 0.23, 0.2, 0.01])
cbar8a = plt.colorbar(ps[7], cax=cbax8a, orientation='horizontal')
cbar8a.ax.tick_params(labelbottom=False, bottom=False, labeltop=True, top=True)
_ = plt.text(0.5, -0.54, 'Depth (m)', fontsize=fstic, ha='center', va='center', transform=cbax8a.transAxes)

fig = plt.figure(figsize=(14,20)) gs = GridSpec(20, 8) proj = ccrs.PlateCarree(central_longitude=205) # Note Robinson can plot incorrectly # See https://github.com/SciTools/cartopy/issues/2457 fslab = 15; fstic = 13 axs = []; levs = []; ps = [] # 1. NO3 climatology # ============================================= axs.append(plt.subplot(gs[0:4,0:4], projection=proj)) levs.append( np.concatenate( (np.arange(0,1,0.1), np.arange(1,10,1), np.arange(10,30.1,2)) ) ) colmap = lighten(cmo.turbid, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point( surno3_om3_reg.mean(dim='month') ) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='max' ) ) # Add boxes around areas of interest axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2) axs[-1].plot((-180.0,180.0),(20,20), transform=ccrs.PlateCarree(), color='k', linewidth=2) axs[-1].annotate( '', # Text to display xy=(190.0,0.0), # The arrowhead location xytext=(190.0,-20.0), # The starting point of the arrow size=2, arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5), transform=ccrs.PlateCarree() ) axs[-1].annotate( '', # Text to display xy=(210.0,0.0), # The arrowhead location xytext=(210.0,20.0), # The starting point of the arrow size=2, arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5), transform=ccrs.PlateCarree() ) # 2. dFe climatology # ============================================= axs.append(plt.subplot(gs[0:4,4:8], projection=proj)) levs.append( np.concatenate( (np.arange(0.1,0.2,0.01), np.arange(0.2,1,0.05)) ) ) colmap = lighten(cmo.turbid, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point( surdfe_om3_reg.mean(dim='month') ) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='both' ) ) # 3. Surface chlorophyll # ============================================= axs.append(plt.subplot(gs[4:8,0:4], projection=proj)) levs.append( np.concatenate( (np.arange(0,0.1,0.01), np.arange(0.1,0.5,0.05), np.arange(0.5,1.51,0.1)) ) ) colmap = lighten(cmo.algae, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point( surchl_om3_reg.mean(dim='month') ) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='max' ) ) # 4. Depth of max Chl # ============================================= axs.append(plt.subplot(gs[4:8,4:8], projection=proj)) levs.append(np.arange(0,151,5)) colmap = lighten(cmo.deep, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point(chmaxz_om3_reg.mean(dim='month')) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='max' ) ) # 5. Depth-integrated NPP (MODIS-based CbPM algorithm) # ============================================= axs.append(plt.subplot(gs[8:12,0:4], projection=proj)) levs.append(np.arange(100,1001,50)) colmap = lighten(cmo.amp, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point(npp_om3_reg.mean(dim='month')) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='both' ) ) # 6. Downward flux of CO2 # ============================================= axs.append(plt.subplot(gs[8:12,4:8], projection=proj)) levs.append(np.arange(-4,4.1,0.4)) colmap = lighten(cmo.balance, 0.8) norm = mcolors.BoundaryNorm(levs[-1], ncolors=256) to_plot = add_cyclic_point(-1 * co2_om3_reg.mean(dim='time')) ps.append( to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap=colmap, norm=norm, levels=levs[-1], extend='both' ) ) # Add box around area of interest axs[-1].plot((-180.0,180.0),(-20,-20), transform=ccrs.PlateCarree(), color='k', linewidth=2) axs[-1].annotate( '', # Text to display xy=(210.0,-35.0), # The arrowhead location xytext=(210.0,-20.0), # The starting point of the arrow size=2, arrowprops=dict(facecolor='k', headwidth=6, shrink=0.0, headlength=5, width=1.5), transform=ccrs.PlateCarree() ) # 7. Downward flux of POC through ocean interior # ============================================= axs.append(plt.subplot(gs[13:15,0:4])) levs.append(None) colmap = cmo.deep ps.append( axs[-1].scatter( partflux['pocflux_obs'], partflux['pocflux_model'], c=partflux['depth'], marker='o', cmap=colmap, alpha=0.25, s=20, edgecolor='grey' ) ) axs[-1].set_xscale('log') axs[-1].set_yscale('log') axs[-1].set_xlim(1e-2, 1e3) axs[-1].set_ylabel('model', fontsize=fstic) axs[-1].set_xlabel('observations', fontsize=fstic) axs[-1].plot((1e-5, 1e3), (1e-5, 1e3), 'k--') # 8. Downward flux of PIC through ocean interior # ============================================= axs.append(plt.subplot(gs[13:15,4:8])) levs.append(None) colmap = cmo.deep ps.append( axs[-1].scatter( partflux['picflux_obs'], partflux['picflux_model'], c=partflux['depth'], marker='o', cmap=colmap, alpha=0.25, s=20, edgecolor='grey' ) ) axs[-1].set_xscale('log') axs[-1].set_yscale('log') axs[-1].set_xlim(1e-2, 1e3) axs[-1].set_xlabel('observations', fontsize=fstic) axs[-1].plot((1e-5, 1e3), (1e-5, 1e3), 'k--') axs[-1].tick_params(right=True, left=False, labelleft=False, labelright=True) # 9. Primary limiting nutrient # ============================================= axs.append(plt.subplot(gs[16:20,0:4], projection=proj)) levs.append(None) colmap = {2.0:'dodgerblue', 3.0:'firebrick', 2.5:'goldenrod', 1.5:'darkorchid'} colors = [colmap[lim] for lim in lim_brown['primary']] to_plot = add_cyclic_point(lim_om3_reg) to_plot.plot.contourf( ax=axs[-1], transform=ccrs.PlateCarree(), add_colorbar=False, cmap="RdBu_r", levels=np.array([0,1e-6,2.1,2.1+1e-6,3.5]), extend='neither', ) ps.append( axs[-1].scatter( lim_brown['lon'], lim_brown['lat'], transform=ccrs.PlateCarree(), c=colors, alpha=0.75, edgecolor='k' ) ) points = [ mlines.Line2D( [], [], color=colmap[2.0], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='N'), mlines.Line2D( [], [], color=colmap[3.0], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='Fe'), mlines.Line2D( [], [], color=colmap[2.5], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='N-Fe'), mlines.Line2D( [], [], color=colmap[1.5], marker='o', markeredgecolor='k', linestyle='None', markersize=10, label='N-P') ] axs[-1].legend( handles=points, loc='lower center', bbox_to_anchor=(0.5,-0.2), ncols=4, frameon=False ) # # 10. Grazing pressure # # ============================================= axs.append(plt.subplot(gs[17:19,4:8])) levs.append(None) axs[-1].scatter( zoograz_long.province[3:], zoograz_long[3:], marker='o', color='k', label="Rohr et al. (2024)" ) ps.append( axs[-1].scatter( zoograz_long_om3.province[3:], zoograz_long_om3[3:], marker='o', color='firebrick', label="ACCESS-OM3-WOMBATlite" ) ) axs[-1].tick_params(axis='x', labelrotation=90) axs[-1].tick_params( left=False, labelleft=False, right=True, labelright=True, labelsize=fstic ) axs[-1].set_ylim(-0.025,0.3) # All plots # ============================================= title_x = 0.5; title_y = 1.05 label_x = 0.05; label_y = 1.06 titles = [ 'Surface NO$_3$ ($\mu$M)', 'Surface dFe (nM)', 'Surface Chlorophyll (mg m$^{-3}$)', 'Depth of Chlorophyll max (m)', 'Depth-integrated NPP (mg C m$^{-2}$ day$^{-1}$)', 'Air-sea CO$_2$ exchange (mol year$^{-1}$)', 'Particle flux (mg C m$^{-2}$ day$^{-1}$)', 'CaCO$_3$ flux (mg CaCO$_3$ m$^{-2}$ day$^{-1}$)', 'Primary limiting nutrient', 'Zooplankton grazing pressure (day$^{-1}$)' ] for idx in range(len(axs)): plt.text( title_x, title_y, titles[idx], ha='center', va='center', fontweight='bold', transform=axs[idx].transAxes ) plt.text( label_x, label_y, string.ascii_lowercase[idx], fontweight='bold', fontsize=fslab+2, ha='center', va='center', transform=axs[idx].transAxes ) if hasattr(axs[idx], 'coastlines'): axs[idx].add_feature(cfeature.LAND, color='grey', zorder=3) axs[idx].add_feature(cfeature.COASTLINE, edgecolor='k', zorder=4) if levs[idx] is not None: plt.colorbar( ps[idx], ax=axs[idx], orientation='horizontal', ticks=levs[idx][::5], fraction=0.035, pad=0.05 ) plt.text( 0.5, 0.96, 'ACCESS-OM3-WOMBATlite model assessment', va='center', ha='center', fontsize=15, transform=fig.transFigure ) plt.subplots_adjust(left=0.1, bottom=0.05, top=0.95, right=0.9, hspace=0.35, wspace=0.25) # Colorbar for particle fluxes # ============================================= cbax8a = fig.add_axes([0.4, 0.23, 0.2, 0.01]) cbar8a = plt.colorbar(ps[7], cax=cbax8a, orientation='horizontal') cbar8a.ax.tick_params(labelbottom=False, bottom=False, labeltop=True, top=True) _ = plt.text(0.5, -0.54, 'Depth (m)', fontsize=fstic, ha='center', va='center', transform=cbax8a.transAxes)

In [47]:

mkmd.savefig(fig,
    "WOMBATlite global assessment metrics.",
    "Model metrics for WOMBATlite global assessment. "
    "[GitHub issue: Evaluation: WOMBATlite global metrics](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/47)",
    dpi=150
)

mkmd.savefig(fig, "WOMBATlite global assessment metrics.", "Model metrics for WOMBATlite global assessment. " "[GitHub issue: Evaluation: WOMBATlite global metrics](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/47)", dpi=150 )

Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global_02.png
Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global.md
Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global.md successfully.

Make Taylor Diagrams of the results¶

In [48]:

def get_southern_ocean(da):
    return da.sel(lat=slice(-90,-30))

def get_equatorial(da):
    return da.sel(lat=slice(-20,20), lon=slice(-180,-80))
    
def flatten(x):
    """Ensure dims are in the same order and flatten"""
    transpose_dims = ("time", "month", "lat", "lon", "province")
    if not isinstance(x, xr.DataArray):
        return np.asarray(x).ravel()
    if extra := set(x.dims) - set(transpose_dims):
        raise ValueError(f"{x.name}: unexpected dims {sorted(extra)}")
    return x.transpose(*(d for d in transpose_dims if d in x.dims)).values.ravel()

def get_southern_ocean(da): return da.sel(lat=slice(-90,-30)) def get_equatorial(da): return da.sel(lat=slice(-20,20), lon=slice(-180,-80)) def flatten(x): """Ensure dims are in the same order and flatten""" transpose_dims = ("time", "month", "lat", "lon", "province") if not isinstance(x, xr.DataArray): return np.asarray(x).ravel() if extra := set(x.dims) - set(transpose_dims): raise ValueError(f"{x.name}: unexpected dims {sorted(extra)}") return x.transpose(*(d for d in transpose_dims if d in x.dims)).values.ravel()

In [49]:

comparisons = {
    "Equatorial Surface NO$_3$ ($\mu$M)": (
        get_equatorial(np.log10(surno3_woa_reg)),
        get_equatorial(np.log10(surno3_om3_reg.clip(1e-3))),
    ),
    "Surface dFe (nM)": (
        surdfe_Huang_reg,
        surdfe_om3_reg
    ),
    "Surface chlorophyll\n(mg m$^{-3}$)": (
        surchl_cop_reg,
        surchl_om3_reg
    ),
    "Depth of the\nchlorophyll maximum (m)": (
        chmaxz_cop_reg,
        chmaxz_om3_reg
    ),
    "Depth-integrated NPP\n(mg C m$^{-2}$ day$^{-1}$)": (
        npp_modis_reg,
        npp_om3_reg
    ),
    "Southern Ocean air-sea\nCO$_2$ flux (mol C m$^{-2}$ year$^{-1}$)": (
        get_southern_ocean(-1*co2_cop_reg),
        get_southern_ocean(-1*co2_om3_reg),
    ),
    "Primary limiting nutrient": (
        lim_brown_reg,
        lim_om3_reg
    ),
    "Zooplankton grazing pressure (day$^{-1}$)": (
        zoograz_long,
        zoograz_long_om3
    ),
    "Sinking POC flux\n(mg C m$^{-2}$ day$^{-1}$))": (
        partflux["pocflux_obs"].values,
        partflux["pocflux_model"].values,
    ),
    'Sinking PIC flux\n(mg C m$^{-2}$ day$^{-1}$))': (
        partflux["picflux_obs"].values,
        partflux["picflux_model"].values
    ),
}

comparisons = {v: (flatten(o), flatten(m)) for v, (o, m) in comparisons.items()}

comparisons = { "Equatorial Surface NO$_3$ ($\mu$M)": ( get_equatorial(np.log10(surno3_woa_reg)), get_equatorial(np.log10(surno3_om3_reg.clip(1e-3))), ), "Surface dFe (nM)": ( surdfe_Huang_reg, surdfe_om3_reg ), "Surface chlorophyll\n(mg m$^{-3}$)": ( surchl_cop_reg, surchl_om3_reg ), "Depth of the\nchlorophyll maximum (m)": ( chmaxz_cop_reg, chmaxz_om3_reg ), "Depth-integrated NPP\n(mg C m$^{-2}$ day$^{-1}$)": ( npp_modis_reg, npp_om3_reg ), "Southern Ocean air-sea\nCO$_2$ flux (mol C m$^{-2}$ year$^{-1}$)": ( get_southern_ocean(-1*co2_cop_reg), get_southern_ocean(-1*co2_om3_reg), ), "Primary limiting nutrient": ( lim_brown_reg, lim_om3_reg ), "Zooplankton grazing pressure (day$^{-1}$)": ( zoograz_long, zoograz_long_om3 ), "Sinking POC flux\n(mg C m$^{-2}$ day$^{-1}$))": ( partflux["pocflux_obs"].values, partflux["pocflux_model"].values, ), 'Sinking PIC flux\n(mg C m$^{-2}$ day$^{-1}$))': ( partflux["picflux_obs"].values, partflux["picflux_model"].values ), } comparisons = {v: (flatten(o), flatten(m)) for v, (o, m) in comparisons.items()}

In [50]:

taylor_stats = {}

for var, data in comparisons.items():
    observed, predicted = data
    isfinite = np.isfinite(observed) & np.isfinite(predicted)
    observed = observed[isfinite]
    predicted = predicted[isfinite]

    # Calculate Taylor statistics
    stats = sm.taylor_statistics(predicted, observed)
    taylor_stats[var] = {k: np.array(v) for k, v in stats.items()}

    # Add bias
    taylor_stats[var]["bias"] = np.array([0.0, np.mean(predicted - observed)])

taylor_stats = {} for var, data in comparisons.items(): observed, predicted = data isfinite = np.isfinite(observed) & np.isfinite(predicted) observed = observed[isfinite] predicted = predicted[isfinite] # Calculate Taylor statistics stats = sm.taylor_statistics(predicted, observed) taylor_stats[var] = {k: np.array(v) for k, v in stats.items()} # Add bias taylor_stats[var]["bias"] = np.array([0.0, np.mean(predicted - observed)])

In [51]:

fig = plt.figure(figsize=(10,22))
gs = GridSpec(5,2)

for idx, (title, stats) in enumerate(taylor_stats.items()):
    ax = plt.subplot(gs[idx])

    vmin = 0 - np.abs(np.round(np.percentile(stats["bias"],0.01),1)*1.2)
    vmax = 0 + np.abs(np.round(np.percentile(stats["bias"],0.01),1)*1.2)
    vmin = 0 - np.abs(np.round(np.min(stats["bias"])))
    vmax = 0 + np.abs(np.round(np.min(stats["bias"])))

    sm.taylor_diagram(
        stats["sdev"],
        stats["crmsd"],
        stats["ccoef"],
        alpha=0.6,
        MarkerDisplayed="colorBar", cmapzdata=stats["bias"], locationColorBar="EastOutside",
        titleColorBar=" ", cmap=lighten(cmo.balance,0.8), cmap_marker='o', markerSize=15,
        cmap_vmax=vmax, cmap_vmin=vmin, 
        colRMS='grey', styleRMS='--', widthRMS=0.75, showlabelsRMS='on', titleRMS='off', rmsLabelFormat='0:.2f',
        colSTD='k', styleSTD=':', widthSTD=0.75, showlabelsSTD='off', titleSTD='off', 
        colCOR='k', styleCOR='-', widthCOR=0.5, showlabelsCOR='on', titleCOR='off', 
        markerOBS='o', colOBS='firebrick',
    )

    plt.text(
        0.5, 1.15, title,
        transform=ax.transAxes, va='bottom', ha='center', fontweight='bold'
    )
    plt.text(
        -0.05, 1.15, string.ascii_lowercase[idx],
        fontweight='bold', fontsize=fslab, ha='center', va='bottom', transform=ax.transAxes
    )

plt.subplots_adjust(hspace=0.5, wspace=0.1, right=0.95, left=0.05, bottom=0.05, top=0.95)

fig = plt.figure(figsize=(10,22)) gs = GridSpec(5,2) for idx, (title, stats) in enumerate(taylor_stats.items()): ax = plt.subplot(gs[idx]) vmin = 0 - np.abs(np.round(np.percentile(stats["bias"],0.01),1)*1.2) vmax = 0 + np.abs(np.round(np.percentile(stats["bias"],0.01),1)*1.2) vmin = 0 - np.abs(np.round(np.min(stats["bias"]))) vmax = 0 + np.abs(np.round(np.min(stats["bias"]))) sm.taylor_diagram( stats["sdev"], stats["crmsd"], stats["ccoef"], alpha=0.6, MarkerDisplayed="colorBar", cmapzdata=stats["bias"], locationColorBar="EastOutside", titleColorBar=" ", cmap=lighten(cmo.balance,0.8), cmap_marker='o', markerSize=15, cmap_vmax=vmax, cmap_vmin=vmin, colRMS='grey', styleRMS='--', widthRMS=0.75, showlabelsRMS='on', titleRMS='off', rmsLabelFormat='0:.2f', colSTD='k', styleSTD=':', widthSTD=0.75, showlabelsSTD='off', titleSTD='off', colCOR='k', styleCOR='-', widthCOR=0.5, showlabelsCOR='on', titleCOR='off', markerOBS='o', colOBS='firebrick', ) plt.text( 0.5, 1.15, title, transform=ax.transAxes, va='bottom', ha='center', fontweight='bold' ) plt.text( -0.05, 1.15, string.ascii_lowercase[idx], fontweight='bold', fontsize=fslab, ha='center', va='bottom', transform=ax.transAxes ) plt.subplots_adjust(hspace=0.5, wspace=0.1, right=0.95, left=0.05, bottom=0.05, top=0.95)

In [52]:

mkmd.savefig(fig,
    "WOMBATlite Taylor diagrams.",
    "Taylor diagrams for WOMBATlite global assessment. "
    "[GitHub issue: Evaluation: WOMBATlite global metrics](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/47)",
    dpi=150
)

mkmd.savefig(fig, "WOMBATlite Taylor diagrams.", "Taylor diagrams for WOMBATlite global assessment. " "[GitHub issue: Evaluation: WOMBATlite global metrics](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/47)", dpi=150 )

Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global_03.png
Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global.md
Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_iaf+wombatlite-test4-d28e0359/notebooks/mkfigs_output_MC_25km_jra_iaf+wombatlite-test4-d28e0359/mkmd/wombatlite_global.md successfully.

In [53]:

client.close()

client.close()

In [ ]: