<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Data loading1.1&nbsp;&nbsp;Tropopause1.2&nbsp;&nbsp;LWA1.3&nbsp;&nbsp;Stat. significance2&nbsp;&nbsp;Figure 3<jupyter_code>import glob
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import shapely.geometry as sgeom
from matplotlib.ticker import FuncFormatter
from palettable.scientific.diverging import Lisbon_20, Tofino_20
%reload_ext watermark
%watermark --iversions -v -m
timescale = 20
DJF_bool = '_DJFonly'
root_path = '/mnt/nas4.meop2/meop40.data.model/CMAM/0A.daily/'
cesta = f'{root_path}composites_woSSW{DJF_bool}/'<jupyter_output><empty_output><jupyter_text>Data loading  Tropopause<jupyter_code>infile_name = '{}{}_{}_comp_*_{}days.nc'.format(cesta, 'ptrop', 'absolute', 20)
infiles = glob.glob(infile_name)
e_types = list(map(lambda x: x.split('comp_')[1].split('_{}days'.format(timescale))[0], infiles))
da_trop =  xr.open_mfdataset(infiles, concat_dim='e_type', combine = 'nested')
da_trop['e_type'] = e_types
da_trop['ptrop'] = da_trop.ptrop/100.
da_trop<jupyter_output><empty_output><jupyter_text>LWA<jupyter_code>var_ls = ['lwatend_logH']
e_type_ls = ['himalayas', 'eastasia', 'westamer']
what_ls = ['anomalies','absolute']
da_ls2 = []

for what in what_ls:
    da_ls = []
    print(what)
    for i,e_type in enumerate(e_type_ls):
        print(e_type)
        ds_ls = []
        for var in var_ls:
            ds = xr.open_dataset(f'{cesta}{var}_{what}_comp_{e_type}_{timescale}days.nc')
            if var not in ['lwa','sink', 'lwatend_logH']:
                ds['plev'] = ds.plev/100.
            else:
                ds = ds.rename({'z': 'plev'})
            ds_ls.append(ds)
        da = xr.merge(ds_ls)

        da_ls.append(da)
    da_conc = xr.concat(da_ls, dim = 'e_type')
    da_conc['e_type'] = e_type_ls
    da_ls2.append(da_conc)
        
da_all = xr.concat(da_ls2, dim = 'what')
da_all['what'] = what_ls
da_all<jupyter_output>anomalies
himalayas
eastasia
westamer
absolute
himalayas
eastasia
westamer<jupyter_text>Stat. significance<jupyter_code>#infile_name = '{3}{0}_pvalues_from{1}_comp_*_{2}days_zm.nc'.format(var, 10000, timescale, cesta)
var = 'lwatend'#'sink'
infile_name = '{3}{0}_pvalues_from{1}_comp_*_{2}days.nc'.format(var, 10000, timescale, cesta)
print(infile_name)
infiles = glob.glob(infile_name)
comp_types = list(map(lambda x: x.split('comp_')[1].split('_{}days'.format(timescale))[0], infiles))
#print(comp_types)
#print(infiles)
#print()
ds_p = xr.open_mfdataset(infiles, concat_dim='e_type', combine ='nested')
ds_p['e_type'] = comp_types
if var not in ['lwa','sink','lwatend']:
    ds_p['plev'] = ds_p.plev/100.
else:
    ds_p = ds_p.rename({'z': 'plev'})
ds_p<jupyter_output>/mnt/nas4.meop2/meop40.data.model/CMAM/0A.daily/composites_woSSW_DJFonly/lwatend_pvalues_from10000_comp_*_20days.nc<jupyter_text>Figure 3<jupyter_code>def thousands(x, pos):
    return '%g' % (x*1e-3)

p0=1013.25
H = 7000
scale = 24*3600
def press2height(x):
    return H*np.log(p0/x)
      
def height2press(x):
    return p0*np.exp(-x/H)
np.linspace(-3,3,41)
#levels= np.linspace(-2.5,2.5,41)
levels = np.arange(-7.2,7.6,0.4)
#levels= np.linspace(-2.5,2.5,41)
levels2=[0.5,1,3]
lev0 = [0]
cbar_label = f'FAWA$_t$ [m/s/day]'
titles = ['Himalayas', 'East Asia', 'Rocky Mountains'] #'West America']
label_ls = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I']
cmap = mpl.cm.PuOr_r#RdBu_r
#cmap = Tofino_20.mpl_colormap
mpl.rcParams.update({'font.size': 18})

plot_arr = da_all['lwatend'].sel(what = what_ls[0], lag = 0, plev = slice(None,60000)).mean('lon')*scale
p = xr.plot.FacetGrid(plot_arr, col = 'e_type', size = 5)

da_trop2 = da_trop.ptrop.sel(lag=0, lat = slice(0,None)).mean('lon')
lats = p.data.lat.values
levs = p.data.plev.values
plot_arr = da_all['lwatend'].sel(what = what_ls[1], lag = 0, plev = slice(None,60000)).mean('lon')*scale

p.axes[0][0].set_xlabel('')
p.axes[0][1].set_xlabel('latitude [degree]')
p.axes[0][2].set_xlabel('')


for i,(ax, sel_dict, label) in enumerate(zip(p.axes.flat, p.name_dicts.flat, label_ls)):
    ax.text(0.05, 0.95, label, transform=ax.transAxes, fontsize = 20, 
              fontweight='bold', va='top', bbox={'alpha': 0.7, 'facecolor': 'white'})
    #da_trop2.sel(**sel_dict).plot.line('k--', ax = ax, lw = 3, zorder = 20)
    
    if i==0:
        ax.set_ylabel('pressure [hPa]', fontsize = 19)
    else:
        ax.set_ylabel('')
    
    ax.set_yscale('log')
    #ax.get_shared_y_axes()
    ax.axis([0, 90, height2press(0), height2press(60000)])
    ax.set_yticks([1000,100,10,1])   
    fmt = mpl.ticker.FormatStrFormatter("%g")
    ax.yaxis.set_major_formatter(fmt)       
    
    ax2 = ax.twinx()
    # plot contourf of lwatend anomalies
    plot_kwargs = dict(ax = ax2, \
                   cmap = cmap, \
                   levels=levels, add_colorbar = False)
    p.data.sel(**sel_dict).plot.contourf(**plot_kwargs)
    # plot contour of lwatend absolute values (positive only)
    plot_kwargs = dict(ax = ax2,\
                   colors = ['grey'], \
                   levels=levels2, add_colorbar=False, alpha =1)
    plot_arr.sel(**sel_dict).plot.contour(**plot_kwargs)
    # plot contour of lwatend absolute values (negative only)
    plot_kwargs['linestyles'] = 'dashed'
    plot_kwargs['levels'] = np.array(levels2)[::-1]*-1
    plot_arr.sel(**sel_dict).plot.contour(**plot_kwargs)
    # plot contour of lwatend absolute values (zero only)
    plot_kwargs['linestyles'] = 'solid'
    plot_kwargs['levels'] = lev0
    plot_kwargs['linewidths'] = 3
    plot_kwargs['alpha'] = 1
    plot_arr.sel(**sel_dict).plot.contour(**plot_kwargs)  
    # plot tropopause line
    press2height(da_trop2.sel(**sel_dict)).plot.line('k--', ax=ax2, lw = 3)
    # plot hotspot lines
    ax2.hlines(press2height(70),-20,90, color='g', ls='--', lw=3)
    
    e_type = sel_dict['e_type']
    if e_type == 'himalayas':
        box = sgeom.box(minx=70, maxx=105, miny=20, maxy=42)
    elif e_type == 'eastasia':
        box = sgeom.box(minx=110, maxx=145, miny=27.5, maxy=52)
    elif e_type == 'westamer':
        box = sgeom.box(minx=-125, maxx=-102.5, miny=27.5, maxy=52)
    
    ax2.vlines(list(box.bounds)[1::2],press2height(90),press2height(50), color='g', ls='--', lw=3)
    # plot stat. significance using hatching
    plot_kwargs = dict(levels = [0,0.05], hatches=['\\\\',None], \
                                colors='none', add_colorbar=False)
    sel_dict['plev'] = slice(None, 60000)
    sel_dict['lat'] = slice(-20,None)
    sel_dict['lag'] = 0
    ds_p['lwatend'].sel(**sel_dict).plot.contourf(**plot_kwargs)
    
    ax2.get_shared_y_axes()
    ax2.axis([0, 90, 0, 60000])
    ax2.set_yticks([10000,30000,50000])
    fmt = FuncFormatter(thousands)
    ax2.yaxis.set_major_formatter(fmt)
    
    if i != 2:
        ax2.set_yticklabels([])
        ax2.set_ylabel('')
    if i == 2:
        ax2.set_ylabel('height [km]', fontsize = 19)
    
    ax2.set_title(f'{titles[i]}', fontdict=dict(fontweight = 'bold', fontsize = 22))  
    
    
p.fig.tight_layout()
p.fig.subplots_adjust(right=0.875)
cbar_ax = p.fig.add_axes([0.075, -0.03, .8, .03])
norm = mpl.colors.BoundaryNorm(levels[1:-1], cmap.N)
label = cbar_label
cb2 = mpl.colorbar.ColorbarBase(cbar_ax, cmap=cmap,  \
                                norm=norm, orientation='horizontal', \
                                boundaries=levels, extend='both',  label = cbar_label, \
                                ticks=[-6,-4,-2,0,2,4,6], spacing='proportional')

outfile='{}_{}_{}_{}days_zm_wsignificance_DJFonly.pdf'.format('lwatend', what_ls[0], 'all', timescale)
print(outfile)
plt.savefig(outfile, dpi = 300, bbox_inches='tight')<jupyter_output>lwatend_anomalies_all_20days_zm_wsignificance_DJFonly.pdf
{'fraction_non_alphanumeric': 0.17451309170643622, 'fraction_numerical': 0.04295111569714949, 'mean_word_length': 3.8553537883531623, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '586190', 'license': 'MIT', 'n_tokens_mistral': 3410, 'n_tokens_neox': 3152, 'n_words': 496, 'path': 'code/GRL_reproduce_Fig3.ipynb', 'repo_name': 'kuchaale/grl_2020', 'size': '7753'}