Horizontal Mean difference and RMS: model versus observations
mom6_tools.horizontalMean collection of functions for computing horizontal mean of difference and rms (model versus obs). This notebook servers as an example on how to compute the following operations:
\[diff(t,z)= A_{TOT}(z)^{-1}\sum_{i=1}^n (y_i(z) - \hat{y_i(x,y,z)}) A_i(z),\]
\[rms(t,z)= [A_{TOT}(z)^{-1}\sum_{i=1}^n (y_i(z) - \hat{y_i}(x,y,z))^2 A_i(z)]^{1/2},\]
where \(y\)(z) is the model output at point \(i\) and level \(z\), \(\hat{y}(z)\) is the observation at point \(i\) and level \(z\), \(n\) is the total number of grid points in the horizontal (i.e., NX x NY), \(A_{i}(z)\) is the area of grid cell \(i\) at level \(z\), and \(A_{TOT}(z) = \sum_{i=1}^n A_i(z)\) is the total ocean area at level z.
Important:
With the porpuses of calculating T and S changes at specific regions, \(A_{i}(z)\) is multiplied by basin masks generated via mom6_tools.m6toolbox.genBasinMasks. See notebook showing how to generate these masks.
[1]:
%load_ext autoreload
%autoreload 2
[2]:
import warnings
warnings.filterwarnings("ignore")
from mom6_tools.MOM6grid import MOM6grid
from mom6_tools.drift import HorizontalMeanDiff_da
from mom6_tools.m6plot import ztplot
from mom6_tools.jobqueue import get_cluster
from mom6_tools.m6toolbox import genBasinMasks
from mom6_tools.m6toolbox import weighted_temporal_mean, add_global_attrs
from mom6_tools.m6toolbox import cime_xmlquery
from IPython.display import display, Markdown, Latex
import yaml, intake, os
import xarray as xr
import matplotlib
import numpy as np
%matplotlib inline
Basemap module not found. Some regional plots may not function properly
[3]:
# Read in the yaml file
diag_config_yml_path = "diag_config.yml"
diag_config_yml = yaml.load(open(diag_config_yml_path,'r'), Loader=yaml.Loader)
[4]:
caseroot = diag_config_yml['Case']['CASEROOT']
casename = cime_xmlquery(caseroot, 'CASE')
DOUT_S = cime_xmlquery(caseroot, 'DOUT_S')
if DOUT_S:
OUTDIR = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ocn/hist/'
else:
OUTDIR = cime_xmlquery(caseroot, 'RUNDIR')
print('Output directory is:', OUTDIR)
print('Casename is:', casename)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[4], line 2
1 caseroot = diag_config_yml['Case']['CASEROOT']
----> 2 casename = cime_xmlquery(caseroot, 'CASE')
3 DOUT_S = cime_xmlquery(caseroot, 'DOUT_S')
4 if DOUT_S:
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/mom6_tools/m6toolbox.py:47, in cime_xmlquery(caseroot, varname)
45 """run CIME's xmlquery for varname in the directory caseroot, return the value"""
46 try:
---> 47 value = subprocess.check_output(
48 ["./xmlquery", "-N", "--value", varname],
49 stderr=subprocess.STDOUT,
50 cwd=caseroot,
51 )
52 except subprocess.CalledProcessError:
53 value = subprocess.check_output(
54 ["./xmlquery", "--value", varname], stderr=subprocess.STDOUT, cwd=caseroot
55 )
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:421, in check_output(timeout, *popenargs, **kwargs)
418 empty = b''
419 kwargs['input'] = empty
--> 421 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
422 **kwargs).stdout
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:503, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
500 kwargs['stdout'] = PIPE
501 kwargs['stderr'] = PIPE
--> 503 with Popen(*popenargs, **kwargs) as process:
504 try:
505 stdout, stderr = process.communicate(input, timeout=timeout)
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:971, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize)
967 if self.text_mode:
968 self.stderr = io.TextIOWrapper(self.stderr,
969 encoding=encoding, errors=errors)
--> 971 self._execute_child(args, executable, preexec_fn, close_fds,
972 pass_fds, cwd, env,
973 startupinfo, creationflags, shell,
974 p2cread, p2cwrite,
975 c2pread, c2pwrite,
976 errread, errwrite,
977 restore_signals,
978 gid, gids, uid, umask,
979 start_new_session)
980 except:
981 # Cleanup if the child failed starting.
982 for f in filter(None, (self.stdin, self.stdout, self.stderr)):
File ~/.asdf/installs/python/3.10.20/lib/python3.10/subprocess.py:1863, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session)
1861 if errno_num != 0:
1862 err_msg = os.strerror(errno_num)
-> 1863 raise child_exception_type(errno_num, err_msg, err_filename)
1864 raise child_exception_type(err_msg)
FileNotFoundError: [Errno 2] No such file or directory: '/glade/work/gmarques/cesm.cases/G/g.e30_a07c_cesm.GJRAv4.TL319_t232_wgx3_hycom1_N75.2025.130/'
[5]:
# The following parameters must be set accordingly
######################################################
# create an empty class object
class args:
pass
# set avg dates
avg = diag_config_yml['Avg']
args.start_date = avg['start_date']
args.end_date = avg['end_date']
args.casename = casename
args.obs = "woa-2018-tx2_3v2-annual-all"
args.z = casename+diag_config_yml['Fnames']['z']
args.static = casename+diag_config_yml['Fnames']['static']
args.geom = casename+diag_config_yml['Fnames']['geom']
args.savefigs = False
args.nw = 6 # requesting 6 workers
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 13
11 args.start_date = avg['start_date']
12 args.end_date = avg['end_date']
---> 13 args.casename = casename
14 args.obs = "woa-2018-tx2_3v2-annual-all"
15 args.z = casename+diag_config_yml['Fnames']['z']
NameError: name 'casename' is not defined
[6]:
# read grid info
geom_file = OUTDIR+'/'+args.geom
if os.path.exists(geom_file):
grd = MOM6grid(OUTDIR+'/'+args.static, geom_file, xrformat=True)
else:
grd = MOM6grid(OUTDIR+'/'+args.static, xrformat=True)
try:
depth = grd.depth_ocean.values
except:
depth = grd.deptho.values
try:
area = grd.area_t.where(grd.wet > 0)
except:
area = grd.areacello.where(grd.wet > 0)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 2
1 # read grid info
----> 2 geom_file = OUTDIR+'/'+args.geom
3 if os.path.exists(geom_file):
4 grd = MOM6grid(OUTDIR+'/'+args.static, geom_file, xrformat=True)
NameError: name 'OUTDIR' is not defined
[7]:
# remote Nan's, otherwise genBasinMasks won't work
depth[np.isnan(depth)] = 0.0
basin_code = genBasinMasks(grd.geolon.values, grd.geolat.values, depth, xda=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 # remote Nan's, otherwise genBasinMasks won't work
----> 2 depth[np.isnan(depth)] = 0.0
3 basin_code = genBasinMasks(grd.geolon.values, grd.geolat.values, depth, xda=True)
NameError: name 'depth' is not defined
[8]:
parallel, cluster, client = get_cluster(args.nw, cluster_class='PBSCluster')
client
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[8], line 1
----> 1 parallel, cluster, client = get_cluster(args.nw, cluster_class='PBSCluster')
2 client
AttributeError: type object 'args' has no attribute 'nw'
[9]:
def preprocess(ds):
if 'thetao' not in ds.variables:
ds["thetao"] = xr.zeros_like(ds.h)
if 'so' not in ds.variables:
ds["so"] = xr.zeros_like(ds.h)
return ds
[10]:
# read dataset
ds = xr.open_mfdataset(OUTDIR+'/'+args.z,
parallel=True,
combine="nested", # concatenate in order of files
concat_dim="time", # concatenate along time
preprocess=preprocess,
).chunk({"time": 12})
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 2
1 # read dataset
----> 2 ds = xr.open_mfdataset(OUTDIR+'/'+args.z,
3 parallel=True,
4 combine="nested", # concatenate in order of files
5 concat_dim="time", # concatenate along time
6 preprocess=preprocess,
7 ).chunk({"time": 12})
NameError: name 'OUTDIR' is not defined
[11]:
# Compute thetao climatologies
var = 'thetao'
attrs = {
'description': 'Annual mean climatology for '+var,
'start_date': args.start_date,
'end_date': args.end_date,
'reduction_method': 'annual mean weighted by days in each month',
'casename': casename
}
thetao_model = weighted_temporal_mean(ds,var)
thetao_model.attrs = attrs
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 8
1 # Compute thetao climatologies
2 var = 'thetao'
3 attrs = {
4 'description': 'Annual mean climatology for '+var,
5 'start_date': args.start_date,
6 'end_date': args.end_date,
7 'reduction_method': 'annual mean weighted by days in each month',
----> 8 'casename': casename
9 }
11 thetao_model = weighted_temporal_mean(ds,var)
12 thetao_model.attrs = attrs
NameError: name 'casename' is not defined
[12]:
# Compute thetao climatologies
var = 'so'
attrs = {
'description': 'Annual mean climatology for '+var,
'start_date': args.start_date,
'end_date': args.end_date,
'reduction_method': 'annual mean weighted by days in each month',
'casename': casename
}
salt_model = weighted_temporal_mean(ds,var)
salt_model.attrs = attrs
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 8
1 # Compute thetao climatologies
2 var = 'so'
3 attrs = {
4 'description': 'Annual mean climatology for '+var,
5 'start_date': args.start_date,
6 'end_date': args.end_date,
7 'reduction_method': 'annual mean weighted by days in each month',
----> 8 'casename': casename
9 }
11 salt_model = weighted_temporal_mean(ds,var)
12 salt_model.attrs = attrs
NameError: name 'casename' is not defined
[13]:
# load obs
catalog = intake.open_catalog(diag_config_yml['oce_cat'])
obs = catalog[args.obs].to_dask()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[13], line 2
1 # load obs
----> 2 catalog = intake.open_catalog(diag_config_yml['oce_cat'])
3 obs = catalog[args.obs].to_dask()
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/__init__.py:186, in open_catalog(uri, **kwargs)
179 raise ValueError(
180 f"Unknown catalog driver '{driver}'. "
181 "Do you need to install a new driver from the plugin directory? "
182 "https://intake.readthedocs.io/en/latest/plugin-directory.html\n"
183 f"Current registry: {list(sorted(registry))}"
184 )
185 try:
--> 186 return registry[driver](uri, **kwargs)
187 except VersionError:
188 # warn that we are switching to V2? The file will be read twice
189 return from_yaml_file(uri, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:617, in YAMLFileCatalog.__init__(self, path, text, autoreload, **kwargs)
615 self.filesystem = kwargs.pop("fs", None)
616 self.access = "name" not in kwargs
--> 617 super(YAMLFileCatalog, self).__init__(**kwargs)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/base.py:128, in Catalog.__init__(self, entries, name, description, metadata, ttl, getenv, getshell, persist_mode, storage_options, user_parameters)
126 self.updated = time.time()
127 self._entries = entries if entries is not None else self._make_entries_container()
--> 128 self.force_reload()
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/base.py:186, in Catalog.force_reload(self)
184 """Imperative reload data now"""
185 self.updated = time.time()
--> 186 self._load()
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:652, in YAMLFileCatalog._load(self, reload)
650 logger.warning("Use of '!template' deprecated - fixing")
651 text = text.replace("!template ", "")
--> 652 self.parse(text)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:730, in YAMLFileCatalog.parse(self, text)
728 # Second, we validate the schema and semantics
729 context = dict(root=self._dir)
--> 730 result = CatalogParser(data, context=context, getenv=self.getenv, getshell=self.getshell)
731 if result.errors:
732 raise exceptions.ValidationError(
733 "Catalog '{}' has validation errors:\n\n{}"
734 "".format(self.path, "\n".join(result.errors)),
735 result.errors,
736 )
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:342, in CatalogParser.__init__(self, data, getenv, getshell, context)
340 self.getenv = getenv
341 self.getshell = getshell
--> 342 self._data = self._parse(data)
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:565, in CatalogParser._parse(self, data)
561 if (data.get("version", None) or data.get("metadata", {}).get("version", None) or 1) > 1:
562 raise VersionError("Not a V1 Catalog; perhaps use intake.open_catalog")
564 return dict(
--> 565 plugin_sources=self._parse_plugins(data),
566 data_sources=self._parse_data_sources(data),
567 metadata=data.get("metadata", {}),
568 name=data.get("name"),
569 description=data.get("description"),
570 )
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/catalog/local.py:400, in CatalogParser._parse_plugins(self, data)
397 elif "module" in plugin_source:
398 import intake
--> 400 intake.import_name(plugin_source["module"])
401 elif "dir" in plugin_source:
402 self.error(
403 "The key 'dir', and in general the feature of registering "
404 "plugins from a directory of Python scripts outside of "
405 "sys.path, is no longer supported. Use 'module'.",
406 plugin_source,
407 )
File ~/checkouts/readthedocs.org/user_builds/mom6-tools/envs/latest/lib/python3.10/site-packages/intake/utils.py:27, in import_name(name)
25 modname = name.split(":", 1)[0]
26 logger.debug("Importing: '%s'" % modname)
---> 27 mod = importlib.import_module(modname)
28 if ":" in name:
29 end = name.split(":")[1]
File ~/.asdf/installs/python/3.10.20/lib/python3.10/importlib/__init__.py:126, in import_module(name, package)
124 break
125 level += 1
--> 126 return _bootstrap._gcd_import(name[level:], package, level)
File <frozen importlib._bootstrap>:1050, in _gcd_import(name, package, level)
File <frozen importlib._bootstrap>:1027, in _find_and_load(name, import_)
File <frozen importlib._bootstrap>:1004, in _find_and_load_unlocked(name, import_)
ModuleNotFoundError: No module named 'intake_xarray'
[14]:
temp_diff = thetao_model - obs.thetao
salt_diff = salt_model - obs.so
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 temp_diff = thetao_model - obs.thetao
2 salt_diff = salt_model - obs.so
NameError: name 'thetao_model' is not defined
Construct a 3D area with land values masked
[15]:
area3d = np.repeat(area.values[np.newaxis, :, :], len(temp_diff.z_l), axis=0)
mask3d = xr.DataArray(area3d, dims=(temp_diff.dims[1:4]), coords= {temp_diff.dims[1]: temp_diff.z_l,
temp_diff.dims[2]: temp_diff.yh,
temp_diff.dims[3]: temp_diff.xh})
area3d_masked = mask3d.where(temp_diff[0,:] == temp_diff[0,:])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 area3d = np.repeat(area.values[np.newaxis, :, :], len(temp_diff.z_l), axis=0)
2 mask3d = xr.DataArray(area3d, dims=(temp_diff.dims[1:4]), coords= {temp_diff.dims[1]: temp_diff.z_l,
3 temp_diff.dims[2]: temp_diff.yh,
4 temp_diff.dims[3]: temp_diff.xh})
5 area3d_masked = mask3d.where(temp_diff[0,:] == temp_diff[0,:])
NameError: name 'area' is not defined
Horizontal Mean difference (model - obs)
[16]:
temp_bias = HorizontalMeanDiff_da(temp_diff,weights=area3d_masked, basins=basin_code)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 temp_bias = HorizontalMeanDiff_da(temp_diff,weights=area3d_masked, basins=basin_code)
NameError: name 'temp_diff' is not defined
[17]:
print('Saving temp_bias...\n')
os.makedirs('ncfiles', exist_ok=True)``
var = 'thetao'
attrs = {'casename': casename,
'description': 'Annual mean bias for '+var,
'obs': args.obs
}
add_global_attrs(temp_bias,attrs)
temp_bias.to_netcdf('ncfiles/'+str(casename)+'_{}_drift.nc'.format(var))
Cell In[17], line 2
os.makedirs('ncfiles', exist_ok=True)``
^
SyntaxError: invalid syntax
Temperature
[18]:
for reg in temp_bias.region:
# remove Nan's
diff_reg = temp_bias.sel(region=reg).dropna('z_l')
if diff_reg.z_l.max() <= 500.0:
splitscale = None
else:
splitscale = [0., -500., -diff_reg.z_l.max()]
ztplot(diff_reg.values, diff_reg.time.values, diff_reg.z_l.values*-1, ignore=np.nan, splitscale=splitscale,
suptitle=casename, contour=True,
title= str(reg.values) + ', Potential Temperature [C], (model - obs)',
extend='both', colormap='dunnePM', autocenter=True, tunits='Year', show=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 for reg in temp_bias.region:
2 # remove Nan's
3 diff_reg = temp_bias.sel(region=reg).dropna('z_l')
4 if diff_reg.z_l.max() <= 500.0:
NameError: name 'temp_bias' is not defined
[19]:
salt_bias = HorizontalMeanDiff_da(salt_diff,weights=area3d_masked, basins=basin_code)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 salt_bias = HorizontalMeanDiff_da(salt_diff,weights=area3d_masked, basins=basin_code)
NameError: name 'salt_diff' is not defined
[20]:
print('Saving salt_bias...\n')
os.makedirs('ncfiles', exist_ok=True)
var = 'so'
attrs = {'casename': casename,
'description': 'Annual mean bias for '+var,
'obs': args.obs
}
add_global_attrs(salt_bias,attrs)
salt_bias.to_netcdf('ncfiles/'+str(casename)+'_{}_drift.nc'.format(var))
Saving salt_bias...
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 5
2 os.makedirs('ncfiles', exist_ok=True)
4 var = 'so'
----> 5 attrs = {'casename': casename,
6 'description': 'Annual mean bias for '+var,
7 'obs': args.obs
8 }
10 add_global_attrs(salt_bias,attrs)
11 salt_bias.to_netcdf('ncfiles/'+str(casename)+'_{}_drift.nc'.format(var))
NameError: name 'casename' is not defined
Salinity
[21]:
for reg in salt_bias.region:
# remove Nan's
diff_reg = salt_bias.sel(region=reg).dropna('z_l')
if diff_reg.z_l.max() <= 500.0:
splitscale = None
else:
splitscale = [0., -500., -diff_reg.z_l.max()]
ztplot(diff_reg.values, diff_reg.time.values, diff_reg.z_l.values*-1, ignore=np.nan, splitscale=splitscale,
suptitle=casename, contour=True,
title= str(reg.values) + ', Salinity [psu], (model - obs)',
extend='both', colormap='dunnePM', autocenter=True, tunits='Year', show=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 1
----> 1 for reg in salt_bias.region:
2 # remove Nan's
3 diff_reg = salt_bias.sel(region=reg).dropna('z_l')
4 if diff_reg.z_l.max() <= 500.0:
NameError: name 'salt_bias' is not defined
Horizontal Mean RMSe (model - obs)
[22]:
# TODO
#temp_rms = HorizontalMeanDiff_da(temp_diff,weights=area3d_masked, basins=basin_code)
Temperature
[23]:
for reg in temp_rms.region:
# remove Nan's
diff_reg = temp_rms.sel(region=reg).dropna('z_l')
if diff_reg.z_l.max() <= 500.0:
splitscale = None
else:
splitscale = [0., -500., -diff_reg.z_l.max()]
ztplot(diff_reg.values, diff_reg.time.values, diff_reg.z_l.values*-1, ignore=np.nan, splitscale=splitscale,
suptitle=dcase._casename, contour=True,
title= str(reg.values) + ', Potential Temperature [C], RMSe',
extend='both', colormap='dunnePM', autocenter=False, tunits='Year', show=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 for reg in temp_rms.region:
2 # remove Nan's
3 diff_reg = temp_rms.sel(region=reg).dropna('z_l')
4 if diff_reg.z_l.max() <= 500.0:
NameError: name 'temp_rms' is not defined
[24]:
# TODO
#salt_rms = HorizontalMeanRmse_da(salt_diff,weights=area3d_masked, basins=basin_code)
Salinity
[25]:
for reg in salt_rms.region:
# remove Nan's
diff_reg = salt_rms.sel(region=reg).dropna('z_l')
if diff_reg.z_l.max() <= 500.0:
splitscale = None
else:
splitscale = [0., -500., -diff_reg.z_l.max()]
ztplot(diff_reg.values, diff_reg.time.values, diff_reg.z_l.values*-1, ignore=np.nan,
splitscale=splitscale, suptitle=dcase._casename, contour=True,
title= str(reg.values) + ', Salinity [psu], RMSe', extend='both',
colormap='dunnePM', autocenter=False, tunits='Year', show=True);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 1
----> 1 for reg in salt_rms.region:
2 # remove Nan's
3 diff_reg = salt_rms.sel(region=reg).dropna('z_l')
4 if diff_reg.z_l.max() <= 500.0:
NameError: name 'salt_rms' is not defined