<jupyter_start><jupyter_text>flag<jupyter_code>ARtag = 'abs'

flag_area = 1000   # minimum size of patches
flag_USstate = 1  # whether to use US west coast 5 states along with land mask. 1 is to use, 0 is to skip
flag_post_adj = 1 # WRF further adjusted, or not (i.e., directly from modified NARR). 1 is further adjusted, 0 for raw

commonAR_thre = 1000

version_tag = 'AR%s_s%d_state%d_post%d_c%d' % (ARtag, flag_area, flag_USstate, flag_post_adj, commonAR_thre)
print(version_tag)<jupyter_output><empty_output><jupyter_text>1. retrieve data<jupyter_code>infile = rootdir + 'data/intermediate_data/AR_stats_separate.%s.mat' % (version_tag)
    
stats_AR_SSTmean = sio.loadmat(infile)['stats_AR_SSTmean']
stats_AR_dist = sio.loadmat(infile)['stats_AR_dist']
stats_AR_landarea = sio.loadmat(infile)['stats_AR_landarea']
stats_AR_IVT = sio.loadmat(infile)['stats_AR_IVT']
stats_AR_IWV = sio.loadmat(infile)['stats_AR_IWV']
stats_AR_IVTs = sio.loadmat(infile)['stats_AR_IVTs']
stats_AR_IWVs = sio.loadmat(infile)['stats_AR_IWVs']
ARday_index = sio.loadmat(infile)['ARday_index']
bg_year = sio.loadmat(infile)['bg_year']
bg_month = sio.loadmat(infile)['bg_month']
commonAR = sio.loadmat(infile)['commonAR'][0]
valid_index = (ARday_index[0])*(ARday_index[1])*(np.abs(stats_AR_SSTmean[0]-stats_AR_SSTmean[1])>=0.1)#*((bg_month[0]>=10)|(bg_month[0]<=3))
valid_index2 = (commonAR==1)*(np.abs(stats_AR_SSTmean[0]-stats_AR_SSTmean[1])>=0.25)

valid_AR_SSTmean = stats_AR_SSTmean[:,valid_index==1]
valid_AR_dist = stats_AR_dist[:, valid_index==1]
valid_AR_landarea = stats_AR_landarea[:, valid_index==1]
valid_AR_IVT = stats_AR_IVT[:, valid_index==1]
valid_AR_IWV = stats_AR_IWV[:, valid_index==1]
valid_AR_IVTs = stats_AR_IVTs[:, valid_index==1]
valid_AR_IWVs = stats_AR_IWVs[:, valid_index==1]

valid_bg_year = bg_year[0, valid_index==1]
valid_bg_month = bg_month[0, valid_index==1]

print('There are %d valid 6-hr snaps with co-occurrence of ARs' % valid_index.sum())
print('alternative: %d' % valid_index2.sum())<jupyter_output><empty_output><jupyter_text>2. difference between ARs<jupyter_code>diff_SST = valid_AR_SSTmean[1]-valid_AR_SSTmean[0]
diff_dist = valid_AR_dist[1]/valid_AR_dist[0]-1
diff_area = valid_AR_landarea[1]/valid_AR_landarea[0]-1
diff_IVT = valid_AR_IVT[1]/valid_AR_IVT[0]-1
diff_IWV = valid_AR_IWV[1]/valid_AR_IWV[0]-1
diff_IVTs = valid_AR_IVTs[1]/valid_AR_IVTs[0]-1
diff_IWVs = valid_AR_IWVs[1]/valid_AR_IWVs[0]-1<jupyter_output><empty_output><jupyter_text>1.1 box plot, grouped by SSTmean difference<jupyter_code>def sub_calc_stats_SSTbin(in_diff_SST, in_stats, xmax=3, xres=0.25):
    for i in np.arange(-1*xmax, xmax+1, xres):
        index = (in_diff_SST>i)*(in_diff_SST<(i+xres))
        bindata = in_stats[index]
        #out_stats = np.mean(bindata)
        out_stats = np.percentile(bindata, 75)
        print('SST diff %.2f-%.2f: %.2f' % (i, i+xres, out_stats))
def compute_regr_coeffs(inx, iny):
    nt = iny.shape[0]
    x = inx
    
    k = ((x*(iny-iny.mean())).sum()) / ((x*(x-x.mean())).sum())
    b = iny.mean() - k * x.mean()
    
    return k, b<jupyter_output><empty_output><jupyter_text>Note:  the following two barplot functions take some parameters from background, so make sure they are set correctly.<jupyter_code>import scipy.stats.distributions as dist


def check_mean_zero_significance(indata, thre_p=0.1):
    n = indata.shape[0]
    se = indata.std()/np.sqrt(n)
    be = indata.mean()
    he = 0
    test_statistic = (be - he)/se
    
    pvalue = 2*dist.norm.cdf(-np.abs(test_statistic))
    
    if np.abs(pvalue)<thre_p:
        sig = 1
    else:
        sig = 0
    
    return sig
def barplot_single_set_with_sig(axis, indata, inlabel, incolors):
    patch_collection_for_legend = []
    for i in np.arange(nbins):
        threb = -1*xmax+i*xres
        threu = -1*xmax+(i+1)*xres
        loc_v = (threb+threu)/2
        index = (diff_SST>threb)*(diff_SST<=threu)
        if index.sum()>20:
            plotdata = indata[index]
            sig = check_mean_zero_significance(plotdata)
            bplot1 = axis.boxplot(plotdata, positions=[loc_v], widths=para_a,
                        showfliers=False, whis=[10,90], patch_artist=True)
            for patch, color in zip(bplot1['boxes'], incolors):
                if sig==1:
                    patch.set_facecolor(color)
                    patch.set_edgecolor(color)
                else:
                    patch.set_facecolor('none')
                    patch.set_edgecolor('black')
                if i==5:  # 7 for abs
                    patch_collection_for_legend.append(patch)
            if i==5: # 7 for abs
                leg = axis.legend(patch_collection_for_legend,
                                 [inlabel], ncol=1, loc='upper right',
                                 frameon=False, fontsize=10)
                for legobj in leg.legendHandles:
                    legobj.set_linewidth(0)
def barplot_double_set_sig(axis, indata, inlabels, incolors):
    patch_collection_for_legend = []
    for i in np.arange(nbins):
        threb = -1*xmax+i*xres
        threu = -1*xmax+(i+1)*xres
        loc_v = (threb+threu)/2
        index = (diff_SST>threb)*(diff_SST<=threu)
        if index.sum()>20:
            plotdata = [indata[0][index], indata[1][index]]
            sigs = [check_mean_zero_significance(plotdata[0]), check_mean_zero_significance(plotdata[1])]
            #print(sigs)
            #bplot1 = ax1.boxplot(plotdata, positions=np.arange(loc_v-1*para_a-1*para_b, loc_v+1*para_a+1*para_b+0.0000001, para_a+para_b), widths=para_a,
            #            showfliers=False, whis=[10,90], patch_artist=True)
            bplot1 = axis.boxplot(plotdata, positions=np.arange(loc_v-0.5*para_a-0.5*para_b, loc_v+0.5*para_a+0.5*para_b+0.0000001, para_a+para_b), widths=para_a,
                        showfliers=False, whis=[10,90], patch_artist=True)
            for patch, color, sig in zip(bplot1['boxes'], incolors, sigs):
                if sig==1:
                    patch.set_facecolor(color)
                    patch.set_edgecolor(color) 
                if sig==0:
                    patch.set_facecolor('none')
                    patch.set_edgecolor('black')
                if i==11: # 7 for abs
                    patch_collection_for_legend.append(patch)
            if i==11: # 7 for abs
                leg = axis.legend(patch_collection_for_legend,
                                 inlabels, ncol=1, loc='upper right',
                                 frameon=False, fontsize=10)
                for legobj in leg.legendHandles:
                    legobj.set_linewidth(0)
para_a = 0.11 # bin width
para_b = 0 # bin interval
colors4 = ['royalblue', 'lightseagreen']
colors_area = ['orchid']

fig2 = plt.figure(figsize=(7,6))
rowspan=9
colspan=10
ax1 = plt.subplot2grid((21,23), (0,0), colspan=colspan, rowspan=rowspan)
ax2 = plt.subplot2grid((21,23), (0,11), colspan=colspan, rowspan=rowspan)
ax3 = plt.subplot2grid((21,23), (10,0), colspan=colspan, rowspan=rowspan)
ax4 = plt.subplot2grid((21,23), (10,11), colspan=colspan, rowspan=rowspan)
axes_full = [ax1, ax2 ,ax3, ax4]

xmax = 3
xres = 0.5
nbins = int(xmax*2/xres)
H_interval=xres

barplot_single_set_with_sig(ax1, diff_area, 'landfalling area', ['orchid'])
barplot_single_set_with_sig(ax2, diff_dist, 'distance', ['orchid'])
barplot_double_set_sig(ax3, [diff_IVT, diff_IWV], ['mean IVT', 'mean IWV'], ['royalblue', 'lightseagreen'])
barplot_double_set_sig(ax4, [diff_IVTs, diff_IWVs], ['int.IVT', 'int.IWV'], ['royalblue', 'lightseagreen'])

# decoration
for i in np.arange(nbins):
    threb = -1*xmax+i*xres
    threu = -1*xmax+(i+1)*xres
    for axis in axes_full:
        # background       
        if i%2==0:
            patch = matplotlib.patches.Rectangle((threb,-100), H_interval, 200, color='white', alpha=1, zorder=0)
        else:
            patch = matplotlib.patches.Rectangle((threb,-100), H_interval, 200, color='lightgray', alpha=0.6, zorder=0)
        axis.add_patch(patch)
         
for axis in axes_full:
    axis.plot(np.arange(-20, 20), np.zeros(40), linestyle='--', linewidth=0.5, color='black')
    axis.set_xlim([-2, xmax])
    axis.set_xticks(np.arange(-2, xmax+0.000001, 1))
    axis.set_xticklabels(np.arange(-2, xmax+xres+0.000001, 1, dtype=int))
    axis.set_ylabel('Relative change', fontsize=12)

ax4.set_ylim([-0.2, 0.4])
ax4.set_yticks(np.arange(-0.2, 0.41, 0.2))
ax4.set_yticklabels(['-20%', '0', '+20%', '+40%'])
    
for axis in [ax1, ax2]:
    axis.set_ylim([-0.15, 0.25])
    axis.set_yticks(np.arange(-0.1, 0.21, 0.1))
    axis.set_yticklabels(['-10%', '0', '+10%', '+20%'])
    
ax3.set_ylim([-0.03, 0.06])
ax3.set_yticks(np.arange(-0.03, 0.06+0.000001, 0.03))
ax3.set_yticklabels(['-3%', '0', '+3%', '+6%'])


for axis in [ax2, ax4]:
    axis.yaxis.set_label_position("right")
    axis.yaxis.tick_right()
    
for axis in [ax1, ax2]:
    axis.set_xlabel('')
    axis.set_xticklabels([])

for axis in [ax3, ax4]:
    axis.set_xlabel('SST perturbation (K)', fontsize=12)

ax1.text(-1.9, 0.25-0.35*0.05, '(a)', ha='left', va='top', fontsize=12)
ax2.text(-1.9, 0.25-0.35*0.05, '(b)', ha='left', va='top', fontsize=12)
ax3.text(-1.9, 0.06-0.09*0.05, '(c)', ha='left', va='top', fontsize=12)
ax4.text(-1.9, 0.4-0.6*0.05, '(d)', ha='left', va='top', fontsize=12)

if ARtag=='abs':
    figname = plotdir + 'fig04.AR_sens.with_sig_test.%s.png'%version_tag
elif ARtag=='p85':
    figname = plotdir + 'fig_S7.AR_sens.with_sig_test.%s.png'%version_tag

#fig2.savefig(figname, dpi=600)

plt.show()
del(fig2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1566189530874517, 'fraction_numerical': 0.04461393793751959, 'mean_word_length': 3.456238361266294, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '579257', 'license': 'MIT', 'n_tokens_mistral': 4061, 'n_tokens_neox': 3811, 'n_words': 614, 'path': 'scripts/fig4.figS7.AR_sen.with_stats_significance.ipynb', 'repo_name': 'lucas-uw/Chen-2020-GRL', 'size': '9571'}