Diagnostic plots - Equatorial Pacific (0-3N, 74W-248W)

Created by Ivan Lima on 2020-05-04

%matplotlib inline
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os, datetime
from cesm_utils import fixtime
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Wed Jul  7 10:28:33 2021

import seaborn as sns
sns.set_theme(context='paper', palette='tab10')
plt.rcParams['figure.dpi'] = 100

Read model output

case = 'cesm2.B1850.f09_g17.1dpop.005'
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
start_yr = 1851 # use "real" year to keep pandas datetime functionality 
ds = xr.open_dataset(os.path.join(datadir, '{}.eq_pac.daily.nc'.format(case)), decode_times=False)
ds = fixtime(ds, start_yr)
ds['time'] = ds.time.to_index().to_datetimeindex() # convert to datetime index for pandas
ds['TAREA'] = ds.TAREA.where(ds.REGION_MASK>0)

# ds_clim = xr.open_dataset(os.path.join(datadir, '../1DPOP/forcing/clim_eqpac.nc'), decode_times=False)
# ds_clim = fixtime(ds_clim,2051)
# ds_clim['time'] = ds_clim.time.to_index().to_datetimeindex() # convert to datetime index for pandas
# for vname in ['dz','TAREA']:
#     ds_clim[vname] = ds[vname]

Time series of mean values

def get_mean(vname, ds=ds):
    return ((ds[vname] * ds.dz * ds.TAREA).sum(dim=['z_t','nlat','nlon']) /
            (ds.dz * ds.TAREA).sum()).to_series()

varlist = ['TEMP','SALT','UVEL','VVEL']
data = {vname: get_mean(vname) for vname in varlist}
df_mean = pd.DataFrame(data)
df_mean = df_mean.stack().reset_index().rename(columns={'level_1':'variable', 0:'value'})

# for vname in varlist:
#     df_mean[vname+'_clim'] = get_mean(vname, ds_clim)

g = sns.relplot(x='time', y='value', data=df_mean, kind='line', col='variable', col_wrap=2,
                height=4, facet_kws={'sharey': False, 'sharex': True})
_ = g.set_xticklabels(rotation=-45, ha='left')
_ = g.set_titles(col_template='mean {col_name}')
for ax, vname in zip(g.axes.flat, varlist):
    ax.set_ylabel('{}'.format(ds[vname].units))
    ax.autoscale(axis='x', tight=True)

Time series of absolute max values

def get_max(vname):
    return np.abs(ds[vname]).max(dim=['z_t','nlat','nlon']).to_series()

data = {vname: get_max(vname) for vname in varlist}
df_max = pd.DataFrame(data)
df_max = df_max.stack().reset_index().rename(columns={'level_1':'variable', 0:'value'})

g = sns.relplot(x='time', y='value', data=df_max, kind='line', col='variable', col_wrap=2,
                height=4, facet_kws={'sharey': False, 'sharex': True})
_ = g.set_xticklabels(rotation=-45, ha='left')
_ = g.set_titles(col_template='absolute max {col_name}')
for ax, vname in zip(g.axes.flat, varlist):
    ax.set_ylabel('{}'.format(ds[vname].units))
    ax.autoscale(axis='x', tight=True)