Currents, barotropic streamfunction, SSH variability¶
Re-create Kiss et al 2020 figs 16, 17, 21, 22, 24
See
- https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/BarotropicStreamfunction/BarotropicStreamfunction.ipynb
- https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/Southwest_Atlantic/Southwest_atlantic_analysis_v3.ipynb
See
In [1]:
Copied!
# 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/25km-iaf-test-for-AK-expt-7df5ef4c/datastore.json"
# esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-5165c0f8/datastore.json"
# esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-00532b88/datastore.json"
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json"
#esm_file = "/g/data/zv30/non-cmip/ACCESS-CM3/cm3-run-03-06-2026/cm3-datastore/cm3-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/25km-iaf-test-for-AK-expt-7df5ef4c/datastore.json"
# esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-5165c0f8/datastore.json"
# esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-00532b88/datastore.json"
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test4-d28e0359/datastore.json"
#esm_file = "/g/data/zv30/non-cmip/ACCESS-CM3/cm3-run-03-06-2026/cm3-datastore/cm3-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]:
Copied!
# Parameters
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/experiment_datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/"
nbname = "Currents_streamfunction_variability.ipynb"
# Parameters
esm_file = "/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/experiment_datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/"
nbname = "Currents_streamfunction_variability.ipynb"
In [3]:
Copied!
import os
if not papermill:
import nci_ipynb # 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)
from exptdata_access import guess_experiment_from_esm_file
expt_key, info = guess_experiment_from_esm_file(esm_file)
model_name = info["model"] # OM3 or CM3
import os
if not papermill:
import nci_ipynb # 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)
from exptdata_access import guess_experiment_from_esm_file
expt_key, info = guess_experiment_from_esm_file(esm_file)
model_name = info["model"] # OM3 or CM3
In [4]:
Copied!
from model_agnostic import patch_broken_conda_env, patch_dask_workers
patch_broken_conda_env() # conda env workaround; no-op on healthy envs
import os
import xarray as xr
import gsw
import intake
from dask.distributed import Client
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from matplotlib.patches import Patch
import cmocean as cm
import cartopy.crs as ccrs
import cartopy.feature as cft
from shapely import geometry
import calendar
import cftime
import textwrap
import numpy as np
from tqdm.notebook import tqdm
import xgcm
from scipy.io import loadmat
from model_agnostic import patch_broken_conda_env, patch_dask_workers
patch_broken_conda_env() # conda env workaround; no-op on healthy envs
import os
import xarray as xr
import gsw
import intake
from dask.distributed import Client
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from matplotlib.patches import Patch
import cmocean as cm
import cartopy.crs as ccrs
import cartopy.feature as cft
from shapely import geometry
import calendar
import cftime
import textwrap
import numpy as np
from tqdm.notebook import tqdm
import xgcm
from scipy.io import loadmat
In [5]:
Copied!
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
In [6]:
Copied!
xr.set_options(keep_attrs=True);
xr.set_options(keep_attrs=True);
In [7]:
Copied!
client = Client(threads_per_worker=1)
patch_dask_workers(client) # patch workers too
print(client.dashboard_link)
client = Client(threads_per_worker=1)
patch_dask_workers(client) # patch workers too
print(client.dashboard_link)
http://127.0.0.1:8787/status
In [8]:
Copied!
client
client
Out[8]:
Client
Client-172cac0a-8672-11f1-a981-00000080fe80
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: http://127.0.0.1:8787/status |
Cluster Info
LocalCluster
10a84eeb
| Dashboard: http://127.0.0.1:8787/status | Workers: 16 |
| Total threads: 16 | Total memory: 188.56 GiB |
| Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-323e34d8-d3ae-46de-8eb8-ae746b4b5c0d
| Comm: tcp://127.0.0.1:34401 | Workers: 0 |
| Dashboard: http://127.0.0.1:8787/status | Total threads: 0 |
| Started: Just now | Total memory: 0 B |
Workers
Worker: 0
| Comm: tcp://127.0.0.1:34877 | Total threads: 1 |
| Dashboard: http://127.0.0.1:36621/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:32961 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-iigjl8l_ | |
Worker: 1
| Comm: tcp://127.0.0.1:33273 | Total threads: 1 |
| Dashboard: http://127.0.0.1:38179/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:37773 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-909k_o2k | |
Worker: 2
| Comm: tcp://127.0.0.1:32863 | Total threads: 1 |
| Dashboard: http://127.0.0.1:46145/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:43199 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-4orwlw6h | |
Worker: 3
| Comm: tcp://127.0.0.1:44179 | Total threads: 1 |
| Dashboard: http://127.0.0.1:40125/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:46419 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-0b3gve4p | |
Worker: 4
| Comm: tcp://127.0.0.1:33011 | Total threads: 1 |
| Dashboard: http://127.0.0.1:45329/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:40911 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-aevf5dxu | |
Worker: 5
| Comm: tcp://127.0.0.1:41361 | Total threads: 1 |
| Dashboard: http://127.0.0.1:33355/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:33497 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-5z037r0m | |
Worker: 6
| Comm: tcp://127.0.0.1:44225 | Total threads: 1 |
| Dashboard: http://127.0.0.1:38795/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:43233 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-69y9j5p2 | |
Worker: 7
| Comm: tcp://127.0.0.1:44679 | Total threads: 1 |
| Dashboard: http://127.0.0.1:34139/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:46531 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-l6_hfdde | |
Worker: 8
| Comm: tcp://127.0.0.1:44171 | Total threads: 1 |
| Dashboard: http://127.0.0.1:42865/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:44287 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-y53_7b5m | |
Worker: 9
| Comm: tcp://127.0.0.1:34729 | Total threads: 1 |
| Dashboard: http://127.0.0.1:44173/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:40663 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-09ynrm_v | |
Worker: 10
| Comm: tcp://127.0.0.1:40073 | Total threads: 1 |
| Dashboard: http://127.0.0.1:36161/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:34481 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-7osjpp4i | |
Worker: 11
| Comm: tcp://127.0.0.1:41639 | Total threads: 1 |
| Dashboard: http://127.0.0.1:33603/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:44573 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-sk8cjlkq | |
Worker: 12
| Comm: tcp://127.0.0.1:45859 | Total threads: 1 |
| Dashboard: http://127.0.0.1:41889/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:45163 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-q2czadlq | |
Worker: 13
| Comm: tcp://127.0.0.1:44139 | Total threads: 1 |
| Dashboard: http://127.0.0.1:39489/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:36375 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-y0q0enot | |
Worker: 14
| Comm: tcp://127.0.0.1:42881 | Total threads: 1 |
| Dashboard: http://127.0.0.1:40609/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:43581 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-_l9qid_p | |
Worker: 15
| Comm: tcp://127.0.0.1:45345 | Total threads: 1 |
| Dashboard: http://127.0.0.1:43355/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:42691 | |
| Local directory: /jobfs/174557261.gadi-pbs/dask-scratch-space/worker-bi7h5k7a | |
In [9]:
Copied!
IAF = esm_file.find('iaf') > 0
IAF
IAF = esm_file.find('iaf') > 0
IAF
Out[9]:
False
In [10]:
Copied!
om3exptname = os.path.basename(os.path.dirname(esm_file))
om3exptname
om3exptname = os.path.basename(os.path.dirname(esm_file))
om3exptname
Out[10]:
'MC_25km_jra_ryf+wombatlite-test3-f4d79e82'
Load ACCESS-OM3 data from ESM datastore¶
In [11]:
Copied!
exptname = os.path.basename(os.path.dirname(esm_file))
exptname = os.path.basename(os.path.dirname(esm_file))
In [12]:
Copied!
catalogs = [
# '/g/data/ol01/access-om3-output/access-om3-025/25km-iaf-test-for-AK-expt-7df5ef4c/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-5165c0f8/datastore.json',
esm_file,
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-COARE-12fab415/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm1-d968c801/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm2-5dc49da6/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm3-da330542/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm4-9fd08880/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm5-9b5dbfa9/datastore.json',
]
catalogs
catalogs = [
# '/g/data/ol01/access-om3-output/access-om3-025/25km-iaf-test-for-AK-expt-7df5ef4c/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-5165c0f8/datastore.json',
esm_file,
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-COARE-12fab415/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm1-d968c801/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm2-5dc49da6/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm3-da330542/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm4-9fd08880/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-gm5-9b5dbfa9/datastore.json',
]
catalogs
Out[12]:
['/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/experiment_datastore.json']
In [13]:
Copied!
om3datastores = { os.path.normpath(c).split(os.sep)[-2]:
intake.open_esm_datastore(c,
columns_with_iterables=[
"variable",
"variable_long_name",
"variable_standard_name",
"variable_cell_methods",
"variable_units"]
)
for c in catalogs }
om3datastores
om3datastores = { os.path.normpath(c).split(os.sep)[-2]:
intake.open_esm_datastore(c,
columns_with_iterables=[
"variable",
"variable_long_name",
"variable_standard_name",
"variable_cell_methods",
"variable_units"]
)
for c in catalogs }
om3datastores
Out[13]:
{'MC_25km_jra_ryf+wombatlite-test3-f4d79e82': <experiment_datastore catalog with 23 dataset(s) from 4495 asset(s)>}
In [14]:
Copied!
fields = { # key is OM3 name, value is OM2 name
"speed": None, # presumably includes high-frequency motion, so differs from sqrt(uo**2 + vo**2) calculated from monthly means
"uo": "u",
"vo": "v",
"umo_2d": None,
"vmo_2d": None,
"zos": None,
"zossq": None,
}
open_kwargs = dict(
# chunks={"xh": -1, "yh": -1},
# chunks={"time": -1},
decode_timedelta=True,
use_cftime=True
)
if model_name == "ACCESS-CM3":
open_kwargs["chunks"] = {"time": 1} # CM3 stores yearly files
data = {
vname: {expt: ds.search(variable=vname, frequency="1mon").to_dask(
xarray_open_kwargs = open_kwargs
)[vname]
for expt, ds in om3datastores.items()
}
for vname in fields
}
# uo/vo are only plotted at 2 depths; select them now (before the time mean)
# so only those levels are read, not all 75
for vname in ["uo", "vo"]:
for expt in data[vname]:
data[vname][expt] = data[vname][expt].sel(z_l=[0, 1000], method="nearest")
# omit latitudes with grid bug in this run https://github.com/ACCESS-NRI/ocean_model_grid_generator/issues/7
for vname, d in data.items():
try:
d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yh=slice(10, None))
except ValueError:
d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yq=slice(10, None))
except KeyError:
pass
if IAF:
for vname, d in data.items():
for expt in d:
try:
d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True)
except KeyError:
pass
fields = { # key is OM3 name, value is OM2 name
"speed": None, # presumably includes high-frequency motion, so differs from sqrt(uo**2 + vo**2) calculated from monthly means
"uo": "u",
"vo": "v",
"umo_2d": None,
"vmo_2d": None,
"zos": None,
"zossq": None,
}
open_kwargs = dict(
# chunks={"xh": -1, "yh": -1},
# chunks={"time": -1},
decode_timedelta=True,
use_cftime=True
)
if model_name == "ACCESS-CM3":
open_kwargs["chunks"] = {"time": 1} # CM3 stores yearly files
data = {
vname: {expt: ds.search(variable=vname, frequency="1mon").to_dask(
xarray_open_kwargs = open_kwargs
)[vname]
for expt, ds in om3datastores.items()
}
for vname in fields
}
# uo/vo are only plotted at 2 depths; select them now (before the time mean)
# so only those levels are read, not all 75
for vname in ["uo", "vo"]:
for expt in data[vname]:
data[vname][expt] = data[vname][expt].sel(z_l=[0, 1000], method="nearest")
# omit latitudes with grid bug in this run https://github.com/ACCESS-NRI/ocean_model_grid_generator/issues/7
for vname, d in data.items():
try:
d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yh=slice(10, None))
except ValueError:
d['25km-iaf-test-for-AK-expt-7df5ef4c'] = d['25km-iaf-test-for-AK-expt-7df5ef4c'].isel(yq=slice(10, None))
except KeyError:
pass
if IAF:
for vname, d in data.items():
for expt in d:
try:
d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True)
except KeyError:
pass
Trim to a consistent time range¶
In [15]:
Copied!
# find largest overlapping time range
# earliest end date in data
datastop = np.min([ data['uo'][e].time.values[-1] for e in data['uo'].keys() ])
# latest start date in data
datastart = np.max([ data['uo'][e].time.values[0] for e in data['uo'].keys() ])
print('common data range:', datastart, ' to ', datastop)
nyears = 30
datastart = np.max([
datastart,
cftime.datetime(cftime.to_tuple(datastop)[0]-nyears+1, 1, 1, calendar=datastop.calendar)
])
# bracket to start/end of year
datestart = cftime.datetime(cftime.to_tuple(datastart)[0], 1, 1, calendar=datastop.calendar)
datestop = cftime.datetime(cftime.to_tuple(datastop)[0], 12, 31, 23, 59, 59, calendar=datastop.calendar)
timeslice = slice(datestart, datestop)
timeslice
# find largest overlapping time range
# earliest end date in data
datastop = np.min([ data['uo'][e].time.values[-1] for e in data['uo'].keys() ])
# latest start date in data
datastart = np.max([ data['uo'][e].time.values[0] for e in data['uo'].keys() ])
print('common data range:', datastart, ' to ', datastop)
nyears = 30
datastart = np.max([
datastart,
cftime.datetime(cftime.to_tuple(datastop)[0]-nyears+1, 1, 1, calendar=datastop.calendar)
])
# bracket to start/end of year
datestart = cftime.datetime(cftime.to_tuple(datastart)[0], 1, 1, calendar=datastop.calendar)
datestop = cftime.datetime(cftime.to_tuple(datastop)[0], 12, 31, 23, 59, 59, calendar=datastop.calendar)
timeslice = slice(datestart, datestop)
timeslice
common data range: 1900-01-16 12:00:00 to 1930-12-16 12:00:00
Out[15]:
slice(cftime.datetime(1901, 1, 1, 0, 0, 0, 0, calendar='noleap', has_year_zero=True), cftime.datetime(1930, 12, 31, 23, 59, 59, 0, calendar='noleap', has_year_zero=True), None)
In [16]:
Copied!
# trim everything to timeslice and time-average
data = { vname: {expt: da.sel(time=timeslice).mean('time') for expt, da in dct.items()} for vname, dct in data.items() }
# trim everything to timeslice and time-average
data = { vname: {expt: da.sel(time=timeslice).mean('time') for expt, da in dct.items()} for vname, dct in data.items() }
Load OM3 grid data¶
In [17]:
Copied!
# regex matches access-om3 and access-cm3 static files
static = xr.open_dataset(om3datastores[exptname].search(filename=".*mom6.static.nc").df.loc[0, "path"]).drop_vars("time")
# regex matches access-om3 and access-cm3 static files
static = xr.open_dataset(om3datastores[exptname].search(filename=".*mom6.static.nc").df.loc[0, "path"]).drop_vars("time")
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:38: UserWarning: An error occurred while loading the ERDDAP data fetcher, it will not be available ! <class 'ImportError'> cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py) warnings.warn( /g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:50: UserWarning: An error occurred while loading the ArgoVis data fetcher, it will not be available ! <class 'ImportError'> cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py) warnings.warn( /g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:66: UserWarning: An error occurred while loading the GDAC data fetcher, it will not be available ! <class 'ImportError'> cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py) warnings.warn( /g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:91: UserWarning: An error occurred while loading the ERDDAP index fetcher, it will not be available ! <class 'ImportError'> cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py) warnings.warn( /g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/argopy/utils/lists.py:107: UserWarning: An error occurred while loading the GDAC index fetcher, it will not be available ! <class 'ImportError'> cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py) warnings.warn( /g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'argo' loading failed: cannot import name '_quote_string_constraints' from 'erddapy.erddapy' (/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/erddapy/erddapy.py) external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
In [18]:
Copied!
# See https://xgcm.readthedocs.io/en/latest/xgcm-examples/03_MOM6/#xgcm-grid-definition
symmetric = True # True for ACCESS-OM3 2025.08.000 and later
if symmetric: # ACCESS-OM3 version 2025.08.000 and later: one more q point than h points
coords={
'X': {'center': 'xh', 'outer': 'xq'},
'Y': {'center': 'yh', 'outer': 'yq'},
# 'Z': { 'inner': 'z_l', 'outer': 'z_i'}
}
else: # ACCESS-OM3 versions prior to 2025.08.000: equal number of q and h points
coords={
'X': {'center': 'xh', 'right': 'xq'},
'Y': {'center': 'yh', 'right': 'yq'},
# 'Z': { 'inner': 'z_l', 'outer': 'z_i'}
}
grid = xgcm.Grid(
static,
coords=coords,
metrics = {
('X',): ['dxt', 'dxCu', 'dxCv'], # X distances
('Y',): ['dyt', 'dyCu', 'dyCv'], # Y distances
('X', 'Y'): ['areacello', 'areacello_cu', 'areacello_cv', 'areacello_bu'] # Areas
},
periodic=['X'],
autoparse_metadata=False
# Note, autoparse_metadata=False is required below for xgcm >= 0.9.0 to prevent xgcm
# from trying to autoparse non-existent grid information from the file metadata,
# leading to the error:
# ValueError: Found two coordinates without `c_grid_axis_shift` attribute for axis X
)
# See https://xgcm.readthedocs.io/en/latest/xgcm-examples/03_MOM6/#xgcm-grid-definition
symmetric = True # True for ACCESS-OM3 2025.08.000 and later
if symmetric: # ACCESS-OM3 version 2025.08.000 and later: one more q point than h points
coords={
'X': {'center': 'xh', 'outer': 'xq'},
'Y': {'center': 'yh', 'outer': 'yq'},
# 'Z': { 'inner': 'z_l', 'outer': 'z_i'}
}
else: # ACCESS-OM3 versions prior to 2025.08.000: equal number of q and h points
coords={
'X': {'center': 'xh', 'right': 'xq'},
'Y': {'center': 'yh', 'right': 'yq'},
# 'Z': { 'inner': 'z_l', 'outer': 'z_i'}
}
grid = xgcm.Grid(
static,
coords=coords,
metrics = {
('X',): ['dxt', 'dxCu', 'dxCv'], # X distances
('Y',): ['dyt', 'dyCu', 'dyCv'], # Y distances
('X', 'Y'): ['areacello', 'areacello_cu', 'areacello_cv', 'areacello_bu'] # Areas
},
periodic=['X'],
autoparse_metadata=False
# Note, autoparse_metadata=False is required below for xgcm >= 0.9.0 to prevent xgcm
# from trying to autoparse non-existent grid information from the file metadata,
# leading to the error:
# ValueError: Found two coordinates without `c_grid_axis_shift` attribute for axis X
)
Load observations¶
In [19]:
Copied!
def xtile(ds_in, x="longitude"):
"""
Tile dataset in longitude to cover MOM6 range.
ds_in: input dataset, assumed 2d with x coordinate second
x: name of x coordinate, defaults to "longitude"
"""
x_tiled = np.tile(ds_in[x], 2)
x_tiled[:ds_in[x].size] -= 360
ds_out = xr.Dataset(
attrs=ds_in.attrs,
coords=dict(ds_in.coords) | {x: x_tiled},
data_vars={k: (v.dims, np.tile(v, (1,2))) for k, v in ds_in.items()}
)
for k, v in ds_out.coords.items():
v.attrs = ds_in[k].attrs
for k, v in ds_out.items():
v.attrs = ds_in[k].attrs
return ds_out
def xtile(ds_in, x="longitude"):
"""
Tile dataset in longitude to cover MOM6 range.
ds_in: input dataset, assumed 2d with x coordinate second
x: name of x coordinate, defaults to "longitude"
"""
x_tiled = np.tile(ds_in[x], 2)
x_tiled[:ds_in[x].size] -= 360
ds_out = xr.Dataset(
attrs=ds_in.attrs,
coords=dict(ds_in.coords) | {x: x_tiled},
data_vars={k: (v.dims, np.tile(v, (1,2))) for k, v in ds_in.items()}
)
for k, v in ds_out.coords.items():
v.attrs = ds_in[k].attrs
for k, v in ds_out.items():
v.attrs = ds_in[k].attrs
return ds_out
Barotropic streamfunction¶
Observational estimate from Colin de Verdiere and Ollitrault (2016).
In [20]:
Copied!
psi_data = loadmat('/g/data/ik11/from_hh5_tmp/cosima/observations_new/original/psi_barotropic.mat')
psi_ds = xtile(xr.Dataset(
coords={"latitude": psi_data["lat_p"][0, :],
"longitude": psi_data["lon_p"][0, :]},
data_vars={"psi": (("latitude", "longitude"), psi_data["psi"]/1e6)} # convert to Sv
))
psi_ds["psi"].attrs={'units': 'Sv', 'long_name': 'Barotropic streamfunction'}
psi_ds.attrs = {'author': 'Observational estimate from Colin de Verdiere and Ollitrault (2016) https://doi.org/10.1175/jpo-d-15-0046.1'}
psi_data = loadmat('/g/data/ik11/from_hh5_tmp/cosima/observations_new/original/psi_barotropic.mat')
psi_ds = xtile(xr.Dataset(
coords={"latitude": psi_data["lat_p"][0, :],
"longitude": psi_data["lon_p"][0, :]},
data_vars={"psi": (("latitude", "longitude"), psi_data["psi"]/1e6)} # convert to Sv
))
psi_ds["psi"].attrs={'units': 'Sv', 'long_name': 'Barotropic streamfunction'}
psi_ds.attrs = {'author': 'Observational estimate from Colin de Verdiere and Ollitrault (2016) https://doi.org/10.1175/jpo-d-15-0046.1'}
In [21]:
Copied!
psi_ds
psi_ds
Out[21]:
<xarray.Dataset> Size: 911kB
Dimensions: (latitude: 157, longitude: 720)
Coordinates:
* latitude (latitude) float64 1kB -75.5 -74.5 -73.5 -72.5 ... 78.5 79.5 80.5
* longitude (longitude) float64 6kB -539.5 -538.5 -537.5 ... 178.5 179.5
Data variables:
psi (latitude, longitude) float64 904kB nan nan nan ... nan nan nan
Attributes:
author: Observational estimate from Colin de Verdiere and Ollitrault (2...In [22]:
Copied!
psi_ds["psi"].plot()
psi_ds["psi"].plot()
Out[22]:
<matplotlib.collections.QuadMesh at 0x148e2b4d5070>
ANDRO deep currents¶
Ollitrault, M. and Rannou, J.-P. (2013). ANDRO: An Argo-based deep displacement dataset. Journal of Atmospheric and Oceanic Technology, 30(4):759–788. http://dx.doi.org/10.1175/JTECH-D-12-00073.1
2026 release: https://doi.org/10.17882/47077 and https://www.seanoe.org/data/00360/47077/
BUG: data seems misaligned with land in places
In [23]:
Copied!
andro_ds = xtile(
xr.open_dataset(
'/g/data/ik11/observations/ANDRO_Argo_displacements/2026/127689.nc'
).transpose('latitude', 'longitude'))
# calculate speed
andro_ds["speed"] = np.sqrt( andro_ds['mean_u']**2 + andro_ds['mean_v']**2 )/100 # convert cm/s to m/s
andro_ds["speed"].attrs = andro_ds["mean_u"].attrs | {'units': 'm s-1', 'long_name': 'Speed at 1000m'}
del andro_ds["speed"].attrs["Units"]
del andro_ds["speed"].attrs["Name"]
andro_ds = xtile(
xr.open_dataset(
'/g/data/ik11/observations/ANDRO_Argo_displacements/2026/127689.nc'
).transpose('latitude', 'longitude'))
# calculate speed
andro_ds["speed"] = np.sqrt( andro_ds['mean_u']**2 + andro_ds['mean_v']**2 )/100 # convert cm/s to m/s
andro_ds["speed"].attrs = andro_ds["mean_u"].attrs | {'units': 'm s-1', 'long_name': 'Speed at 1000m'}
del andro_ds["speed"].attrs["Units"]
del andro_ds["speed"].attrs["Name"]
In [24]:
Copied!
andro_ds["speed"]
andro_ds["speed"]
Out[24]:
<xarray.DataArray 'speed' (latitude: 545, longitude: 1440)> Size: 6MB
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]], shape=(545, 1440))
Coordinates:
* latitude (latitude) float64 4kB -77.01 -76.9 -76.78 ... 89.49 89.69 89.9
* longitude (longitude) float64 12kB -540.0 -539.5 -539.0 ... 179.0 179.5
Attributes:
units: m s-1
long_name: Speed at 1000mIn [25]:
Copied!
andro_ds["speed"].plot()
andro_ds["speed"].plot()
Out[25]:
<matplotlib.collections.QuadMesh at 0x148e202bbef0>
SSH variance¶
Drifters¶
Drifter surface speed data from Kiss et al fig 24 was from Laurindo et al 2017, downloaded in 2017.
Here we use an updated 2025 version from http://www.aoml.noaa.gov/phod/dac/dac_meanvel.php instead.
In [26]:
Copied!
drifters_ds = xtile(
xr.open_dataset(
'/g/data/ik11/observations/Laurindo_et_al_surface_current/v3.11/drifter_annualmeans.nc'
).transpose('latitude', 'longitude'))
# calculate speed
drifters_ds["speed"] = np.sqrt( drifters_ds['U']**2 + drifters_ds['V']**2 ).compute()
drifters_ds["speed"].attrs = drifters_ds["U"].attrs | {'units': 'm s-1', 'long_name': 'Surface speed'}
drifters_ds = xtile(
xr.open_dataset(
'/g/data/ik11/observations/Laurindo_et_al_surface_current/v3.11/drifter_annualmeans.nc'
).transpose('latitude', 'longitude'))
# calculate speed
drifters_ds["speed"] = np.sqrt( drifters_ds['U']**2 + drifters_ds['V']**2 ).compute()
drifters_ds["speed"].attrs = drifters_ds["U"].attrs | {'units': 'm s-1', 'long_name': 'Surface speed'}
In [27]:
Copied!
drifters_ds["speed"]
drifters_ds["speed"]
Out[27]:
<xarray.DataArray 'speed' (latitude: 632, longitude: 2880)> Size: 15MB
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]], shape=(632, 2880))
Coordinates:
* latitude (latitude) float64 5kB -72.88 -72.62 -72.38 ... 84.38 84.62 84.88
* longitude (longitude) float64 23kB -539.9 -539.6 -539.4 ... 179.6 179.9
Attributes:
colorBarMaximum: 1.5
colorBarMinimum: -1.5
long_name: Surface speed
units: m s-1
valid_max: 1.2255
valid_min: -1.2494In [28]:
Copied!
drifters_ds["speed"].plot()
drifters_ds["speed"].plot()
Out[28]:
<matplotlib.collections.QuadMesh at 0x148e20373200>
In [29]:
Copied!
regions = { # [minx, maxx, miny, maxy], using model longitude range (-280 to 80)
"Agulhas region": [0, 80, -60, -15], # fig 16 in Kiss et al 2020
"EAC region": [-220, -160, -52, -20], # fig 17 in Kiss et al 2020
"Kuroshio region": [-240, -120, 10, 65], # fig 21 in Kiss et al 2020
"Gulf Stream region": [-95, -5, 10, 65], # fig 22 in Kiss et al 2020
"Argentine Basin": [-60, -30, -60, -30], # fig 24 in Kiss et al 2020
# "Global": [-280, 80, -90, 90],
# "Arctic": [-280, 80, 65, 90],
# "Southern Ocean": [-280, 80, -82, -63],
# "ACC": [-280, 80, -63, -45],
# "Southern Pacific": [-210, -70, -45, -20],
# "Tropical Pacific": [-240, -100, -20, 20],
# "North Pacific": [-240, -100, 20, 65],
# "South Atlantic": [-60, 20, -45, -20],
# "Tropical Atlantic": [-70, 20, -20, 20],
# "North Atlantic": [-100, 0, 20, 65],
# "Indian": [30, 120, -45, 20],
# "Aegean Sea": [18, 27.5, 34, 44],
# "Black Sea": [27.5, 43, 40.5, 48],
# "Baltic Sea": [13, 30, 53, 58],
# "Mediterranean Sea": [0, 35, 31, 41],
# "Red Sea": [33, 44, 12, 29],
# "Persian Gulf": [47, 56, 24, 31],
# "White Sea": [31, 41, 63, 68],
}
regions = {k: dict(zip(["minx", "maxx", "miny", "maxy"], v)) for k, v in regions.items()} # convert to dicts
for r, d in regions.items():
for k, x in d.items():
if k in ["minx", "maxx"] and x != max(-280, min(x, 80)):
raise ValueError(f"{r} {k} = {x} is outside the range -280 to 80")
fig = plt.figure(figsize=(10, 4))
colors = mpl.color_sequences['tab20']
try: # set Global color to white
colors[list(regions.keys()).index("Global")] = (1, 1, 1)
except ValueError:
pass
ax = plt.axes(position=[0.05,0.05,0.7,0.9], projection=ccrs.Robinson())
ax.coastlines(resolution="110m")
ax.gridlines(draw_labels=False)
legend_elements = []
for i, (region, limits) in enumerate(regions.items()):
ax.add_geometries([geometry.box(**limits)], crs=ccrs.PlateCarree(), color=colors[i], alpha=0.3)
legend_elements.append(Patch(color=colors[i], alpha=0.3, label=region))
ax.legend(handles=legend_elements, bbox_to_anchor=(1.2, 1.00))
plt.title("Regions (BUG: incorrect shapes north of 65°N)")
regions = { # [minx, maxx, miny, maxy], using model longitude range (-280 to 80)
"Agulhas region": [0, 80, -60, -15], # fig 16 in Kiss et al 2020
"EAC region": [-220, -160, -52, -20], # fig 17 in Kiss et al 2020
"Kuroshio region": [-240, -120, 10, 65], # fig 21 in Kiss et al 2020
"Gulf Stream region": [-95, -5, 10, 65], # fig 22 in Kiss et al 2020
"Argentine Basin": [-60, -30, -60, -30], # fig 24 in Kiss et al 2020
# "Global": [-280, 80, -90, 90],
# "Arctic": [-280, 80, 65, 90],
# "Southern Ocean": [-280, 80, -82, -63],
# "ACC": [-280, 80, -63, -45],
# "Southern Pacific": [-210, -70, -45, -20],
# "Tropical Pacific": [-240, -100, -20, 20],
# "North Pacific": [-240, -100, 20, 65],
# "South Atlantic": [-60, 20, -45, -20],
# "Tropical Atlantic": [-70, 20, -20, 20],
# "North Atlantic": [-100, 0, 20, 65],
# "Indian": [30, 120, -45, 20],
# "Aegean Sea": [18, 27.5, 34, 44],
# "Black Sea": [27.5, 43, 40.5, 48],
# "Baltic Sea": [13, 30, 53, 58],
# "Mediterranean Sea": [0, 35, 31, 41],
# "Red Sea": [33, 44, 12, 29],
# "Persian Gulf": [47, 56, 24, 31],
# "White Sea": [31, 41, 63, 68],
}
regions = {k: dict(zip(["minx", "maxx", "miny", "maxy"], v)) for k, v in regions.items()} # convert to dicts
for r, d in regions.items():
for k, x in d.items():
if k in ["minx", "maxx"] and x != max(-280, min(x, 80)):
raise ValueError(f"{r} {k} = {x} is outside the range -280 to 80")
fig = plt.figure(figsize=(10, 4))
colors = mpl.color_sequences['tab20']
try: # set Global color to white
colors[list(regions.keys()).index("Global")] = (1, 1, 1)
except ValueError:
pass
ax = plt.axes(position=[0.05,0.05,0.7,0.9], projection=ccrs.Robinson())
ax.coastlines(resolution="110m")
ax.gridlines(draw_labels=False)
legend_elements = []
for i, (region, limits) in enumerate(regions.items()):
ax.add_geometries([geometry.box(**limits)], crs=ccrs.PlateCarree(), color=colors[i], alpha=0.3)
legend_elements.append(Patch(color=colors[i], alpha=0.3, label=region))
ax.legend(handles=legend_elements, bbox_to_anchor=(1.2, 1.00))
plt.title("Regions (BUG: incorrect shapes north of 65°N)")
Out[29]:
Text(0.5, 1.0, 'Regions (BUG: incorrect shapes north of 65°N)')
Calculate speed at particular depths¶
Interpolate onto NE corner points¶
In [30]:
Copied!
# new dict containing interpolated time-mean dataarrays
interp_axis = {"uo": "Y", "vo": "X"} # interpolate to NE corner to match B-grid OM2 data
coord = {"uo": "yh", "vo": "xh"}
for vname in interp_axis:
data[vname+"_ne"] = dict()
for expt, da in data[vname].items():
v = grid.interp(da.chunk(chunks={coord[vname]: -1}),
interp_axis[vname], boundary="extend")
v = v.assign_coords({
"z_l": da.z_l,
"geolon_c": static.geolon_c,
"geolat_c": static.geolat_c})
data[vname+"_ne"][expt] = v
# new dict containing interpolated time-mean dataarrays
interp_axis = {"uo": "Y", "vo": "X"} # interpolate to NE corner to match B-grid OM2 data
coord = {"uo": "yh", "vo": "xh"}
for vname in interp_axis:
data[vname+"_ne"] = dict()
for expt, da in data[vname].items():
v = grid.interp(da.chunk(chunks={coord[vname]: -1}),
interp_axis[vname], boundary="extend")
v = v.assign_coords({
"z_l": da.z_l,
"geolon_c": static.geolon_c,
"geolat_c": static.geolat_c})
data[vname+"_ne"][expt] = v
In [31]:
Copied!
# calculate speed at given depths, using uo, vo interpolated to NE corner of cells (B-grid vel point)
import dask
depths = [0, 1000]
for depth in depths:
data[f"speed_{depth}"] = dict()
for expt in om3datastores:
for depth in depths:
u = data['uo_ne'][expt].sel(z_l=depth, method='nearest')
v = data['vo_ne'][expt].sel(z_l=depth, method='nearest')
data[f"speed_{depth}"][expt] = np.sqrt(u**2 + v**2 )
# compute all speeds in ONE pass (shares the uo/vo time-mean graph),
# otherwise every region plot re-reads the full uo/vo record
(computed,) = dask.compute({f"speed_{d}": dict(data[f"speed_{d}"]) for d in depths})
for d in depths:
data[f"speed_{d}"] = computed[f"speed_{d}"]
# calculate speed at given depths, using uo, vo interpolated to NE corner of cells (B-grid vel point)
import dask
depths = [0, 1000]
for depth in depths:
data[f"speed_{depth}"] = dict()
for expt in om3datastores:
for depth in depths:
u = data['uo_ne'][expt].sel(z_l=depth, method='nearest')
v = data['vo_ne'][expt].sel(z_l=depth, method='nearest')
data[f"speed_{depth}"][expt] = np.sqrt(u**2 + v**2 )
# compute all speeds in ONE pass (shares the uo/vo time-mean graph),
# otherwise every region plot re-reads the full uo/vo record
(computed,) = dask.compute({f"speed_{d}": dict(data[f"speed_{d}"]) for d in depths})
for d in depths:
data[f"speed_{d}"] = computed[f"speed_{d}"]
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.07/lib/python3.12/site-packages/distributed/client.py:3387: UserWarning: Sending large graph of size 51.87 MiB. This may cause some slowdown. Consider loading the data with Dask directly or using futures or delayed objects to embed the data into the graph without repetition. See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information. warnings.warn(
Plot¶
In [32]:
Copied!
blue_marble = plt.imread('/g/data/ik11/grids/BlueMarble.tiff')
blue_marble_extent = (-180, 180, -90, 90)
blue_marble = plt.imread('/g/data/ik11/grids/BlueMarble.tiff')
blue_marble_extent = (-180, 180, -90, 90)
In [33]:
Copied!
# plot surface current speed
margin = 10
cmp = plt.get_cmap('nipy_spectral')
clev1 = np.arange(0, 1.0, 0.01)
for region, bounds in regions.items():
# if region == "Kuroshio region":
# break
# if region != "Argentine Basin":
# continue
minx, maxx, miny, maxy = bounds.values()
print(region, minx, maxx, miny, maxy)
xslice = slice(minx-margin, maxx+margin)
yslice = slice(miny-margin, maxy+margin)
nrows = 1 # TODO
ncols = len(data["speed_0"]) + 1
proj = ccrs.PlateCarree(central_longitude=(minx+maxx)/2)
fig, axes = plt.subplots(nrows, ncols, figsize=(6*ncols, 3*nrows),
subplot_kw=dict(projection=proj,
frame_on=False))
# plot model
for i, (expt, da) in enumerate(data["speed_0"].items(), start=1):
ax = plt.subplot(nrows, ncols, i, projection=proj)
d = da.sel(xq=xslice, yq=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title(f"{expt}\n{datestart.strftime('%Y-%m')} - {datestop.strftime('%Y-%m')} mean")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
if i > 1:
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
# plot obs
ax = plt.subplot(nrows, ncols, ncols, projection=proj)
d = drifters_ds["speed"].sel(longitude=xslice, latitude=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title("Drifter observations")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
ax1 = plt.axes([0.92, 0.2, 0.01, 0.6])
cb = plt.colorbar(p1, cax=ax1, orientation='vertical', pad=.05)
cb.ax.set_ylabel('Surface current speed (m/s)')
mkmd.savefig(fig, f"Current Speed: {region}", "Model vs drifter current speed (surface). Region: "+region+" GitHub issues: [#15](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/15) Regional SLA & barotropic streamfunction, [#18](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/18) SW Atlantic, [#41](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/41) Western boundary current timeseries", dpi=150)
# plot surface current speed
margin = 10
cmp = plt.get_cmap('nipy_spectral')
clev1 = np.arange(0, 1.0, 0.01)
for region, bounds in regions.items():
# if region == "Kuroshio region":
# break
# if region != "Argentine Basin":
# continue
minx, maxx, miny, maxy = bounds.values()
print(region, minx, maxx, miny, maxy)
xslice = slice(minx-margin, maxx+margin)
yslice = slice(miny-margin, maxy+margin)
nrows = 1 # TODO
ncols = len(data["speed_0"]) + 1
proj = ccrs.PlateCarree(central_longitude=(minx+maxx)/2)
fig, axes = plt.subplots(nrows, ncols, figsize=(6*ncols, 3*nrows),
subplot_kw=dict(projection=proj,
frame_on=False))
# plot model
for i, (expt, da) in enumerate(data["speed_0"].items(), start=1):
ax = plt.subplot(nrows, ncols, i, projection=proj)
d = da.sel(xq=xslice, yq=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title(f"{expt}\n{datestart.strftime('%Y-%m')} - {datestop.strftime('%Y-%m')} mean")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
if i > 1:
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
# plot obs
ax = plt.subplot(nrows, ncols, ncols, projection=proj)
d = drifters_ds["speed"].sel(longitude=xslice, latitude=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title("Drifter observations")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
ax1 = plt.axes([0.92, 0.2, 0.01, 0.6])
cb = plt.colorbar(p1, cax=ax1, orientation='vertical', pad=.05)
cb.ax.set_ylabel('Surface current speed (m/s)')
mkmd.savefig(fig, f"Current Speed: {region}", "Model vs drifter current speed (surface). Region: "+region+" GitHub issues: [#15](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/15) Regional SLA & barotropic streamfunction, [#18](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/18) SW Atlantic, [#41](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/41) Western boundary current timeseries", dpi=150)
Agulhas region 0 80 -60 -15
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_01.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md
Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. EAC region -220 -160 -52 -20
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_02.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. Kuroshio region -240 -120 10 65
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_03.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. Gulf Stream region -95 -5 10 65
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_04.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. Argentine Basin -60 -30 -60 -30
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_05.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully.
In [34]:
Copied!
# plot 1000m current speed
margin = 1
cmp = plt.get_cmap('nipy_spectral')
clev1 = np.arange(0, 1.0, 0.01) * 0.4
for region, bounds in regions.items():
# if region != "Argentine Basin":
# continue
minx, maxx, miny, maxy = bounds.values()
print(region, minx, maxx, miny, maxy)
xslice = slice(minx-margin, maxx+margin)
yslice = slice(miny-margin, maxy+margin)
nrows = 1 # TODO
ncols = len(data["speed_0"]) + 1
proj = ccrs.PlateCarree(central_longitude=(minx+maxx)/2)
fig, axes = plt.subplots(nrows, ncols, figsize=(6*ncols, 3*nrows),
subplot_kw=dict(projection=proj,
frame_on=False))
# plot model
for i, (expt, da) in enumerate(data["speed_1000"].items(), start=1):
ax = plt.subplot(nrows, ncols, i, projection=proj)
d = da.sel(xq=xslice, yq=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title(f"{expt}\n{datestart.strftime('%Y-%m')} - {datestop.strftime('%Y-%m')} mean")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
if i > 1:
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
# plot obs
ax = plt.subplot(nrows, ncols, ncols, projection=proj)
d = andro_ds["speed"].sel(longitude=xslice, latitude=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title("ANDRO Argo observations")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
ax1 = plt.axes([0.92, 0.2, 0.01, 0.6])
cb = plt.colorbar(p1, cax=ax1, orientation='vertical', pad=.05)
cb.ax.set_ylabel('Current speed at 1000m (m/s)')
mkmd.savefig(fig, f"Current Speed: {region}", "Model vs drifter current speed (1000m). Region: "+region+" GitHub issues: [#15](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/15) Regional SLA & barotropic streamfunction, [#18](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/18) SW Atlantic, [#41](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/41) Western boundary current timeseries", dpi=150)
# plot 1000m current speed
margin = 1
cmp = plt.get_cmap('nipy_spectral')
clev1 = np.arange(0, 1.0, 0.01) * 0.4
for region, bounds in regions.items():
# if region != "Argentine Basin":
# continue
minx, maxx, miny, maxy = bounds.values()
print(region, minx, maxx, miny, maxy)
xslice = slice(minx-margin, maxx+margin)
yslice = slice(miny-margin, maxy+margin)
nrows = 1 # TODO
ncols = len(data["speed_0"]) + 1
proj = ccrs.PlateCarree(central_longitude=(minx+maxx)/2)
fig, axes = plt.subplots(nrows, ncols, figsize=(6*ncols, 3*nrows),
subplot_kw=dict(projection=proj,
frame_on=False))
# plot model
for i, (expt, da) in enumerate(data["speed_1000"].items(), start=1):
ax = plt.subplot(nrows, ncols, i, projection=proj)
d = da.sel(xq=xslice, yq=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title(f"{expt}\n{datestart.strftime('%Y-%m')} - {datestop.strftime('%Y-%m')} mean")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
if i > 1:
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
# plot obs
ax = plt.subplot(nrows, ncols, ncols, projection=proj)
d = andro_ds["speed"].sel(longitude=xslice, latitude=yslice)
p1 = d.plot.contourf(cmap=cmp,
levels=clev1,
add_colorbar=False,
extend='max',
transform=ccrs.PlateCarree(),
)
p1.cmap.set_over(color=cmp(255), alpha=None)
plt.title("ANDRO Argo observations")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='gray', alpha=0.5, linestyle='-')
gl.top_labels = False
gl.right_labels = False
gl.left_labels = False
gl.xlocator = mticker.FixedLocator(np.arange(-270, 180, 10))
plt.imshow(
blue_marble,
extent=blue_marble_extent,
transform=ccrs.PlateCarree(),
origin="upper")
ax.coastlines(resolution="50m")
ax.set_extent(bounds.values())
ax1 = plt.axes([0.92, 0.2, 0.01, 0.6])
cb = plt.colorbar(p1, cax=ax1, orientation='vertical', pad=.05)
cb.ax.set_ylabel('Current speed at 1000m (m/s)')
mkmd.savefig(fig, f"Current Speed: {region}", "Model vs drifter current speed (1000m). Region: "+region+" GitHub issues: [#15](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/15) Regional SLA & barotropic streamfunction, [#18](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/18) SW Atlantic, [#41](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/41) Western boundary current timeseries", dpi=150)
Agulhas region 0 80 -60 -15
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_06.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md This title already exists in the notebook markdown – appending an additional figure. Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. EAC region -220 -160 -52 -20
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_07.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md This title already exists in the notebook markdown – appending an additional figure. Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. Kuroshio region -240 -120 10 65
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_08.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md This title already exists in the notebook markdown – appending an additional figure. Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. Gulf Stream region -95 -5 10 65
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_09.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md This title already exists in the notebook markdown – appending an additional figure. Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully. Argentine Basin -60 -30 -60 -30
Saved /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability_10.png Adding entry to per-notebook markdown: /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md This title already exists in the notebook markdown – appending an additional figure. Lines appended to /g/data/tm70/cyb561/access-om3-paper-1-runs/MC_25km_jra_ryf+wombatlite-test3-f4d79e82/notebooks/mkfigs_output_MC_25km_jra_ryf+wombatlite-test3-f4d79e82/mkmd/Currents_streamfunction_variability.md successfully.
