Overview¶
The availability of several macronutrients controls production in most of the ocean: nitrate, phosphate, and silicate. Here we take a look at maps and depth profiles of these nutrients, and compare them to an observational dataset.
- General setup
- Subsetting
- Processing - means in time and space
- Compare to World Ocean Atlas data
- Make depth profiles
Prerequisites¶
| Concepts | Importance | Notes |
|---|---|---|
| Matplotlib | Necessary | |
| Intro to Cartopy | Necessary | |
| Dask Cookbook | Helpful | |
| Intro to Xarray | Helpful |
- Time to learn: 30 min
Imports¶
import xarray as xr
import glob
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import pop_tools
from dask.distributed import LocalCluster
import s3fs
from module import adjust_pop_grid/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/pop_tools/__init__.py:4: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
from pkg_resources import DistributionNotFound, get_distribution
General setup (see intro notebooks for explanations)¶
Connect to cluster¶
cluster = LocalCluster()
client = cluster.get_client()/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/node.py:187: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 41067 instead
warnings.warn(
cluster.scale(20)clientLoading...
Bring in POP grid utilities¶
ds_grid = pop_tools.get_grid('POP_gx1v7')
lons = ds_grid.TLONG
lats = ds_grid.TLAT
depths = ds_grid.z_t * 0.01Downloading file 'inputdata/ocn/pop/gx1v7/grid/horiz_grid_20010402.ieeer8' from 'https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/ocn/pop/gx1v7/grid/horiz_grid_20010402.ieeer8' to '/home/runner/.pop_tools'.
/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/urllib3/connectionpool.py:1097: InsecureRequestWarning: Unverified HTTPS request is being made to host 'svn-ccsm-inputdata.cgd.ucar.edu'. Adding certificate verification is strongly advised. See: https://urllib3.readthedocs.io/en/latest/advanced-usage.html#tls-warnings
warnings.warn(
Downloading file 'inputdata/ocn/pop/gx1v7/grid/region_mask_20151008.ieeei4' from 'https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/ocn/pop/gx1v7/grid/region_mask_20151008.ieeei4' to '/home/runner/.pop_tools'.
/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/urllib3/connectionpool.py:1097: InsecureRequestWarning: Unverified HTTPS request is being made to host 'svn-ccsm-inputdata.cgd.ucar.edu'. Adding certificate verification is strongly advised. See: https://urllib3.readthedocs.io/en/latest/advanced-usage.html#tls-warnings
warnings.warn(
Load the data¶
jetstream_url = 'https://js2.jetstream-cloud.org:8001/'
s3 = s3fs.S3FileSystem(anon=True, client_kwargs=dict(endpoint_url=jetstream_url))
# Generate a list of all files in CESM folder
s3path = 's3://pythia/ocean-bgc/cesm/g.e22.GOMIPECOIAF_JRA-1p4-2018.TL319_g17.4p2z.002branch/ocn/proc/tseries/month_1/*'
remote_files = s3.glob(s3path)
s3.invalidate_cache()
# Open all files from folder
fileset = [s3.open(file) for file in remote_files]
# Open with xarray
ds = xr.open_mfdataset(fileset, data_vars="minimal", coords='minimal', compat="override", parallel=True,
drop_variables=["transport_components", "transport_regions", 'moc_components'], decode_times=True)Fetching long content....
2025-09-05 01:45:13,779 - distributed.protocol.core - CRITICAL - Failed to deserialize
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/core.py", line 175, in loads
return msgpack.loads(
~~~~~~~~~~~~~^
frames[0], object_hook=_decode_default, use_list=False, **msgpack_opts
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "msgpack/_unpacker.pyx", line 194, in msgpack._cmsgpack.unpackb
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/core.py", line 159, in _decode_default
return merge_and_deserialize(
sub_header, sub_frames, deserializers=deserializers
)
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/contextlib.py", line 85, in inner
return func(*args, **kwds)
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 525, in merge_and_deserialize
return deserialize(header, merged_frames, deserializers=deserializers)
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 452, in deserialize
return loads(header, frames)
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 195, in serialization_error_loads
raise TypeError(msg)
TypeError: Could not serialize object of type method
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/pickle.py", line 63, in dumps
result = pickle.dumps(x, **dump_kwargs)
TypeError: cannot pickle 'module' object
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/pickle.py", line 68, in dumps
pickler.dump(x)
~~~~~~~~~~~~^^^
TypeError: cannot pickle 'module' object
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 366, in serialize
header, frames = dumps(x, context=context) if wants_context else dumps(x)
~~~~~^^^^^^^^^^^^^^^^^^^^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 78, in pickle_dumps
frames[0] = pickle.dumps(
~~~~~~~~~~~~^
x,
^^
buffer_callback=buffer_callback,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
protocol=context.get("pickle-protocol", None) if context else None,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/pickle.py", line 80, in dumps
result = cloudpickle.dumps(x, **dump_kwargs)
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/cloudpickle/cloudpickle.py", line 1537, in dumps
cp.dump(obj)
~~~~~~~^^^^^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/cloudpickle/cloudpickle.py", line 1303, in dump
return super().dump(obj)
~~~~~~~~~~~~^^^^^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/h5py/_hl/base.py", line 369, in __getnewargs__
raise TypeError("h5py objects cannot be pickled")
TypeError: h5py objects cannot be pickled
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[6], line 13
10 fileset = [s3.open(file) for file in remote_files]
12 # Open with xarray
---> 13 ds = xr.open_mfdataset(fileset, data_vars="minimal", coords='minimal', compat="override", parallel=True,
14 drop_variables=["transport_components", "transport_regions", 'moc_components'], decode_times=True)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/xarray/backends/api.py:1812, in open_mfdataset(paths, chunks, concat_dim, compat, preprocess, engine, data_vars, coords, combine, parallel, join, attrs_file, combine_attrs, errors, **kwargs)
1807 datasets = [preprocess(ds) for ds in datasets]
1809 if parallel:
1810 # calling compute here will return the datasets/file_objs lists,
1811 # the underlying datasets will still be stored as dask arrays
-> 1812 datasets, closers = dask.compute(datasets, closers)
1814 # Combine all datasets, closing them in case of a ValueError
1815 try:
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/dask/base.py:681, in compute(traverse, optimize_graph, scheduler, get, *args, **kwargs)
678 expr = expr.optimize()
679 keys = list(flatten(expr.__dask_keys__()))
--> 681 results = schedule(expr, keys, **kwargs)
683 return repack(results)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/utils_comm.py:416, in retry_operation(coro, operation, *args, **kwargs)
410 retry_delay_min = parse_timedelta(
411 dask.config.get("distributed.comm.retry.delay.min"), default="s"
412 )
413 retry_delay_max = parse_timedelta(
414 dask.config.get("distributed.comm.retry.delay.max"), default="s"
415 )
--> 416 return await retry(
417 partial(coro, *args, **kwargs),
418 count=retry_count,
419 delay_min=retry_delay_min,
420 delay_max=retry_delay_max,
421 operation=operation,
422 )
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/utils_comm.py:395, in retry(coro, count, delay_min, delay_max, jitter_fraction, retry_on_exceptions, operation)
393 delay *= 1 + random.random() * jitter_fraction
394 await asyncio.sleep(delay)
--> 395 return await coro()
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/core.py:1259, in PooledRPCCall.__getattr__.<locals>.send_recv_from_rpc(**kwargs)
1257 prev_name, comm.name = comm.name, "ConnectionPool." + key
1258 try:
-> 1259 return await send_recv(comm=comm, op=key, **kwargs)
1260 finally:
1261 self.pool.reuse(self.addr, comm)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/core.py:1018, in send_recv(comm, reply, serializers, deserializers, **kwargs)
1016 await comm.write(msg, serializers=serializers, on_error="raise")
1017 if reply:
-> 1018 response = await comm.read(deserializers=deserializers)
1019 else:
1020 response = None
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/comm/tcp.py:248, in TCP.read(self, deserializers)
246 else:
247 try:
--> 248 msg = await from_frames(
249 frames,
250 deserialize=self.deserialize,
251 deserializers=deserializers,
252 allow_offload=self.allow_offload,
253 )
254 except EOFError:
255 # Frames possibly garbled or truncated by communication error
256 self.abort()
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/comm/utils.py:78, in from_frames(frames, deserialize, deserializers, allow_offload)
76 res = await offload(_from_frames)
77 else:
---> 78 res = _from_frames()
80 return res
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/comm/utils.py:61, in from_frames.<locals>._from_frames()
59 def _from_frames():
60 try:
---> 61 return protocol.loads(
62 frames, deserialize=deserialize, deserializers=deserializers
63 )
64 except EOFError:
65 if size > 1000:
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/core.py:175, in loads(frames, deserialize, deserializers)
172 return pickle.loads(sub_header["pickled-obj"], buffers=sub_frames)
173 return msgpack_decode_default(obj)
--> 175 return msgpack.loads(
176 frames[0], object_hook=_decode_default, use_list=False, **msgpack_opts
177 )
179 except Exception:
180 logger.critical("Failed to deserialize", exc_info=True)
File msgpack/_unpacker.pyx:194, in msgpack._cmsgpack.unpackb()
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/core.py:159, in loads.<locals>._decode_default(obj)
157 if "compression" in sub_header:
158 sub_frames = decompress(sub_header, sub_frames)
--> 159 return merge_and_deserialize(
160 sub_header, sub_frames, deserializers=deserializers
161 )
162 else:
163 return Serialized(sub_header, sub_frames)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/contextlib.py:85, in ContextDecorator.__call__.<locals>.inner(*args, **kwds)
82 @wraps(func)
83 def inner(*args, **kwds):
84 with self._recreate_cm():
---> 85 return func(*args, **kwds)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py:525, in merge_and_deserialize(header, frames, deserializers)
521 merged = host_array_from_buffers(subframes)
523 merged_frames.append(merged)
--> 525 return deserialize(header, merged_frames, deserializers=deserializers)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py:452, in deserialize(header, frames, deserializers)
447 raise TypeError(
448 "Data serialized with %s but only able to deserialize "
449 "data with %s" % (name, str(list(deserializers)))
450 )
451 dumps, loads, wants_context = families[name]
--> 452 return loads(header, frames)
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py:195, in serialization_error_loads(header, frames)
193 def serialization_error_loads(header, frames):
194 msg = "\n".join([codecs.decode(frame, "utf8") for frame in frames])
--> 195 raise TypeError(msg)
TypeError: Could not serialize object of type method
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/pickle.py", line 63, in dumps
result = pickle.dumps(x, **dump_kwargs)
TypeError: cannot pickle 'module' object
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/pickle.py", line 68, in dumps
pickler.dump(x)
~~~~~~~~~~~~^^^
TypeError: cannot pickle 'module' object
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 366, in serialize
header, frames = dumps(x, context=context) if wants_context else dumps(x)
~~~~~^^^^^^^^^^^^^^^^^^^^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/serialize.py", line 78, in pickle_dumps
frames[0] = pickle.dumps(
~~~~~~~~~~~~^
x,
^^
buffer_callback=buffer_callback,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
protocol=context.get("pickle-protocol", None) if context else None,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/distributed/protocol/pickle.py", line 80, in dumps
result = cloudpickle.dumps(x, **dump_kwargs)
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/cloudpickle/cloudpickle.py", line 1537, in dumps
cp.dump(obj)
~~~~~~~^^^^^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/cloudpickle/cloudpickle.py", line 1303, in dump
return super().dump(obj)
~~~~~~~~~~~~^^^^^
File "/home/runner/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.13/site-packages/h5py/_hl/base.py", line 369, in __getnewargs__
raise TypeError("h5py objects cannot be pickled")
TypeError: h5py objects cannot be pickled
Subsetting¶
Make our dataset smaller so it has just a couple of macronutrient variables we’re interested in.
variables =['PO4','NO3','SiO3']keep_vars=['z_t','z_t_150m','dz','time_bound','time','TAREA','TLAT','TLONG'] + variables
ds = ds.drop_vars([v for v in ds.variables if v not in keep_vars])Let’s take a quick look at nitrate to make sure that things look okay...
ds.NO3.isel(time=0,z_t=0).plot(cmap="viridis")Transforming from monthly to annual data¶
We can’t just use xarray’s regular mean() function because months have different numbers of days in them, so we have to weight by that to ensure the annual mean is accurate. See this ESDS blog post for a more detailed explanation with examples!
def year_mean(ds):
"""
Properly convert monthly data to annual means, taking into account month lengths.
Source: https://ncar.github.io/esds/posts/2021/yearly-averages-xarray/
"""
# Make a DataArray with the number of days in each month, size = len(time)
month_length = ds.time.dt.days_in_month
# Calculate the weights by grouping by 'time.year'
weights = (
month_length.groupby("time.year") / month_length.groupby("time.year").sum()
)
# Test that the sum of the year for each season is 1.0
np.testing.assert_allclose(weights.groupby("time.year").sum().values, np.ones((len(ds.groupby("time.year")), )))
# Calculate the weighted average
return (ds * weights).groupby("time.year").sum(dim="time")ds_annual = year_mean(ds)
ds_annualNote that our time coordinate is now called year instead, and has only years now. We can select specific years to plot:
ds_annual['NO3'].sel(year=2010).isel(z_t=0).plot()Let’s make a nicer-looking map¶
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(1,1,1, projection=ccrs.Robinson(central_longitude=305.0))
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)
ax.set_title('CESM surface NO$_3$', fontsize=10)
lon, lat, field = adjust_pop_grid(lons, lats, ds_annual.NO3.sel(year=2010).isel(z_t=0))
pc1=ax.pcolormesh(lon, lat,field, vmin=0, vmax=20, cmap='Greens',
transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(pc1, ax=ax,extend='max',label='NO$_3$ (mmol m$^{-3}$)')Compare long-term mean to World Ocean Atlas 2018¶
- About the World Ocean Atlas
- Data access to WOA18
We’ve already regridded the WOA data to be on the same grid as the CESM data, so we don’t need to worry about that step. However, if you wanted to compare to a dataset that’s on a different grid, you’d need to go through the regridding process, which is beyond the scope of this cookbook.
This dataset has also already had a time mean taken, so there’s no time coordinate.
You might notice that there are three coordinates: z_t, z_w, and z_w_bot. Each of these are different versions of the same vertical coordinate - z_t represents the midpoint of a depth layer, z_w the top, and z_w_bot the bottom. We use z_t in this demonstration.
woa_file_path = 's3://pythia/ocean-bgc/obs/WOA2018_POPgrid.nc'
woa_file = s3.open(woa_file_path)
ds_woa = xr.load_dataset(woa_file, decode_times=False, decode_coords=False)
ds_woa['z_t'] = ds.z_t
ds_woaNow that we’re doing more involved calculations, we’re going to just take a mean over a couple years (2010-2011) to make the computational load a bit lighter. For a more accurate analysis, we’d want to include more years than this.
ds_annual_subset = ds_annual.sel(year=[2010,2011])
ds_mean = ds_annual_subset.mean("year").compute()NO3_diff = ds_mean.NO3 - ds_woa.NO3
PO4_diff = ds_mean.PO4 - ds_woa.PO4
SiO3_diff = ds_mean.SiO3 - ds_woa.SiO3Surface comparison¶
We choose to set up a dictionary with some parameters for each plot we want to make, to cut down on repetition in the actual plotting code block. This could be condensed even further, but there’s a tradeoff between conciseness and readability! We specify the variables we want to plot (in this case different nutrients) and things like the colormaps and normalization. In addition to plotting each nutrient from the modeled data and observations, we also plot the bias, which is the difference between the two datasets. This helps us see how the model differs from observations.
ds_dict_surf = {'CESMNO3': {'title': 'CESM surface NO$_3$',
'label': 'NO$_3$ (mmol m$^{-3}$)',
'cmap': 'Greens',
'vmin': 0, 'vmax': 20,
'ds': ds_mean.NO3},
'WOANO3': {'title': 'WOA surface NO$_3$',
'label': 'NO$_3$ (mmol m$^{-3}$)',
'cmap': 'Greens',
'vmin': 0, 'vmax': 20,
'ds': ds_woa.NO3},
'DIFFNO3': {'title': 'Surface NO$_3$ model bias',
'label': 'NO$_3$ bias (mmol m$^{-3}$)',
'cmap': 'bwr',
'vmin': -10, 'vmax': 10,
'ds': ds_mean.NO3 - ds_woa.NO3},
'CESMPO4': {'title': 'CESM surface PO$_4$',
'label': 'PO$_4$ (mmol m$^{-3}$)',
'cmap': 'Oranges',
'vmin': 0, 'vmax': 2,
'ds': ds_mean.PO4},
'WOAPO4': {'title': 'WOA surface PO$_4$',
'label': 'PO$_4$ (mmol m$^{-3}$)',
'cmap': 'Oranges',
'vmin': 0, 'vmax': 2,
'ds': ds_woa.PO4},
'DIFFPO4': {'title': 'Surface PO$_4$ model bias',
'label': 'PO$_4$ bias (mmol m$^{-3}$)',
'cmap': 'bwr',
'vmin': -1, 'vmax': 1,
'ds': ds_mean.PO4 - ds_woa.PO4},
'CESMSiO3': {'title': 'CESM surface SiO$_3$',
'label': 'SiO$_3$ (mmol m$^{-3}$)',
'cmap': 'Blues',
'vmin': 0, 'vmax': 30,
'ds': ds_mean.SiO3},
'WOASiO3': {'title': 'WOA surface SiO$_3$',
'label': 'SiO$_3$ (mmol m$^{-3}$)',
'cmap': 'Blues',
'vmin': 0, 'vmax': 30,
'ds': ds_woa.SiO3},
'DIFFSiO3': {'title': 'Surface SiO$_3$ model bias',
'label': 'SiO$_3$ bias (mmol m$^{-3}$)',
'cmap': 'bwr',
'vmin': -15, 'vmax': 15,
'ds': ds_mean.SiO3 - ds_woa.SiO3}
}
Here we pull from the above dictionary to actually make the plots.
fig = plt.figure(figsize=(18,10))
plot_count = 1
for key, item in ds_dict_surf.items():
ax = fig.add_subplot(3,3,plot_count, projection=ccrs.Robinson(central_longitude=305.0))
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)
ax.set_title(item['title'], fontsize=10)
lon, lat, field = adjust_pop_grid(lons, lats, item['ds'].isel(z_t=0))
pc=ax.pcolormesh(lon, lat,field, vmin=item['vmin'], vmax=item['vmax'], cmap=item['cmap'],
transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(pc, ax=ax,label=item['label'])
plot_count += 1
Comparison at 100m¶
Similar to above, but at a depth of 100m rather than at the surface.
ds_dict_100m = {'CESMNO3': {'title': 'CESM 100m NO$_3$',
'label': 'NO$_3$ (mmol m$^{-3}$)',
'cmap': 'Greens',
'vmin': 0, 'vmax': 20,
'ds': ds_mean.NO3},
'WOANO3': {'title': 'WOA 100m NO$_3$',
'label': 'NO$_3$ (mmol m$^{-3}$)',
'cmap': 'Greens',
'vmin': 0, 'vmax': 20,
'ds': ds_woa.NO3},
'DIFFNO3': {'title': '100m NO$_3$ model bias',
'label': 'NO$_3$ bias (mmol m$^{-3}$)',
'cmap': 'bwr',
'vmin': -10, 'vmax': 10,
'ds': ds_mean.NO3 - ds_woa.NO3},
'CESMPO4': {'title': 'CESM 100m PO$_4$',
'label': 'PO$_4$ (mmol m$^{-3}$)',
'cmap': 'Oranges',
'vmin': 0, 'vmax': 2,
'ds': ds_mean.PO4},
'WOAPO4': {'title': 'WOA 100m PO$_4$',
'label': 'PO$_4$ (mmol m$^{-3}$)',
'cmap': 'Oranges',
'vmin': 0, 'vmax': 2,
'ds': ds_woa.PO4},
'DIFFPO4': {'title': '100m PO$_4$ model bias',
'label': 'PO$_4$ bias (mmol m$^{-3}$)',
'cmap': 'bwr',
'vmin': -1, 'vmax': 1,
'ds': ds_mean.PO4 - ds_woa.PO4},
'CESMSiO3': {'title': 'CESM 100m SiO$_3$',
'label': 'SiO$_3$ (mmol m$^{-3}$)',
'cmap': 'Blues',
'vmin': 0, 'vmax': 30,
'ds': ds_mean.SiO3},
'WOASiO3': {'title': 'WOA 100m SiO$_3$',
'label': 'SiO$_3$ (mmol m$^{-3}$)',
'cmap': 'Blues',
'vmin': 0, 'vmax': 30,
'ds': ds_woa.SiO3},
'DIFFSiO3': {'title': '100m SiO$_3$ model bias',
'label': 'SiO$_3$ bias (mmol m$^{-3}$)',
'cmap': 'bwr',
'vmin': -15, 'vmax': 15,
'ds': ds_mean.SiO3 - ds_woa.SiO3}
}
fig = plt.figure(figsize=(18,10))
plot_count = 1
for key, item in ds_dict_100m.items():
ax = fig.add_subplot(3,3,plot_count, projection=ccrs.Robinson(central_longitude=305.0))
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)
ax.set_title(item['title'], fontsize=10)
lon, lat, field = adjust_pop_grid(lons, lats, item['ds'].isel(z_t=10))
pc=ax.pcolormesh(lon, lat,field, vmin=item['vmin'], vmax=item['vmax'], cmap=item['cmap'],
transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(pc, ax=ax,label=item['label'])
plot_count += 1Global mean macronutrient profiles¶
Let’s write a function to take a global mean of the variables we’re interested in, so that we can look at some depth profiles rather than maps. Also remember that we already took a mean over the whole time range (and the WOA dataset already had this mean taken), so this is a mean in time as well. Like the above maps, we also plot a bias panel to directly compare the difference between the datasets.
def global_mean(ds, ds_grid, compute_vars, include_ms=False):
"""
Compute the global mean on a POP dataset.
Return computed quantity in conventional units.
"""
dso = xr.Dataset({v: ds_grid[v] for v in ['z_t']})
for var in compute_vars:
area_depth = np.full([384,320,60],np.nan)
var_profile = np.full([60],np.nan)
for z in np.arange(0,60,1):
if include_ms: # marginal seas
area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.KMT > 0).where(ds[var].isel(z_t=z) > 0)
else:
area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)
area_depth = xr.DataArray(area_depth,dims=('nlat','nlon','z_t'))
for z in np.arange(0,60,1):
var_profile[z] = (ds[var].isel(z_t=z) * area_depth.isel(z_t=z)).sum(dim=('nlon','nlat')) / area_depth.isel(z_t=z).sum(dim=('nlon','nlat'))
dso[var] = var_profile
return dsods_glb = global_mean(ds_mean, ds_grid, ['NO3','PO4','SiO3']).compute()ds_glb_woa = global_mean(ds_woa, ds_grid, ['NO3','PO4','SiO3']).compute()Rather than setting up a dictionary of parameters, here we choose to make the plots inline since there aren’t as many.
fig = plt.figure(figsize=(6,10))
plt.suptitle('Global mean macronutrient profiles', fontsize=14)
### Row 1 - NO3
ax = fig.add_subplot(3,2,1)
ax.set_title('Global mean NO$_3$')
ax.plot(ds_glb_woa['NO3'].values, depths, label='WOA', linewidth=3, color='lightgreen')
ax.plot(ds_glb['NO3'].values, depths, label='CESM', linewidth=3, color='green')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='NO$_3$ (mmol m$^{-3}$)')
plt.gca().invert_yaxis()
# Bias plot
ax = fig.add_subplot(3,2,2)
ax.plot(ds_glb['NO3'].values - ds_glb_woa['NO3'].values, depths, label='bias', linewidth=3, color='black')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='NO$_3$ bias (mmol m$^{-3}$)')
ax.axvline(x=0,color='black',linestyle='--',linewidth=0.5)
plt.gca().invert_yaxis()
### Row 2 - PO4
ax = fig.add_subplot(3,2,3)
ax.set_title('Global mean PO$_4$')
ax.plot(ds_glb_woa['PO4'].values, depths, label='WOA', linewidth=3, color='peachpuff')
ax.plot(ds_glb['PO4'].values, depths, label='CESM', linewidth=3, color='orange')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='PO$_4$ (mmol m$^{-3}$)')
plt.gca().invert_yaxis()
# Bias plot
ax = fig.add_subplot(3,2,4)
ax.plot(ds_glb['PO4'].values - ds_glb_woa['PO4'].values, depths, label='bias', linewidth=3, color='black')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='PO$_4$ bias (mmol m$^{-3}$)')
ax.axvline(x=0,color='black',linestyle='--',linewidth=0.5)
plt.gca().invert_yaxis()
### Row 3 - SiO3
ax = fig.add_subplot(3,2,5)
ax.set_title('Global mean SiO$_3$')
ax.plot(ds_glb_woa['SiO3'].values, depths, label='WOA', linewidth=3, color='lightblue')
ax.plot(ds_glb['SiO3'].values, depths, label='CESM', linewidth=3, color='blue')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='SiO$_3$ (mmol m$^{-3}$)')
plt.gca().invert_yaxis()
# Bias plot
ax = fig.add_subplot(3,2,6)
ax.plot(ds_glb['SiO3'].values - ds_glb_woa['SiO3'].values, depths, label='bias', linewidth=3, color='black')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='SiO$_3$ bias (mmol m$^{-3}$)')
ax.axvline(x=0,color='black',linestyle='--',linewidth=0.5)
plt.gca().invert_yaxis()
fig.tight_layout()
And close the Dask cluster we spun up at the beginning.
cluster.close()Summary¶
You’ve learned how to plot and evaluate the distribution of some key ocean nutrients in CESM output.
Resources and references¶
- (2013). In Ocean Biogeochemical Dynamics (pp. 102–172). Princeton University Press. 10.2307/j.ctt3fgxqx.7