<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Finding forced trends in oceanic oxygen1.1&nbsp;&nbsp;Learning Objectives1.2&nbsp;&nbsp;Create Dask Cluster1.3&nbsp;&nbsp;Open an intake-esm catalog for CESM1-LE data stored on Glade and Tape1.4&nbsp;&nbsp;Load oxygen Data From Control Run into an Xarray Dataset and Compute its Annual Mean1.5&nbsp;&nbsp;Load oxygen Data From 20C, RCP85 Runs into an Xarray Dataset and Compute its Annual Mean1.6&nbsp;&nbsp;Compute Linear Trend1.7&nbsp;&nbsp;Plot All Trends - Total1.8&nbsp;&nbsp;Plot All Trends - Internal1.9&nbsp;&nbsp;Going Further  Finding forced trends in oceanic oxygen  Learning Objectives-<jupyter_code>%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
from ncar_jobqueue import NCARCluster
from dask.distributed import Client
import esmlab
import xarray as xr
import numpy as np
import dask
import util<jupyter_output><empty_output><jupyter_text>Create Dask Cluster<jupyter_code>cluster = NCARCluster(memory="100GB", cores=4, processes=1)
cluster.adapt(minimum=3, maximum=12,
              wait_count=60)  # Adaptively scale between 3 and 12 dask workers
client = Client(cluster)  # Connect this local process to remote workers
cluster<jupyter_output><empty_output><jupyter_text>(👆don't forget to look at the dashboard!)  Open an intake-esm catalog for CESM1-LE data stored on Glade and Tape<jupyter_code>import intake
from intake_esm import config
with config.set({'database-directory': '/glade/u/home/abanihi/.intake_esm/collections'}):
    col = intake.open_esm_metadatastore(collection_name='CESM1-LE')

col.df.head()
keep_vars = ['O2', 'z_t', 'KMT', 'TAREA', 'TLONG',
             'TLAT', 'time', 'time_bound', 'member_id']
yr0, yr1 = 2006, 2055
time_slice = slice(f'{yr0:04d}', f'{yr1:04d}')
ctrl_time_slice = slice(f'{yr0-1448:04d}', f'{yr1-1448:04d}')<jupyter_output><empty_output><jupyter_text>Load oxygen Data From Control Run into an Xarray Dataset and Compute its Annual Mean<jupyter_code>%%time
dsets = col.search(experiment=['CTRL'], variable=['O2'])\
        .to_xarray(decode_times=False, chunks={'time': 240})
_, dd = dsets.popitem()
ctrl = util.clean_ds(dd, keep_vars)
ctrl = util.sel_time(ctrl, ctrl_time_slice)
ctrl = esmlab.resample(ctrl, freq='ann')
print(ctrl)<jupyter_output><empty_output><jupyter_text>Load oxygen Data From 20C, RCP85 Runs into an Xarray Dataset and Compute its Annual Mean<jupyter_code>%%time
dsets = col.search(experiment=['20C', 'RCP85'], variable=['O2'])\
        .to_xarray(decode_times=False, chunks={'time': 240})
print(dsets.keys())
ds_20c = dsets['ocn.20C.pop.h']
ds_rcp85 = dsets['ocn.RCP85.pop.h']
dd = xr.concat([ds_20c, ds_rcp85], dim='time')
ds = util.clean_ds(dd, keep_vars)
ds = util.sel_time(ds, time_slice)
ds = esmlab.resample(ds, freq='ann')
print(ds)<jupyter_output><empty_output><jupyter_text>Compute Linear Trend<jupyter_code>def linear_trend(da, dim='time'):
    da_chunk = da.chunk({dim: -1})
    trend = xr.apply_ufunc(calc_slope,
                           da_chunk,
                           vectorize=True,
                           input_core_dims=[[dim]],
                           output_core_dims=[[]],
                           output_dtypes=[np.float],
                           dask='parallelized')
    return trend


def calc_slope(y):
    """ufunc to be used by linear_trend"""
    x = np.arange(len(y))
    return np.polyfit(x, y, 1)[0]
%%time
ctrl_trend = linear_trend(ctrl.O2.chunk({
    'nlat': 30,
    'nlon': 30,
    'time': -1
})).load()
ctrl_trend = ctrl_trend * len(ctrl.time)
ctrl_trend.attrs['units'] = f'mmol m$^{-3}$ ({len(ctrl.time)} yr)$^{-1}$'
ctrl_trend.plot()
%%time
ds_trend = linear_trend(ds.O2.chunk({'nlat': 30, 'nlon': 30, 'time': -1}))\
                       .load()
ds_trend = ds_trend * len(ds.time)
ds_trend.attrs['units'] = f'mmol m$^{{-3}}$ ({len(ds.time)} yr)$^{{-1}}$'
ds_trend.isel(member_id=3).plot()
## 
npac_trend = ds_trend.where((10 < ds.TLAT) & (ds.TLAT < 65) & (120 < ds.TLONG)
                            & (ds.TLONG < 260))
npac_trend = esmlab.weighted_mean(npac_trend,
                                  dim=('nlat', 'nlon'),
                                  weights=ds.TAREA).load()
npac_trend
## 
member_id_pick = [
    npac_trend.where(npac_trend == npac_trend.max(),
                     drop=True).member_id.values.astype('int')[0],
    npac_trend.where(npac_trend == npac_trend.min(),
                     drop=True).member_id.values.astype('int')[0]
]
member_id_pick
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import cmocean
import cartopy.crs as ccrs
import cartopy
levels = np.arange(-100, 45, 5)
norm = colors.DivergingNorm(vmin=levels[0], vmax=levels[-1], vcenter=0.)

dsp = util.pop_add_cyclic(
    xr.Dataset({
        'trend': ds_trend,
        'TLAT': ds.TLAT[0].compute(),
        'TLONG': ds.TLONG[0].compute()
    }))

forced = dsp.trend.mean('member_id').compute()
internal = [(dsp.trend.sel(member_id=i) - forced).compute()
            for i in member_id_pick]

extent = [120, 260, 10, 65]
prj = ccrs.Mercator(central_longitude=np.mean(extent[0:2]),
                    min_latitude=extent[2],
                    max_latitude=extent[3])
def one_plot(da, lines=True):
    # filled contours
    cf = ax.contourf(dsp.TLONG,
                     dsp.TLAT,
                     da,
                     levels=levels,
                     norm=norm,
                     cmap=cmocean.cm.curl_r,
                     extend='both',
                     transform=ccrs.PlateCarree())

    # contour lines
    cs = ax.contour(dsp.TLONG,
                    dsp.TLAT,
                    da,
                    colors='k',
                    levels=levels,
                    linewidths=0.5,
                    transform=ccrs.PlateCarree())

    if lines:
        # add contour labels
        lb = plt.clabel(cs, fontsize=6, inline=True, fmt='%r')

    # land
    land = ax.add_feature(
        cartopy.feature.NaturalEarthFeature('physical',
                                            'land',
                                            '110m',
                                            facecolor='lightgray'))

    ax.coastlines(linewidth=0.5)

    return cf
fig = plt.figure(figsize=(12, 6))

axs = []

# plot total
ax = fig.add_subplot(2, 3, 1, projection=prj)
ax.set_extent(extent)
one_plot(dsp.trend.sel(member_id=member_id_pick[0]))
ax.text(235., 60, f'{member_id_pick[0]:03d}', transform=ccrs.PlateCarree())
ax.set_title('Total')
axs.append(ax)

ax = fig.add_subplot(2, 3, 4, projection=prj)
ax.set_extent(extent)
one_plot(dsp.trend.sel(member_id=member_id_pick[1]))
ax.text(235., 60, f'{member_id_pick[1]:03d}', transform=ccrs.PlateCarree())
axs.append(ax)

# plot internal variability
ax = fig.add_subplot(2, 3, 2, projection=prj)
ax.set_extent(extent)
one_plot(internal[0])
ax.text(235., 60, f'{member_id_pick[0]:03d}', transform=ccrs.PlateCarree())
ax.set_title('Internal')
axs.append(ax)

ax = fig.add_subplot(2, 3, 5, projection=prj)
ax.set_extent(extent)
one_plot(internal[1])
ax.text(235., 60, f'{member_id_pick[1]:03d}', transform=ccrs.PlateCarree())
axs.append(ax)

# plot forced
ax = fig.add_subplot(2, 3, 3, projection=prj)
ax.set_extent(extent)
one_plot(forced)
ax.set_title('Forced')
axs.append(ax)

ax = fig.add_subplot(2, 3, 6, projection=prj)
ax.set_extent(extent)
cf = one_plot(forced)
axs.append(ax)

# add colorbar
plt.subplots_adjust(hspace=0.02, wspace=0.02)

# colorbar and labels
cb = plt.colorbar(cf, shrink=0.5, orientation='horizontal', ax=axs, pad=0.075)
cb.ax.set_title(dsp.trend.units)
plt.savefig('trend-decomp-O2-200m-NPac.png', dpi=300, bbox_inches='tight')
field = [
    dsp.trend.isel(member_id=i).compute() for i in range(0, len(ds.member_id))
]<jupyter_output><empty_output><jupyter_text>Plot All Trends - Total<jupyter_code>fig = plt.figure(figsize=(10, 10))

axs = []

# plot total
for i in range(0, len(ds.member_id)):
    ax = fig.add_subplot(7, 5, i + 1, projection=prj)
    ax.set_extent(extent)
    cf = one_plot(field[i], lines=False)
    axs.append(ax)

plt.subplots_adjust(hspace=0.04, wspace=0.02)

# colorbar and labels
cb = plt.colorbar(cf, shrink=0.5, orientation='horizontal', ax=axs, pad=0.075)
cb.ax.set_title(dsp.trend.units)

plt.savefig('all-trends-O2-200m-NPac.png', dpi=300, bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>Plot All Trends - Internal<jupyter_code>fig = plt.figure(figsize=(10, 10))

axs = []

# plot total
for i in range(0, len(ds.member_id)):
    ax = fig.add_subplot(7, 5, i + 1, projection=prj)
    ax.set_extent(extent)
    cf = one_plot(field[i] - forced, lines=False)
    axs.append(ax)

plt.subplots_adjust(hspace=0.04, wspace=0.02)

# colorbar and labels
cb = plt.colorbar(cf, shrink=0.5, orientation='horizontal', ax=axs, pad=0.075)
cb.ax.set_title(dsp.trend.units)

plt.savefig('all-trends-internal-O2-200m-NPac.png',
            dpi=300,
            bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14895752034786486, 'fraction_numerical': 0.03277957408852715, 'mean_word_length': 3.4471988101140307, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '200681', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3585, 'n_tokens_neox': 3349, 'n_words': 611, 'path': 'notebooks/workflows/cesm/oxygen-trend-computation/O2_forced_unforced.ipynb', 'repo_name': 'zbruick/ncar-python-tutorial', 'size': '8969'}