Regional diagnostics - annual averages¶

Created by Ivan Lima on Mon Oct 11 2021 08:54:21 -0400

In [1]:
import xarray as xr
import pandas as pd
import numpy as np
import holoviews as hv
import os, datetime, my_bokeh_themes
import hvplot.pandas
from cesm_utils import fixtime, spd
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Wed Jul 20 11:41:17 2022
In [2]:
from bokeh.plotting import show
hv.extension('bokeh')

Annual Means¶

In [3]:
from pencil_cases import *
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
# datadir = '/mnt/chromeos/removable/WD/Data/Postproc/1DPOP_CESM2/'
start_yr = 1851 # use "real" year to keep pandas datetime functionality

case_list = ['cesm2.B1850.f09_g17.1dpop.010', 'cesm2.B1850.f09_g17.1dpop.011', 'cesm2.B1850.f09_g17.1dpop.012',
             'cesm2.B1850.f09_g17.1dpop.013', 'cesm2.B1850.f09_g17.1dpop.014', 'cesm2.B1850.f09_g17.1dpop.015']
for case in case_list:
    print('{}: {}'.format(case, case_info[case]))

def plot_tseries(region, varname):
    # read time series data
    data_dict = {}
    cols = []
    for case in case_list:
        col_name = '.'.join(case.split('.')[-2:])
        cols.append(col_name)
        with xr.open_dataset(os.path.join(datadir, '{}.ts_{}.nc'.format(case, region)), 
                             decode_times=False) as ds:
            ds = fixtime(ds, start_yr)
            ds['time'] = ds.time.to_index().to_datetimeindex()
            # convert to annual averages
            data_dict[col_name] = ds[varname].resample(time='A').mean().to_series()
        
    df = pd.DataFrame(data_dict)
    df.columns.name = 'Case'
    ds_ctrl = xr.open_dataset(os.path.join(datadir, 'b.e21.B1850.f09_g17.CMIP6-piControl.001.ts_{}.nc'.format(region)),
                          decode_times=False)
    mean = ds_ctrl[varname].mean(dim='time').values.item()
    std  = ds_ctrl[varname].std(dim='time').values.item()
    df['pi control'] = mean
    df['lo'] = mean - std
    df['hi'] = mean + std
    df['min'] = ds_ctrl[varname].min(dim='time').values.item()
    df['max'] = ds_ctrl[varname].max(dim='time').values.item()
    # create plots
    aspect, height = 1.1, 350
    p_tseries = df[cols + ['pi control']].hvplot(aspect=aspect, frame_height=height, grid=True, shared_axes=False,
                                                 xlabel='time', ylabel='{}'.format(ds[varname].units),title=varname)
    p_std = df[['lo','hi']].hvplot.area(y='lo', y2='hi', alpha=0.2, shared_axes=False, aspect=aspect, frame_height=height,
                                        grid=True, color='green')
    p_minmax = df[['min','max']].hvplot.area(y='min', y2='max', alpha=0.2, shared_axes=False, aspect=aspect, frame_height=height,
                                             grid=True,color='lightgreen')
    return p_tseries * p_std * p_minmax

ds = xr.open_dataset(os.path.join(datadir, '{}.ts_near_global.nc'.format(case_list[0])), decode_times=False)

cesm2.B1850.f09_g17.1dpop.010: 1D + T, S, U & V restoring (tau = 15 days) + EBM off + start year = 1161
cesm2.B1850.f09_g17.1dpop.011: 1D + G-terms (corrected) + SSS restoring (tau = 30 days) + EBM off + start year = 1161
cesm2.B1850.f09_g17.1dpop.012: 1D + G-terms (corrected) + SSS restoring (tau = 365 days) + EBM off + start year = 1161
cesm2.B1850.f09_g17.1dpop.013: 1D + G-terms (corrected) + SSS restoring (tau = 30 days) + EBM off + start year = 6 restoring run
cesm2.B1850.f09_g17.1dpop.014: 1D + G-terms (corrected) + SSS restoring (tau = 365 days) + EBM off + start year = 6 restoring run
cesm2.B1850.f09_g17.1dpop.015: 1D + T, S, U & V restoring (tau = 30 days) + EBM off + start year = 1161

The green shaded areas in the plots below represent the pi control mean plus/minus one standard deviation: $\bar{x} \pm \sigma$. The light green shaded area represents the pi control min and max range.

Near global (60S - 60N)¶

In [4]:
plist_ng = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ng += plot_tseries('near_global',vname)

show(hv.render(plist_ng.cols(2).opts(shared_axes=False)))

North Atlantic¶

In [5]:
plist_na = hv.Layout()
for vname in ds.data_vars.keys():
    plist_na += plot_tseries('north_atl',vname)
    
show(hv.render(plist_na.cols(2).opts(shared_axes=False)))

Tropical Pacific¶

In [6]:
plist_tp = hv.Layout()
for vname in ds.data_vars.keys():
    plist_tp += plot_tseries('trop_pac',vname)
    
show(hv.render(plist_tp.cols(2).opts(shared_axes=False)))

Labrador Sea¶

In [7]:
plist_ls = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ls += plot_tseries('lab_sea',vname)
    
show(hv.render(plist_ls.cols(2).opts(shared_axes=False)))

Arctic¶

In [8]:
plist_ar = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ar += plot_tseries('arctic',vname)
    
show(hv.render(plist_ar.cols(2).opts(shared_axes=False)))

Southern Ocean¶

In [9]:
plist_so = hv.Layout()
for vname in ds.data_vars.keys():
    plist_so += plot_tseries('socean',vname)
    
show(hv.render(plist_so.cols(2).opts(shared_axes=False)))

Mean Annual Cycle¶

In [10]:
def plot_annual_cycle(region, varname):
    # read time series data
    data_dict = {}
    cols = []
    for case in case_list:
        col_name = '.'.join(case.split('.')[-2:])
        cols.append(col_name)
        with xr.open_dataset(os.path.join(datadir, '{}.ts_{}.nc'.format(case, region)), 
                             decode_times=False) as ds:
            ds = fixtime(ds, start_yr)
            ds['time'] = ds.time.to_index().to_datetimeindex()
            # convert to monthly climatology
            data_dict[col_name] = ds[varname].groupby(ds.time.dt.month).mean().to_series()
    
    df = pd.DataFrame(data_dict)
    # create plots
    aspect, height = 1.1, 350
    p_tseries = df.hvplot(aspect=aspect, frame_height=height, grid=True, shared_axes=False, xlabel='month',
                          ylabel='{}'.format(ds[varname].units),title=varname)
    return p_tseries

Near global (60S - 60N)¶

In [11]:
plist_ng = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ng += plot_annual_cycle('near_global',vname)

show(hv.render(plist_ng.cols(2).opts(shared_axes=False)))

North Atlantic¶

In [12]:
plist_na = hv.Layout()
for vname in ds.data_vars.keys():
    plist_na += plot_annual_cycle('north_atl',vname)
    
show(hv.render(plist_na.cols(2).opts(shared_axes=False)))

Tropical Pacific¶

In [13]:
plist_tp = hv.Layout()
for vname in ds.data_vars.keys():
    plist_tp += plot_annual_cycle('trop_pac',vname)
    
show(hv.render(plist_tp.cols(2).opts(shared_axes=False)))

Labrador Sea¶

In [14]:
plist_ls = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ls += plot_annual_cycle('lab_sea',vname)
    
show(hv.render(plist_ls.cols(2).opts(shared_axes=False)))

Arctic¶

In [15]:
plist_ar = hv.Layout()
for vname in ds.data_vars.keys():
    plist_ar += plot_annual_cycle('arctic',vname)
    
show(hv.render(plist_ar.cols(2).opts(shared_axes=False)))

Southern Ocean¶

In [17]:
plist_so = hv.Layout()
for vname in ds.data_vars.keys():
    plist_so += plot_annual_cycle('socean',vname)
    
show(hv.render(plist_so.cols(2).opts(shared_axes=False)))
In [ ]:
# import panel as pn
# row = pn.Row(plist_so.cols(2).opts(shared_axes=False))
# row.embed()
# display(plist_so.cols(2).opts(shared_axes=False))