Analysis of Forcing Fields
mom6_tools.MOM6grid returns an object with MOM6 grid data.
mom6_tools.latlon_analysis has a collection of tools used to perform spatial analysis (e.g., time averages and spatial mean).
The goal of this notebook is the following:
serve as an example of how to post-process CESM/MOM6 output;
create time averages of forcing fields;
create time-series of globally-averaged forcing fields;
compare model results vs. observations (TODO);
compare model results vs. another model results (TODO).
[1]:
%load_ext autoreload
%autoreload 2
[2]:
%matplotlib inline
import warnings, yaml, os
warnings.filterwarnings("ignore")
from mom6_tools.MOM6grid import MOM6grid
from mom6_tools.m6toolbox import genBasinMasks, cime_xmlquery
from mom6_tools.latlon_analysis import time_mean_latlon
import matplotlib.pyplot as plt
import numpy as np
Basemap module not found. Some regional plots may not function properly
ERROR 1: PROJ: proj_create_from_database: Open of /glade/work/gmarques/conda-envs/mom6-tools/share/proj failed
[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')
rundir = cime_xmlquery(caseroot, 'RUNDIR')
if DOUT_S:
OUTDIR = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ocn/hist/'
else:
OUTDIR = cime_xmlquery(caseroot, 'RUNDIR')
print('Rundir directory is:', rundir)
print('Casename is:', casename)
Rundir directory is: /glade/derecho/scratch/gmarques/g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079/run
Casename is: g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079
[5]:
# create an empty class object
class args:
pass
args.case_name = casename
args.start_date = ''
args.end_date = ''
args.savefigs = False
args.time_series = True
# set avg dates
avg = diag_config_yml['Avg']
if not args.start_date : args.start_date = avg['start_date']
if not args.end_date : args.end_date = avg['end_date']
# array with 2D forcing variables to be processed
args.variables = ['hfds','wfo', 'prlq']
args.nw = 6 # number of workers
args.static = casename+diag_config_yml['Fnames']['static']
args.native = casename+diag_config_yml['Fnames']['native']
args.geom = casename+diag_config_yml['Fnames']['geom']
args.infile = OUTDIR+args.native
[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.values
except:
area = grd.areacello.values
MOM6 grid successfully loaded...
[7]:
# remove Nan's, otherwise genBasinMasks won't work
depth[np.isnan(depth)] = 0.0
basin_code = genBasinMasks(grd.geolon.values, grd.geolat.values, depth, verbose=False)
[8]:
%matplotlib inline
# plot time averages. If variables is NOT specified, all 2D variables will be plotted.
#time_mean_latlon_plot(args,grd, variables=args.variables)
time_mean_latlon(args,grd,args.variables)
Requesting 6 workers...
/proxy/43571/status
About to plot time-average for Surface ocean heat flux from SW+LW+latent+sensible+masstransfer+frazil+seaice_melt_heat (hfds)...
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
/glade/derecho/scratch/gmarques/tmp/ipykernel_21702/3597705672.py in ?()
1 get_ipython().run_line_magic('matplotlib', 'inline')
2 # plot time averages. If variables is NOT specified, all 2D variables will be plotted.
3 #time_mean_latlon_plot(args,grd, variables=args.variables)
----> 4 time_mean_latlon(args,grd,args.variables)
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/mom6_tools/latlon_analysis.py in ?(args, grd, variables)
171 title=r'%s, [%s] averaged over years %s-%s'%(var,units,args.start_date,args.end_date),
172 extend='both',
173 save=filename)
174 else:
--> 175 m6plot.xyplot( data , grd.geolon, grd.geolat, area=area,
176 suptitle=args.case_name,
177 title=r'%s, [%s] averaged over years %s-%s'%(var,units,args.start_date,args.end_date),
178 extend='both',
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/mom6_tools/m6plot.py in ?(field, x, y, area, xlabel, xunits, ylabel, yunits, title, suptitle, clim, colormap, extend, centerlabels, nbins, landcolor, axis, add_cbar, aspect, resolution, sigma, annotate, ignore, save, debug, show, interactive, logscale)
600
601 # Diagnose statistics
602 if ignore is not None: maskedField = numpy.ma.masked_array(field, mask=[field==ignore])
603 else: maskedField = field.copy()
--> 604 sMin, sMax, sMean, sStd, sRMS = myStats(maskedField, area, debug=debug)
605 xLims = boundaryStats(xCoord)
606 yLims = boundaryStats(yCoord)
607
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/mom6_tools/m6plot.py in ?(s, area, debug)
1307 if debug: print('myStats: max(s) =',sMax)
1308 if area is None: return sMin, sMax, None, None, None
1309 weight = area.copy()
1310 if debug: print('myStats: sum(area) =',numpy.ma.sum(weight))
-> 1311 if not numpy.ma.getmask(s).any()==numpy.ma.nomask: weight[s.mask] = 0.
1312 sumArea = numpy.ma.sum(weight)
1313 if debug: print('myStats: sum(area) =',sumArea,'after masking')
1314 if debug: print('myStats: sum(s) =',numpy.ma.sum(s))
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, key, value)
881 else:
882 # Coordinates in key, value and self[key] should be consistent.
883 # TODO Coordinate consistency in key is checked here, but it
884 # causes unnecessary indexing. It should be optimized.
--> 885 obj = self[key]
886 if isinstance(value, DataArray):
887 assert_coordinate_consistent(value, obj.coords.variables)
888 value = value.variable
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, key)
872 if isinstance(key, str):
873 return self._getitem_coord(key)
874 else:
875 # xarray-style array indexing
--> 876 return self.isel(indexers=self._item_key_to_dict(key))
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/xarray/core/dataarray.py in ?(self, indexers, drop, missing_dims, **indexers_kwargs)
1497
1498 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
1499
1500 if any(is_fancy_indexer(idx) for idx in indexers.values()):
-> 1501 ds = self._to_temp_dataset()._isel_fancy(
1502 indexers, drop=drop, missing_dims=missing_dims
1503 )
1504 return self._from_temp_dataset(ds)
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/xarray/core/dataset.py in ?(self, indexers, drop, missing_dims)
3012 *,
3013 drop: bool,
3014 missing_dims: ErrorOptionsWithWarn = "raise",
3015 ) -> Self:
-> 3016 valid_indexers = dict(self._validate_indexers(indexers, missing_dims))
3017
3018 variables: dict[Hashable, Variable] = {}
3019 indexes, index_variables = isel_indexes(self.xindexes, valid_indexers)
/glade/work/gmarques/conda-envs/mom6-tools/lib/python3.11/site-packages/xarray/core/dataset.py in ?(self, indexers, missing_dims)
2784 elif isinstance(index, CFTimeIndex):
2785 v = _parse_array_of_cftime_strings(v, index.date_type)
2786
2787 if v.ndim > 1:
-> 2788 raise IndexError(
2789 "Unlabeled multi-dimensional array cannot be "
2790 f"used for indexing: {k}"
2791 )
IndexError: Unlabeled multi-dimensional array cannot be used for indexing: yh
[ ]: