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---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 7
3 import numpy as np
4 import matplotlib.pyplot as plt
5 import cartopy
6 import cartopy.crs as ccrs
----> 7 import pop_tools
8 from dask.distributed import LocalCluster
9 import s3fs
10
File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.14/site-packages/pop_tools/__init__.py:4
1 """Top-level module for pop_tools"""
----> 4 from pkg_resources import DistributionNotFound, get_distribution
6 from .calc import cfc11sol, cfc12sol
7 from .datasets import DATASETS
ModuleNotFoundError: No module named 'pkg_resources'General setup (see intro notebooks for explanations)¶
Connect to cluster¶
cluster = LocalCluster()
client = cluster.get_client()cluster.scale(20)clientBring 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.01Load 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)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
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