An Exception was encountered at 'In [16]'.
ACCESS-OM3 and ACCESS-OM2 timeseries¶
Includes both scalar variables and spatial averages of 2d variables.
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-27-07-2026-PD-control/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-27-07-2026-PD-control/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 = "/jobfs/181533135.gadi-pbs/MC_8km_jra_ryf-test2-expt-c349f800/datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-live/access-om3-paper-1/notebooks/mkfigs_output_MC_8km_jra_ryf-test2-expt-c349f800/"
nbname = "timeseries.ipynb"
# Parameters
esm_file = "/jobfs/181533135.gadi-pbs/MC_8km_jra_ryf-test2-expt-c349f800/datastore.json"
papermill = True
cwd = "/g/data/tm70/cyb561/access-om3-live/access-om3-paper-1/notebooks/mkfigs_output_MC_8km_jra_ryf-test2-expt-c349f800/"
nbname = "timeseries.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, get_exptname_from_path
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, get_exptname_from_path
expt_key, info = guess_experiment_from_esm_file(esm_file)
model_name = info["model"] # OM3 or CM3
In [4]:
Copied!
%matplotlib inline
from model_agnostic import patch_broken_conda_env, patch_dask_workers
patch_broken_conda_env() # conda env workaround; no-op on healthy envs
import pandas as pd
import intake
import xarray as xr
import numpy as np
import cf_xarray as cfxr
import pint_xarray
from pint import application_registry as ureg
import cf_xarray.units
import cftime
import xgcm
import os
import matplotlib.pyplot as plt
import cmocean as cm
import cartopy.crs as ccrs
import cartopy.feature as cft
from dask.distributed import Client
%matplotlib inline
from model_agnostic import patch_broken_conda_env, patch_dask_workers
patch_broken_conda_env() # conda env workaround; no-op on healthy envs
import pandas as pd
import intake
import xarray as xr
import numpy as np
import cf_xarray as cfxr
import pint_xarray
from pint import application_registry as ureg
import cf_xarray.units
import cftime
import xgcm
import os
import matplotlib.pyplot as plt
import cmocean as cm
import cartopy.crs as ccrs
import cartopy.feature as cft
from dask.distributed import Client
In [5]:
Copied!
xr.set_options(keep_attrs=True); # cf_xarray works best when xarray keeps attributes by default
xr.set_options(keep_attrs=True); # cf_xarray works best when xarray keeps attributes by default
In [6]:
Copied!
client = Client(threads_per_worker=1)
patch_dask_workers(client) # patch workers too
client
client = Client(threads_per_worker=1)
patch_dask_workers(client) # patch workers too
client
Out[6]:
Client
Client-88dae4bb-cb46-11f1-b7c1-00000081fe80
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: http://127.0.0.1:8787/status |
Cluster Info
LocalCluster
a4fa9276
| 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-b9344652-c504-41d5-9c78-00e9207d5c29
| Comm: tcp://127.0.0.1:37157 | 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:38665 | Total threads: 1 |
| Dashboard: http://127.0.0.1:46093/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:42503 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-7s0e1fwp | |
Worker: 1
| Comm: tcp://127.0.0.1:37807 | Total threads: 1 |
| Dashboard: http://127.0.0.1:41547/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:39797 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-m9hdhosq | |
Worker: 2
| Comm: tcp://127.0.0.1:46715 | Total threads: 1 |
| Dashboard: http://127.0.0.1:34745/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:38433 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-i_gkbas1 | |
Worker: 3
| Comm: tcp://127.0.0.1:36977 | Total threads: 1 |
| Dashboard: http://127.0.0.1:45869/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:44901 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-qbwpg5ch | |
Worker: 4
| Comm: tcp://127.0.0.1:46733 | Total threads: 1 |
| Dashboard: http://127.0.0.1:36577/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:34015 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-sgugx6dx | |
Worker: 5
| Comm: tcp://127.0.0.1:40345 | Total threads: 1 |
| Dashboard: http://127.0.0.1:46337/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:42123 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-5gdijwgi | |
Worker: 6
| Comm: tcp://127.0.0.1:37957 | Total threads: 1 |
| Dashboard: http://127.0.0.1:43389/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:35723 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-lq7jfbmq | |
Worker: 7
| Comm: tcp://127.0.0.1:46119 | Total threads: 1 |
| Dashboard: http://127.0.0.1:45181/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:36745 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-xxb65u2n | |
Worker: 8
| Comm: tcp://127.0.0.1:43917 | Total threads: 1 |
| Dashboard: http://127.0.0.1:39809/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:33771 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-v46hfnib | |
Worker: 9
| Comm: tcp://127.0.0.1:43183 | Total threads: 1 |
| Dashboard: http://127.0.0.1:33065/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:40749 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-ndgyk69x | |
Worker: 10
| Comm: tcp://127.0.0.1:37295 | Total threads: 1 |
| Dashboard: http://127.0.0.1:37609/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:42705 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-rw0qnusk | |
Worker: 11
| Comm: tcp://127.0.0.1:41959 | Total threads: 1 |
| Dashboard: http://127.0.0.1:44197/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:36925 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-62wkdjfp | |
Worker: 12
| Comm: tcp://127.0.0.1:33167 | Total threads: 1 |
| Dashboard: http://127.0.0.1:43575/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:40107 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-lgizwer3 | |
Worker: 13
| Comm: tcp://127.0.0.1:42781 | Total threads: 1 |
| Dashboard: http://127.0.0.1:45689/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:45517 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-n113t0c9 | |
Worker: 14
| Comm: tcp://127.0.0.1:40685 | Total threads: 1 |
| Dashboard: http://127.0.0.1:33185/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:42113 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-k5n52l9m | |
Worker: 15
| Comm: tcp://127.0.0.1:40737 | Total threads: 1 |
| Dashboard: http://127.0.0.1:45403/status | Memory: 11.78 GiB |
| Nanny: tcp://127.0.0.1:35597 | |
| Local directory: /jobfs/181533135.gadi-pbs/dask-scratch-space/worker-icrdyig8 | |
In [7]:
Copied!
thisdir = './'
thisdir = './'
In [8]:
Copied!
IAF = esm_file.find('iaf') > 0
IAF
IAF = esm_file.find('iaf') > 0
IAF
Out[8]:
False
Define data sources¶
In [9]:
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',
'/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-00532b88/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-COARE-12fab415/datastore.json',
esm_file, # the CM3/OM3 experiment under test -- already includes CM3 when esm_file is a CM3 path
"/g/data/ol01/outputs/access-om3-100km/MC_100km_jra_ryf+wombatlite-1e74abf-11f9df5c/experiment_datastore.json", # 100km RYF spin-up
]
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',
'/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-00532b88/datastore.json',
# '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-COARE-12fab415/datastore.json',
esm_file, # the CM3/OM3 experiment under test -- already includes CM3 when esm_file is a CM3 path
"/g/data/ol01/outputs/access-om3-100km/MC_100km_jra_ryf+wombatlite-1e74abf-11f9df5c/experiment_datastore.json", # 100km RYF spin-up
]
catalogs
Out[9]:
['/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf-1.0-beta-5165c0f8/datastore.json', '/g/data/ol01/outputs/access-om3-25km/MC_25km_jra_iaf+wombatlite-test3v2-00532b88/datastore.json', '/jobfs/181533135.gadi-pbs/MC_8km_jra_ryf-test2-expt-c349f800/datastore.json', '/g/data/ol01/outputs/access-om3-100km/MC_100km_jra_ryf+wombatlite-1e74abf-11f9df5c/experiment_datastore.json']
In [10]:
Copied!
datastores = { get_exptname_from_path(c):
intake.open_esm_datastore(c, columns_with_iterables=['variable'])
for c in catalogs }
datastores
datastores = { get_exptname_from_path(c):
intake.open_esm_datastore(c, columns_with_iterables=['variable'])
for c in catalogs }
datastores
Out[10]:
{'MC_25km_jra_iaf-1.0-beta-5165c0f8': <datastore catalog with 22 dataset(s) from 6042 asset(s)>,
'MC_25km_jra_iaf+wombatlite-test3v2-00532b88': <datastore catalog with 24 dataset(s) from 32927 asset(s)>,
'MC_8km_jra_ryf-test2-expt-c349f800': <datastore catalog with 26 dataset(s) from 1948 asset(s)>,
'MC_100km_jra_ryf+wombatlite-1e74abf-11f9df5c': <experiment_datastore catalog with 15 dataset(s) from 8540 asset(s)>}
In [11]:
Copied!
if model_name == "ACCESS-CM3":
# ACCESS-CM2 cj877 control run reference (as in other CM3 notebooks).
# Its scalar diagnostics use the same variable names as the OM2 catalog.
om2expt = 'cj877'
import glob
_cm2_files = sorted(glob.glob('/g/data/p73/archive/CMIP6/non-CMIP/ACCESS-CM2/cj877/history/ocn/ocean-scalar-1-daily-ym_*.nc'))
om2cat = xr.open_mfdataset(_cm2_files, combine="by_coords", parallel=True, chunks={"time": 181},
decode_timedelta=False, use_cftime=True,
data_vars="minimal", coords="minimal", compat="override")
else:
# ACCESS-OM2-025 experiment
# use control - see https://forum.access-hive.org.au/t/access-om2-control-experiments/258
if IAF:
om2expt = '025deg_jra55_iaf_omip2_cycle1'
else:
om2expt = '025deg_jra55_ryf9091_gadi'
om2cat = intake.cat.access_nri[om2expt]
if model_name == "ACCESS-CM3":
# ACCESS-CM2 cj877 control run reference (as in other CM3 notebooks).
# Its scalar diagnostics use the same variable names as the OM2 catalog.
om2expt = 'cj877'
import glob
_cm2_files = sorted(glob.glob('/g/data/p73/archive/CMIP6/non-CMIP/ACCESS-CM2/cj877/history/ocn/ocean-scalar-1-daily-ym_*.nc'))
om2cat = xr.open_mfdataset(_cm2_files, combine="by_coords", parallel=True, chunks={"time": 181},
decode_timedelta=False, use_cftime=True,
data_vars="minimal", coords="minimal", compat="override")
else:
# ACCESS-OM2-025 experiment
# use control - see https://forum.access-hive.org.au/t/access-om2-control-experiments/258
if IAF:
om2expt = '025deg_jra55_iaf_omip2_cycle1'
else:
om2expt = '025deg_jra55_ryf9091_gadi'
om2cat = intake.cat.access_nri[om2expt]
Load grid data¶
This is needed for xgcm so we can get spatial averages of non-scalar data.
In [12]:
Copied!
exptname = get_exptname_from_path(esm_file)
# filename regex matches om3 and cm3 static files
static = xr.open_dataset(datastores[exptname].search(filename=".*mom6.static.nc").df.loc[0, "path"]).drop_vars("time")
exptname = get_exptname_from_path(esm_file)
# filename regex matches om3 and cm3 static files
static = xr.open_dataset(datastores[exptname].search(filename=".*mom6.static.nc").df.loc[0, "path"]).drop_vars("time")
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/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 'requests.exceptions.ConnectionError'>
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bcbaf753d0>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/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 'requests.exceptions.ConnectionError'>
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bbf8d23020>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/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 'requests.exceptions.ConnectionError'>
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bbf9321070>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/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 'requests.exceptions.ConnectionError'>
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bbf919d340>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/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 'requests.exceptions.ConnectionError'>
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bbf8d233e0>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/lib/python3.12/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'argo' loading failed:
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bbf917bef0>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/lib/python3.12/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'erddapy' loading failed:
HTTPSConnectionPool(host='raw.githubusercontent.com', port=443): Max retries exceeded with url: /IrishMarineInstitute/awesome-erddap/master/erddaps.json (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x14bbf3a348c0>: Failed to establish a new connection: [Errno 101] Network is unreachable'))
external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
WARNING:pint.util:Redefining 'degrees_north' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degrees_N' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degreesN' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degree_north' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degree_N' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degreeN' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degrees_east' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degrees_E' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degreesE' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degree_east' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degree_E' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'degreeE' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
WARNING:pint.util:Redefining 'gpm' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
In [13]:
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': 'zl', 'outer': 'zi'}
}
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': 'zl', 'outer': 'zi'}
}
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': 'zl', 'outer': 'zi'}
}
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': 'zl', 'outer': 'zi'}
}
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
)
Define list of variables and load them¶
In [14]:
Copied!
fields = { # key is OM3 name, value is OM2 name
# 'mlotst': 'mld',
# keys from ncdump -h ctrl_run_25km_0.5/archive/output009/access-om3.mom6.scalar.1day.snap.1962.nc | grep double
'masso': None, # Mass of liquid ocean
'soga': 'salt_global_ave',
'thetaoga': 'temp_global_ave', # BUG: is OM2 potential temperature? https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/3#issuecomment-3736872326
'tosga': 'temp_surface_ave', # BUG: is OM2 potential temperature? https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/3#issuecomment-3736872326
'sosga': 'salt_surface_ave',
'speed_max': None,
'mlotst_max': None,
'tos_max': None,
'tos_min': None,
'sos_max': None,
'sos_min': None,
'zos_max': None,
'zos_min': None,
'total_salt_Flux_Added': None,
'total_salt_Flux_In': None,
'total_salt_flux': None,
'net_fresh_water_global_adjustment': None,
'salt_flux_global_restoring_adjustment': None,
'total_wfo': None,
'total_evs': None,
'total_fsitherm': None,
'total_precip': None,
'total_prsn': None,
'total_lprec': None,
'total_ficeberg': None,
'total_friver': None,
'total_net_massout': None,
'total_net_massin': None,
}
fields
fields = { # key is OM3 name, value is OM2 name
# 'mlotst': 'mld',
# keys from ncdump -h ctrl_run_25km_0.5/archive/output009/access-om3.mom6.scalar.1day.snap.1962.nc | grep double
'masso': None, # Mass of liquid ocean
'soga': 'salt_global_ave',
'thetaoga': 'temp_global_ave', # BUG: is OM2 potential temperature? https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/3#issuecomment-3736872326
'tosga': 'temp_surface_ave', # BUG: is OM2 potential temperature? https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/3#issuecomment-3736872326
'sosga': 'salt_surface_ave',
'speed_max': None,
'mlotst_max': None,
'tos_max': None,
'tos_min': None,
'sos_max': None,
'sos_min': None,
'zos_max': None,
'zos_min': None,
'total_salt_Flux_Added': None,
'total_salt_Flux_In': None,
'total_salt_flux': None,
'net_fresh_water_global_adjustment': None,
'salt_flux_global_restoring_adjustment': None,
'total_wfo': None,
'total_evs': None,
'total_fsitherm': None,
'total_precip': None,
'total_prsn': None,
'total_lprec': None,
'total_ficeberg': None,
'total_friver': None,
'total_net_massout': None,
'total_net_massin': None,
}
fields
Out[14]:
{'masso': None,
'soga': 'salt_global_ave',
'thetaoga': 'temp_global_ave',
'tosga': 'temp_surface_ave',
'sosga': 'salt_surface_ave',
'speed_max': None,
'mlotst_max': None,
'tos_max': None,
'tos_min': None,
'sos_max': None,
'sos_min': None,
'zos_max': None,
'zos_min': None,
'total_salt_Flux_Added': None,
'total_salt_Flux_In': None,
'total_salt_flux': None,
'net_fresh_water_global_adjustment': None,
'salt_flux_global_restoring_adjustment': None,
'total_wfo': None,
'total_evs': None,
'total_fsitherm': None,
'total_precip': None,
'total_prsn': None,
'total_lprec': None,
'total_ficeberg': None,
'total_friver': None,
'total_net_massout': None,
'total_net_massin': None}
In [15]:
Copied!
# from https://github.com/COSIMA/cosima-recipes/blob/main/Cooking-Lessons-101-Tutorials/ACCESS-NRI_Intake_Catalog.ipynb
# Speeding up opening your datasets
# Try passing the following argument to your to_dask or to_dataset_dict call:
# See the xarray documentation on Reading multi-file datasets for more details about these arguments.
# https://docs.xarray.dev/en/stable/user-guide/io.html#reading-multi-file-datasets
xarray_combine_by_coords_kwargs=dict(
compat="override",
data_vars="minimal",
coords="minimal"
)
# from https://github.com/COSIMA/cosima-recipes/blob/main/Cooking-Lessons-101-Tutorials/ACCESS-NRI_Intake_Catalog.ipynb
# Speeding up opening your datasets
# Try passing the following argument to your to_dask or to_dataset_dict call:
# See the xarray documentation on Reading multi-file datasets for more details about these arguments.
# https://docs.xarray.dev/en/stable/user-guide/io.html#reading-multi-file-datasets
xarray_combine_by_coords_kwargs=dict(
compat="override",
data_vars="minimal",
coords="minimal"
)
Execution using papermill encountered an exception here and stopped:
In [16]:
Copied!
# load data, with lots of workarounds for disambiguation
data = dict()
for k, om2 in fields.items():
print(k)
# d = [ ds.search(variable=k).to_dataset_dict(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for ds in datastores ]
try:
d = { n: ds.search(variable=k) for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except ValueError:
try:
d = { n: ds.search(variable=k, frequency='1mon', variable_cell_methods='.*time: mean.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: mean.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: point.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: min.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: max.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
print(f'{k} failed')
continue
if om2 is None: # add reference model data, with variable name changed to match om3
d2 = {}
elif model_name == "ACCESS-CM3":
cm2_da = om2cat[[om2]].isel(scalar_axis=0, drop=True).rename({om2: k}) # drop dummy scalar_axis dim
# relabel cj877 (cftime, model years) onto CM3's own datetime64 time axis,
# matching elapsed time since each run's start (both are daily)
_cm3t = d[exptname].time
n = min(cm2_da.sizes["time"], _cm3t.sizes["time"])
cm2_da = cm2_da.isel(time=slice(None, n)).assign_coords(time=("time", _cm3t.isel(time=slice(None, n)).values))
d2 = {om2expt: cm2_da}
else:
# d[om2expt] = om2cat.search(variable=om2).to_dask().rename({om2: k})
d2 = {om2expt: om2cat.search(variable=om2, frequency=d[list(d.keys())[0]].attrs['intake_esm_attrs:frequency']).to_dask().rename({om2: k})}
# if d:
# data[k] = xr.concat(d, 'time')
# else:
# print(f'no data for {k}')
# break
data[k] = d2 | d
# break
if IAF:
for vname, d in data.items():
for expt in d.keys():
try:
d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True)
except KeyError:
pass
# RYF spin-up runs use a synthetic calendar (e.g. starting ~1900s) unrelated to real
# forcing dates, so their time axis does not overlap the IAF experiments at all.
# Relabel onto the primary experiment's real time axis so the RYF run plots on the
# same timescale as the IAF runs. Only applied to the RYF experiment.
_ryf_expt = "MC_100km_jra_ryf+wombatlite-1e74abf-11f9df5c"
for _vname, _d in data.items():
if _ryf_expt not in _d or exptname not in _d:
continue
_ref = _d[exptname].time
_rds = _d[_ryf_expt]
_n = min(_rds.sizes["time"], _ref.sizes["time"])
_d[_ryf_expt] = (_rds.isel(time=slice(None, _n))
.assign_coords(time=("time", _ref.isel(time=slice(None, _n)).values)))
data
# load data, with lots of workarounds for disambiguation
data = dict()
for k, om2 in fields.items():
print(k)
# d = [ ds.search(variable=k).to_dataset_dict(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for ds in datastores ]
try:
d = { n: ds.search(variable=k) for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except ValueError:
try:
d = { n: ds.search(variable=k, frequency='1mon', variable_cell_methods='.*time: mean.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: mean.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: point.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: min.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
try:
d = { n: ds.search(variable=k, variable_cell_methods='.*time: max.*') for n, ds in datastores.items() }
d = { n: ds.to_dask(xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs) for n, ds in d.items() if ds } # avoid .to_dask() for empty datasets
except:
print(f'{k} failed')
continue
if om2 is None: # add reference model data, with variable name changed to match om3
d2 = {}
elif model_name == "ACCESS-CM3":
cm2_da = om2cat[[om2]].isel(scalar_axis=0, drop=True).rename({om2: k}) # drop dummy scalar_axis dim
# relabel cj877 (cftime, model years) onto CM3's own datetime64 time axis,
# matching elapsed time since each run's start (both are daily)
_cm3t = d[exptname].time
n = min(cm2_da.sizes["time"], _cm3t.sizes["time"])
cm2_da = cm2_da.isel(time=slice(None, n)).assign_coords(time=("time", _cm3t.isel(time=slice(None, n)).values))
d2 = {om2expt: cm2_da}
else:
# d[om2expt] = om2cat.search(variable=om2).to_dask().rename({om2: k})
d2 = {om2expt: om2cat.search(variable=om2, frequency=d[list(d.keys())[0]].attrs['intake_esm_attrs:frequency']).to_dask().rename({om2: k})}
# if d:
# data[k] = xr.concat(d, 'time')
# else:
# print(f'no data for {k}')
# break
data[k] = d2 | d
# break
if IAF:
for vname, d in data.items():
for expt in d.keys():
try:
d[expt] = d[expt].convert_calendar("proleptic_gregorian", use_cftime=True)
except KeyError:
pass
# RYF spin-up runs use a synthetic calendar (e.g. starting ~1900s) unrelated to real
# forcing dates, so their time axis does not overlap the IAF experiments at all.
# Relabel onto the primary experiment's real time axis so the RYF run plots on the
# same timescale as the IAF runs. Only applied to the RYF experiment.
_ryf_expt = "MC_100km_jra_ryf+wombatlite-1e74abf-11f9df5c"
for _vname, _d in data.items():
if _ryf_expt not in _d or exptname not in _d:
continue
_ref = _d[exptname].time
_rds = _d[_ryf_expt]
_n = min(_rds.sizes["time"], _ref.sizes["time"])
_d[_ryf_expt] = (_rds.isel(time=slice(None, _n))
.assign_coords(time=("time", _ref.isel(time=slice(None, _n)).values)))
data
masso
soga
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[16], line 44 40 cm2_da = cm2_da.isel(time=slice(None, n)).assign_coords(time=("time", _cm3t.isel(time=slice(None, n)).values)) 41 d2 = {om2expt: cm2_da} 42 else: 43 # d[om2expt] = om2cat.search(variable=om2).to_dask().rename({om2: k}) ---> 44 d2 = {om2expt: om2cat.search(variable=om2, frequency=d[list(d.keys())[0]].attrs['intake_esm_attrs:frequency']).to_dask().rename({om2: k})} 45 # if d: 46 # data[k] = xr.concat(d, 'time') 47 # else: File /g/data/xp65/public/apps/med_conda/envs/analysis3-26.08/lib/python3.12/site-packages/intake_esm/core.py:926, in esm_datastore.to_dask(self, **kwargs) 924 if len(self) != 1: # quick check to fail more quickly if there are many results 925 lens = self._aggr_lengths --> 926 raise ValueError( 927 f'Expected exactly one dataset. Received {len(self)} datasets. Please refine your search on {", ".join(lens.keys())} or use `.to_dataset_dict()`.' 928 ) 929 res = self.to_dataset_dict(**{**kwargs, 'progressbar': False}) 930 if len(res) != 1: # extra check in case kwargs did modify something ValueError: Expected exactly one dataset. Received 0 datasets. Please refine your search on or use `.to_dataset_dict()`.
In [ ]:
Copied!
data.keys()
data.keys()
In [ ]:
Copied!
data
data
Plot timeseries¶
In [ ]:
Copied!
# set date range
# trange = slice(cftime.DatetimeNoLeap(1900, 1, 1, 0, 0, 0, 0),
# cftime.DatetimeNoLeap(1960, 1, 1, 0, 0, 0, 0))
trange = slice(None, None) # all available
# set date range
# trange = slice(cftime.DatetimeNoLeap(1900, 1, 1, 0, 0, 0, 0),
# cftime.DatetimeNoLeap(1960, 1, 1, 0, 0, 0, 0))
trange = slice(None, None) # all available
In [ ]:
Copied!
for k, d in data.items():
print(k)
fname = thisdir+k+'_'+list(d.values())[-1][k].attrs['long_name'].replace(' ', '_').replace('/','_')+'.png'
# print(fname)
if False: #os.path.isfile(fname):
print(f'---- skipping existing file {fname}')
else:
plt.figure(figsize=(10,5))
ax = plt.subplot(111)
for n, ds in d.items():
dat = ds[k].sel(time=trange)
if dat.cell_methods.find('time: mean') >= 0:
# calculate spatial average if needed
if 'depth' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y', 'Z'])
elif 'longitude' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y'])
elif dat.cell_methods.find('time: max') >= 0:
dat = dat.cf.max(dim=['longitude', 'latitude'])
elif dat.cell_methods.find('time: min') >= 0:
dat = dat.cf.min(dim=['longitude', 'latitude'])
# dat = dat - dat.isel(time=0)
dat.load()
# if int((dat.time[1]-dat.time[0]).values/1e9/60/60/24) == 1:
# label = 'daily'
# else:
# label = 'monthly mean' # possible BUG: plausible guess
# p = dat.plot(label=' '.join([n, label]), linewidth=1)
p = dat.plot(linewidth=.7)
dat.cf.resample(time='1YE').mean('time').plot(
color=p[-1].get_color(),
linewidth=2,
label=n
)
# dat.cf.rolling(time=12, center=True).mean('time').plot()
# dat.cf.rolling(time='1YE', center=True).mean('time').plot()
# NB: title from first item (OM3), but y axis label from last item (may be OM2)
# plt.title(k+': '+list(d.values())[0][k].attrs['long_name'])
plt.title(f'{k}: {dat.attrs["long_name"]}')
plt.xlabel('')
plt.autoscale(enable=True, axis='x', tight=True)
plt.legend()
# # Shrink current axis's height from the bottom
# box = ax.get_position()
# ax.set_position([box.x0, box.y0 + box.height * 0.4,
# box.width, box.height * 0.6])
# # Put a legend below current axis
# ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.08))
# # try:
# plt.savefig(fname, dpi=150)
# print(f'saved {fname}')
# except FileNotFoundError:
# print(f'*** FileNotFoundError when saving {fname}')
# break
mkmd.savefig(plt.gcf(), "Timeseries", "ACCESS-OM3 global model timeseries. [GitHub issue: Global average annual mean time series](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/3)")
for k, d in data.items():
print(k)
fname = thisdir+k+'_'+list(d.values())[-1][k].attrs['long_name'].replace(' ', '_').replace('/','_')+'.png'
# print(fname)
if False: #os.path.isfile(fname):
print(f'---- skipping existing file {fname}')
else:
plt.figure(figsize=(10,5))
ax = plt.subplot(111)
for n, ds in d.items():
dat = ds[k].sel(time=trange)
if dat.cell_methods.find('time: mean') >= 0:
# calculate spatial average if needed
if 'depth' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y', 'Z'])
elif 'longitude' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y'])
elif dat.cell_methods.find('time: max') >= 0:
dat = dat.cf.max(dim=['longitude', 'latitude'])
elif dat.cell_methods.find('time: min') >= 0:
dat = dat.cf.min(dim=['longitude', 'latitude'])
# dat = dat - dat.isel(time=0)
dat.load()
# if int((dat.time[1]-dat.time[0]).values/1e9/60/60/24) == 1:
# label = 'daily'
# else:
# label = 'monthly mean' # possible BUG: plausible guess
# p = dat.plot(label=' '.join([n, label]), linewidth=1)
p = dat.plot(linewidth=.7)
dat.cf.resample(time='1YE').mean('time').plot(
color=p[-1].get_color(),
linewidth=2,
label=n
)
# dat.cf.rolling(time=12, center=True).mean('time').plot()
# dat.cf.rolling(time='1YE', center=True).mean('time').plot()
# NB: title from first item (OM3), but y axis label from last item (may be OM2)
# plt.title(k+': '+list(d.values())[0][k].attrs['long_name'])
plt.title(f'{k}: {dat.attrs["long_name"]}')
plt.xlabel('')
plt.autoscale(enable=True, axis='x', tight=True)
plt.legend()
# # Shrink current axis's height from the bottom
# box = ax.get_position()
# ax.set_position([box.x0, box.y0 + box.height * 0.4,
# box.width, box.height * 0.6])
# # Put a legend below current axis
# ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.08))
# # try:
# plt.savefig(fname, dpi=150)
# print(f'saved {fname}')
# except FileNotFoundError:
# print(f'*** FileNotFoundError when saving {fname}')
# break
mkmd.savefig(plt.gcf(), "Timeseries", "ACCESS-OM3 global model timeseries. [GitHub issue: Global average annual mean time series](https://github.com/ACCESS-Community-Hub/access-om3-paper-1/issues/3)")
In [ ]:
Copied!
# just plot soga, omitting 25km-iaf-test-for-AK-expt-7df5ef4c
for k, d in data.items():
if k != "soga":
continue
print(k)
fname = thisdir+k+'_'+list(d.values())[-1][k].attrs['long_name'].replace(' ', '_').replace('/','_')+'.png'
# print(fname)
if False: #os.path.isfile(fname):
print(f'---- skipping existing file {fname}')
else:
plt.figure(figsize=(10,5))
ax = plt.subplot(111)
for n, ds in d.items():
if n == "25km-iaf-test-for-AK-expt-7df5ef4c":
continue
dat = ds[k].sel(time=trange)
if dat.cell_methods.find('time: mean') >= 0:
# calculate spatial average if needed
if 'depth' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y', 'Z'])
elif 'longitude' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y'])
elif dat.cell_methods.find('time: max') >= 0:
dat = dat.cf.max(dim=['longitude', 'latitude'])
elif dat.cell_methods.find('time: min') >= 0:
dat = dat.cf.min(dim=['longitude', 'latitude'])
# dat = dat - dat.isel(time=0)
dat.load()
# if int((dat.time[1]-dat.time[0]).values/1e9/60/60/24) == 1:
# label = 'daily'
# else:
# label = 'monthly mean' # possible BUG: plausible guess
# p = dat.plot(label=' '.join([n, label]), linewidth=1)
p = dat.plot(linewidth=.7)
dat.cf.resample(time='1YE').mean('time').plot(
color=p[-1].get_color(),
linewidth=2,
label=n
)
# dat.cf.rolling(time=12, center=True).mean('time').plot()
# dat.cf.rolling(time='1YE', center=True).mean('time').plot()
# NB: title from first item (OM3), but y axis label from last item (may be OM2)
# plt.title(k+': '+list(d.values())[0][k].attrs['long_name'])
plt.title(f'{k}: {dat.attrs["long_name"]}')
plt.xlabel('')
plt.autoscale(enable=True, axis='x', tight=True)
plt.legend()
# # Shrink current axis's height from the bottom
# box = ax.get_position()
# ax.set_position([box.x0, box.y0 + box.height * 0.4,
# box.width, box.height * 0.6])
# # Put a legend below current axis
# ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.08))
# # try:
# plt.savefig(fname, dpi=150)
# print(f'saved {fname}')
# except FileNotFoundError:
# print(f'*** FileNotFoundError when saving {fname}')
# break
# just plot soga, omitting 25km-iaf-test-for-AK-expt-7df5ef4c
for k, d in data.items():
if k != "soga":
continue
print(k)
fname = thisdir+k+'_'+list(d.values())[-1][k].attrs['long_name'].replace(' ', '_').replace('/','_')+'.png'
# print(fname)
if False: #os.path.isfile(fname):
print(f'---- skipping existing file {fname}')
else:
plt.figure(figsize=(10,5))
ax = plt.subplot(111)
for n, ds in d.items():
if n == "25km-iaf-test-for-AK-expt-7df5ef4c":
continue
dat = ds[k].sel(time=trange)
if dat.cell_methods.find('time: mean') >= 0:
# calculate spatial average if needed
if 'depth' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y', 'Z'])
elif 'longitude' in dat.cf.coords: # https://xgcm.readthedocs.io/en/latest/grid_metrics.html?highlight=average#Grid-aware-(weighted)-average
dat = grid.average(dat, ['X', 'Y'])
elif dat.cell_methods.find('time: max') >= 0:
dat = dat.cf.max(dim=['longitude', 'latitude'])
elif dat.cell_methods.find('time: min') >= 0:
dat = dat.cf.min(dim=['longitude', 'latitude'])
# dat = dat - dat.isel(time=0)
dat.load()
# if int((dat.time[1]-dat.time[0]).values/1e9/60/60/24) == 1:
# label = 'daily'
# else:
# label = 'monthly mean' # possible BUG: plausible guess
# p = dat.plot(label=' '.join([n, label]), linewidth=1)
p = dat.plot(linewidth=.7)
dat.cf.resample(time='1YE').mean('time').plot(
color=p[-1].get_color(),
linewidth=2,
label=n
)
# dat.cf.rolling(time=12, center=True).mean('time').plot()
# dat.cf.rolling(time='1YE', center=True).mean('time').plot()
# NB: title from first item (OM3), but y axis label from last item (may be OM2)
# plt.title(k+': '+list(d.values())[0][k].attrs['long_name'])
plt.title(f'{k}: {dat.attrs["long_name"]}')
plt.xlabel('')
plt.autoscale(enable=True, axis='x', tight=True)
plt.legend()
# # Shrink current axis's height from the bottom
# box = ax.get_position()
# ax.set_position([box.x0, box.y0 + box.height * 0.4,
# box.width, box.height * 0.6])
# # Put a legend below current axis
# ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.08))
# # try:
# plt.savefig(fname, dpi=150)
# print(f'saved {fname}')
# except FileNotFoundError:
# print(f'*** FileNotFoundError when saving {fname}')
# break
Water balance checks¶
In [ ]:
Copied!
esm_name = get_exptname_from_path(esm_file)
esm_name
esm_name = get_exptname_from_path(esm_file)
esm_name
In [ ]:
Copied!
masked_areacello = static['areacello'] + 0*data['tos_max'][esm_name].tos_max.isel(time=0)
masked_areacello = static['areacello'] + 0*data['tos_max'][esm_name].tos_max.isel(time=0)
In [ ]:
Copied!
masked_areacello.plot()
masked_areacello.plot()
In [ ]:
Copied!
total_ocean_area = masked_areacello.sum().load()
total_ocean_area = masked_areacello.sum().load()
In [ ]:
Copied!
total_fsitherm = data['total_fsitherm'][esm_name]['total_fsitherm'].load()
total_fsitherm = data['total_fsitherm'][esm_name]['total_fsitherm'].load()
In [ ]:
Copied!
total_fsitherm.plot()
total_fsitherm.plot()
In [ ]:
Copied!
try:
net_fresh_water_global_adjustment = (data['net_fresh_water_global_adjustment'][esm_name]['net_fresh_water_global_adjustment'] * total_ocean_area).load()
net_fresh_water_global_adjustment.plot()
(net_fresh_water_global_adjustment - total_fsitherm).plot()
except KeyError:
pass
try:
net_fresh_water_global_adjustment = (data['net_fresh_water_global_adjustment'][esm_name]['net_fresh_water_global_adjustment'] * total_ocean_area).load()
net_fresh_water_global_adjustment.plot()
(net_fresh_water_global_adjustment - total_fsitherm).plot()
except KeyError:
pass