Diagnostic plots 2 - daily output

Created by Ivan Lima on Wed Apr 14 2021 16:09:59 -0400

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

Last updated on Wed Jul  7 10:45:26 2021

warnings.filterwarnings('ignore')
plt.rcParams['figure.max_open_warning'] = 50
plt.rcParams['figure.dpi'] = 100

Time series of absolute max values

case = 'cesm2.B1850.f09_g17.1dpop.006'
datadir = '/home/ivan/Data/Postproc/1DPOP_CESM2'
start_yr = 1851 # use "real" year to keep pandas datetime functionality 

ds_max = xr.open_dataset(os.path.join(datadir, '{}.abs.max.nc'.format(case)), decode_times=False)
ds_max = fixtime(ds_max, start_yr)
ds_max['time'] = ds_max.time.to_index().to_datetimeindex()

data_dict = {vname: ds_max[vname].to_series() for vname in ['HMXL','TEMP','SALT','UVEL','VVEL']}
df = pd.DataFrame(data_dict)
df = df.stack().reset_index().rename(columns={'level_1':'variable', 0:'value'})

with sns.axes_style('darkgrid'):
    g = sns.relplot(x='time', y='value', data=df, kind='line', col='variable', col_wrap=3,
                    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, ['HMXL','TEMP','SALT','UVEL','VVEL']):
        ax.set_ylabel('{}'.format(ds_max[vname].units))
        ax.autoscale(axis='x', tight=True)

Maps of absolute max values

ds_2D_max = xr.open_dataset(os.path.join(datadir, '{}.abs.max_2D.nc'.format(case)), decode_times=False, chunks='auto')
ds_2D_max = fixtime(ds_2D_max, start_yr)

The maps below show the absolute max value at each grid point across depth and time.

from cartopy import crs as ccrs
from cartopy import feature as cfeature

tlat, tlon = ds_2D_max.TLAT, ds_2D_max.TLONG
ulat, ulon = ds_2D_max.ULAT, ds_2D_max.ULONG
tlat, tlon = np.c_[tlat,tlat[:,0]], np.c_[tlon,tlon[:,0]]
ulat, ulon = np.c_[ulat,ulat[:,0]], np.c_[ulon,ulon[:,0]]

def surf_map(vname, cmap=plt.cm.viridis, log=False):
    fig, ax = plt.subplots(subplot_kw={'projection':ccrs.Robinson(central_longitude=210)}, figsize=(9,6.4))
    ax.set_global()
#     ax.add_feature(cfeature.LAND, facecolor='#b0b0b0')
    ax.coastlines(linewidth=0.5)
    if log:
        norm = norm=colors.LogNorm()
    else:
        norm = None    
    data = ds_2D_max[vname].max(dim=['time'])
    pm = ax.pcolormesh(tlon, tlat, data, transform=ccrs.PlateCarree(), norm=norm, cmap=cmap)
    _ = ax.set_title('{} absolute max'.format(vname))
    cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.05)
    cb.set_label(ds_2D_max[vname].units)

surf_map('TEMP', cmocean.cm.thermal)
surf_map('SALT', cmocean.cm.haline)
surf_map('HMXL', cmocean.cm.deep, log=True)
surf_map('IFRAC', plt.cm.BuPu)
for vname in ['UVEL','VVEL']:
    surf_map(vname, plt.cm.cubehelix_r)

lonmin, lonmax = 98, 167
latmin, latmax =  -15, 15
lonc = (lonmin + lonmax)/2.
latc = (latmin + latmax)/2.
# proj = ccrs.LambertConformal(central_longitude=lonc,central_latitude=latc)
proj = ccrs.Mercator(central_longitude=lonc)

for vname in ['UVEL','VVEL']:
    fig, ax = plt.subplots(subplot_kw={'projection':proj}, figsize=(9,6.4))
    ax.set_extent([lonmin,lonmax,latmin,latmax])
    ax.coastlines(linewidth=0.5)
    gl = ax.gridlines(xlocs=np.arange(90,180,10),ylocs=np.arange(-20,30,10),draw_labels=True)
    gl.xlabels_top=False
    gl.ylabels_right=False
    data = ds_2D_max[vname].max(dim=['time'])
    pm = ax.pcolormesh(tlon, tlat, data, transform=ccrs.PlateCarree(), cmap=plt.cm.cubehelix_r)
    _ = ax.set_title('{} absolute max'.format(vname))
    cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.075)
    cb.set_label(ds_2D_max[vname].units)

lonmin, lonmax = 360-170, 360-78
latmin, latmax =  -14, 14
lonc = (lonmin + lonmax)/2.
latc = (latmin + latmax)/2.
proj = ccrs.Mercator(central_longitude=lonc)

for vname in ['UVEL','VVEL']:
    fig, ax = plt.subplots(subplot_kw={'projection':proj}, figsize=(9.5,6))
    ax.set_extent([lonmin,lonmax,latmin,latmax])
    ax.coastlines(linewidth=0.5)
    gl = ax.gridlines(xlocs=np.arange(-180,190,10),ylocs=np.arange(-20,30,10),draw_labels=True)
    gl.xlabels_top=False
    gl.ylabels_right=False
    data = ds_2D_max[vname].max(dim=['time'])
    pm = ax.pcolormesh(tlon, tlat, data, transform=ccrs.PlateCarree(), cmap=plt.cm.cubehelix_r)
    _ = ax.set_title('{} absolute max'.format(vname))
    cb = fig.colorbar(pm, ax=ax, orientation='horizontal', pad=0.075)
    cb.set_label(ds_2D_max[vname].units)