<jupyter_start><jupyter_text>flags<jupyter_code>ARtag = 'p85'
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. data 1.1 background data, SST<jupyter_code>def sub_compute_SSTmean(in_SST, ocean_mask):
    mask_full = ocean_mask==1
    out_SSTmean_all = in_SST[mask_full].mean()
    
    mask_north = (ocean_mask==1)*(wrf_lats>=40)
    out_SSTmean_north = in_SST[mask_north].mean()
    
    mask_south = (ocean_mask==1)*(wrf_lats<40)
    out_SSTmean_south = in_SST[mask_south].mean()
    
    return np.array([out_SSTmean_all, out_SSTmean_north, out_SSTmean_south])
def compute_SSTmean_series():
    reffile = rootdir + 'data/common_ref/SERDP6km.dist_to_coastal.nc'
    dist_to_coast = xr.open_dataset(reffile).dist_to_coast.values
    dist_to_coast[dist_to_coast==9999] = 0
    ocean_mask = np.zeros((450,450))
    ocean_mask[dist_to_coast==0] = 1
    
    infile = rootdir + 'data/HIST/SST_background/NARR_TS.SERDP6km.6hourly.monsum.200310-201512.nc'
    in_SSTmonmean_full = xr.open_dataset(infile).var11.values[0:144]
    ts_monthly = pd.period_range(start='2003-10', end='2015-09', freq='M')

    infile = rootdir + 'data/HIST/SST_background/NARR_TS.SERDP6km.2000.10.01.00.nc'
    fSST_full = xr.open_dataset(infile).var11.values[0]
    
    SSTmean_mseries = np.zeros((144,3))
    for i in np.arange(144):
        SSTmean_mseries[i] = sub_compute_SSTmean(in_SSTmonmean_full[i], ocean_mask)
        
    fSSTmean_mseries = sub_compute_SSTmean(fSST_full, ocean_mask)
        
    return SSTmean_mseries, fSSTmean_mseries
SSTmean_mseries, fSSTmean_mseries = compute_SSTmean_series()
ts_monthly = pd.period_range(start='2003-10', end='2015-09', freq='M')<jupyter_output><empty_output><jupyter_text>1.2 AR data<jupyter_code>tmpfile = rootdir + 'data/intermediate_data/ARstats.monthly_count.%s.mat' % version_tag
count_HIST1 = sio.loadmat(tmpfile)['count_HIST1']
count_HIST2 = sio.loadmat(tmpfile)['count_HIST2']
count_fSST1 = sio.loadmat(tmpfile)['count_fSST1']
count_fSST2 = sio.loadmat(tmpfile)['count_fSST2']

count_HIST = count_HIST1 + count_HIST2
count_fSST = count_fSST1 + count_fSST2

count_subE = count_fSST1 + count_HIST2
count_subF = count_HIST1 + count_fSST2<jupyter_output><empty_output><jupyter_text>2. plotsgiven the same SSTmean change, how the histogram of frequencies change?<jupyter_code>def compute_freq_hist_given_SSTchg(in_dSST, region='all', month='all', dt=2):
   
    if region=='all':
        ii = 0  # used in SSTmean array to identify regions
        sub_index = USstate==0
    elif region=='north':
        ii = 1
        sub_index = (wrf_lats>=40)*(USstate==0)
    elif region=='south':
        ii = 2
        sub_index = (wrf_lats<40)*(USstate==0)
        
    ngrids = sub_index.sum()
        
        
    dSST_full = fSSTmean_mseries[ii]-SSTmean_mseries[:,ii]
    
    if month=='all':
        valid_mindex = (dSST_full>in_dSST[0]) * (dSST_full<=in_dSST[1]) # SST difference range
    elif month=='winter':
        valid_mindex = (dSST_full>in_dSST[0]) * (dSST_full<=in_dSST[1])*((ts_monthly.month>=10)|(ts_monthly.month<=3))
    elif month=='summer':
        valid_mindex = (dSST_full>in_dSST[0]) * (dSST_full<=in_dSST[1])*((ts_monthly.month>=4)&(ts_monthly.month<=9))
    else:
        valid_mindex = (dSST_full>in_dSST[0]) * (dSST_full<=in_dSST[1])*(ts_monthly.month==month)
    
    
    tmp_subset_HIST = []
    tmp_subset_fSST = []
    tmp_subset_subE = []
    tmp_subset_subF = []
    for t in np.arange(144):
        if valid_mindex[t]==1:
            tmp_subset_HIST = np.concatenate((tmp_subset_HIST, count_HIST[t][sub_index]))
            tmp_subset_fSST = np.concatenate((tmp_subset_fSST, count_fSST[t][sub_index]))
            tmp_subset_subE = np.concatenate((tmp_subset_subE, count_subE[t][sub_index]))
            tmp_subset_subF = np.concatenate((tmp_subset_subF, count_subF[t][sub_index]))
            
    
    hist_HIST, bin_edges = np.histogram(tmp_subset_HIST, bins=np.arange(0, 60.1, dt))
    hist_fSST, _ = np.histogram(tmp_subset_fSST, bins=np.arange(0, 60.1, dt))
    hist_subE, _ = np.histogram(tmp_subset_subE, bins=np.arange(0, 60.1, dt))
    hist_subF, _ = np.histogram(tmp_subset_subF, bins=np.arange(0, 60.1, dt))
    
    hist_HIST[hist_HIST<=50] = 0
    hist_fSST[hist_fSST<=50] = 0
    hist_subE[hist_subE<=50] = 0
    hist_subF[hist_subF<=50] = 0
    
    bin_centers = (bin_edges[0:(bin_edges.shape[0]-1)]+bin_edges[1:bin_edges.shape[0]])/2

    return hist_HIST, hist_fSST, hist_subE, hist_subF, bin_edges, bin_centers
dt = 4
nlen = int(60/dt)

histogram_HIST_stack = np.zeros((8, nlen))
histogram_fSST_stack = np.zeros((8, nlen))
histogram_subE_stack = np.zeros((8, nlen))
histogram_subF_stack = np.zeros((8, nlen))

for i in np.arange(8):
    d_range = [i-2, i-1]
    tmp_HIST, tmp_fSST, tmp_subE, tmp_subF, bin_edges, bin_centers = compute_freq_hist_given_SSTchg(in_dSST=d_range, region='all', 
                                                                                month='winter', dt=dt)
    #print(histogram_fSST_stack[i].shape, tmp_fSST.shape)
    histogram_HIST_stack[i] = tmp_HIST
    histogram_fSST_stack[i] = tmp_fSST
    histogram_subE_stack[i] = tmp_subE
    histogram_subF_stack[i] = tmp_subF
se = 0
se_pdf = 0
nlen = histogram_HIST_stack.shape[1]

pdf_HIST = np.zeros((8, nlen-se))
pdf_fSST = np.zeros((8, nlen-se))
pdf_subE = np.zeros((8, nlen-se))
pdf_subF = np.zeros((8, nlen-se))

for i in np.arange(8):
    pdf_HIST[i] = histogram_HIST_stack[i, se:nlen]/histogram_HIST_stack[i, se_pdf:nlen].sum()
    pdf_fSST[i] = histogram_fSST_stack[i, se:nlen]/histogram_HIST_stack[i, se_pdf:nlen].sum()
    pdf_subE[i] = histogram_subE_stack[i, se:nlen]/histogram_HIST_stack[i, se_pdf:nlen].sum()
    pdf_subF[i] = histogram_subF_stack[i, se:nlen]/histogram_HIST_stack[i, se_pdf:nlen].sum()
def derive_color(in_value, cmapname='temp_19lev'):
    cmap_linecolor = generate_NCL_cmap(cmapname)
    if in_value<0:
        return cmap_linecolor((in_value*3+6)/12)
    else:
        return cmap_linecolor((in_value+6)/12)  # so for -6 and 6<jupyter_output><empty_output><jupyter_text>2.1 by each SST warming K<jupyter_code>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), (11,0), colspan=colspan, rowspan=rowspan)
ax4 = plt.subplot2grid((21,23), (11,11), colspan=colspan, rowspan=rowspan)
axes_full = [ax1, ax2 ,ax3, ax4]

for i in np.arange(8):
    
    ax1.bar(bin_edges[se:nlen]+i*0.4+0.6, pdf_HIST[i], width=0.4, 
            color=[derive_color(i-1.5, cmapname='BlueRed')], label='%dK-%dK'%(i-2, i-1))
    
    ax2.bar(bin_edges[se:nlen]+i*0.4+0.6, pdf_fSST[i]-pdf_HIST[i], width=0.4, 
            color=[derive_color(i-1.5, cmapname='BlueRed')], label='%dK-%dK'%(i-2, i-1))
    
    ax3.bar(bin_edges[se:nlen]+i*0.4+0.6, pdf_subE[i]-pdf_HIST[i], width=0.4, 
            color=[derive_color(i-1.5, cmapname='BlueRed')], label='%dK-%dK'%(i-2, i-1))
    
    ax4.bar(bin_edges[se:nlen]+i*0.4+0.6, pdf_subF[i]-pdf_HIST[i], width=0.4, 
            color=[derive_color(i-1.5, cmapname='BlueRed')], label='%dK-%dK'%(i-2, i-1))

    # background
    for axis in axes_full:
        if i%2==0:
            patch = matplotlib.patches.Rectangle((bin_edges[i],-1), 4, 2, color='white', alpha=1, zorder=0)
        else:
            patch = matplotlib.patches.Rectangle((bin_edges[i],-1), 4, 2, color='lightgray', alpha=0.6, zorder=0)
        axis.add_patch(patch)
    
ax1.legend(ncol=2, frameon=False, fontsize=9)

for axis in axes_full:
    axis.plot(np.arange(40), np.zeros(40), linestyle='--', linewidth=0.5, color='black')
    axis.set_xlim([0, 36])
    axis.set_xticks(np.arange(0, 36.1, 8))
    axis.set_xticklabels(np.arange(0, 36.1, 8, dtype=int)*6)
    
ax1.set_ylim([0, 1])
for axis in [ax2, ax3, ax4]:
    if ARtag=='p85':
        axis.set_ylim([-0.075, 0.05])
        axis.set_yticks([-0.075, -0.05, -0.025, 0, 0.025, 0.05])
    elif ARtag=='abs':
        axis.set_ylim([-0.04, 0.02])
        axis.set_yticks(np.arange(-0.04, 0.02+0.00001, 0.02))
    axis.set_ylabel('Change of fraction', fontsize=12)
    
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('Hours of AR per month', fontsize=12)

ax1.set_title('(a) Dist. of AR frequency: CTRL', fontsize=11)
ax1.set_ylabel('Fraction', fontsize=12)
ax2.set_title('(b) Change in dist. (fSST-CTRL)', fontsize=11)
ax3.set_title('(c) Contributions from AR ext.', fontsize=11)
ax4.set_title('(d) Contributions from AR freq.', fontsize=11)

if ARtag=='abs':
    figname = plotdir + 'fig_S2.land_perspect.AR_freq_change.%s.png'%version_tag
elif ARtag=='p85':
    figname = plotdir + 'fig_S6.land_perspect.AR_freq_change.%s.png'%version_tag


#fig2.savefig(figname, dpi=600)

plt.show()
plt.close()
del(fig2)<jupyter_output><empty_output><jupyter_text>2.2 just warm and cold<jupyter_code>pdf_HIST_cold = np.ma.masked_array(pdf_HIST[0:2], mask=np.isnan(pdf_HIST[0:2])).mean(axis=0)
pdf_HIST_warm = np.ma.masked_array(pdf_HIST[3:8], mask=np.isnan(pdf_HIST[3:8])).mean(axis=0)

pdf_fSST_cold = np.ma.masked_array(pdf_fSST[0:2], mask=np.isnan(pdf_fSST[0:2])).mean(axis=0)
pdf_fSST_warm = np.ma.masked_array(pdf_fSST[3:8], mask=np.isnan(pdf_fSST[3:8])).mean(axis=0)

pdf_subE_cold = np.ma.masked_array(pdf_subE[0:2], mask=np.isnan(pdf_subE[0:2])).mean(axis=0)
pdf_subE_warm = np.ma.masked_array(pdf_subE[3:8], mask=np.isnan(pdf_subE[3:8])).mean(axis=0)

pdf_subF_cold = np.ma.masked_array(pdf_subF[0:2], mask=np.isnan(pdf_subF[0:2])).mean(axis=0)
pdf_subF_warm = np.ma.masked_array(pdf_subF[3:8], mask=np.isnan(pdf_subF[3:8])).mean(axis=0)
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), (11,0), colspan=colspan, rowspan=rowspan)
ax4 = plt.subplot2grid((21,23), (11,11), colspan=colspan, rowspan=rowspan)
axes_full = [ax1, ax2 ,ax3, ax4]

ax1.bar(bin_edges[se:nlen]+1.5, pdf_HIST_cold, width=1, color='royalblue', label='fSST is colder')
ax1.bar(bin_edges[se:nlen]+2.5, pdf_HIST_warm, width=1, color='orchid', label='fSST is warmer')

ax2.bar(bin_edges[se:nlen]+1.5, pdf_fSST_cold-pdf_HIST_cold, width=1, color='royalblue', label='cold')
ax2.bar(bin_edges[se:nlen]+2.5, pdf_fSST_warm-pdf_HIST_warm, width=1, color='orchid', label='warm')

ax3.bar(bin_edges[se:nlen]+1.5, pdf_subE_cold-pdf_HIST_cold, width=1, color='royalblue', label='cold')
ax3.bar(bin_edges[se:nlen]+2.5, pdf_subE_warm-pdf_HIST_warm, width=1, color='orchid', label='warm')

ax4.bar(bin_edges[se:nlen]+1.5, pdf_subF_cold-pdf_HIST_cold, width=1, color='royalblue', label='cold')
ax4.bar(bin_edges[se:nlen]+2.5, pdf_subF_warm-pdf_HIST_warm, width=1, color='orchid', label='warm')

for i in np.arange(8):
    # background
    for axis in axes_full:
        if i%2==0:
            patch = matplotlib.patches.Rectangle((bin_edges[i],-1), 4, 2, color='white', alpha=1, zorder=0)
        else:
            patch = matplotlib.patches.Rectangle((bin_edges[i],-1), 4, 2, color='lightgray', alpha=0.6, zorder=0)
        axis.add_patch(patch)
    
ax1.legend(ncol=1, frameon=False, fontsize=10)

for axis in axes_full:
    axis.plot(np.arange(40), np.zeros(40), linestyle='--', linewidth=0.5, color='black')
    axis.set_xlim([0, 36])
    axis.set_xticks(np.arange(0, 36.1, 8))
    axis.set_xticklabels(np.arange(0, 36.1, 8, dtype=int)*6)
    
ax1.set_ylim([0, 1])
for axis in [ax2, ax3, ax4]:
    if ARtag=='p85':
        axis.set_ylim([-0.075, 0.05])
        axis.set_yticks([-0.075, -0.05, -0.025, 0, 0.025, 0.05])
    elif ARtag=='abs':
        axis.set_ylim([-0.04, 0.02])
        axis.set_yticks(np.arange(-0.04, 0.02+0.00001, 0.02))
    axis.set_ylabel('Change of fraction', fontsize=12)
    
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('Hours of AR per month', fontsize=12)

ax1.set_title('(a) Dist. of AR frequency: CTRL', fontsize=11)
ax1.set_ylabel('Fraction', fontsize=12)
ax2.set_title('(b) Change in dist. (fSST-CTRL)', fontsize=11)
ax3.set_title('(c) Contributions from AR ext.', fontsize=11)
ax4.set_title('(d) Contributions from AR freq.', fontsize=11)


#fig2.savefig(plotdir + 'fig03.land_perspect.AR_freq_change.red_form.%s.png'%version_tag, dpi=600)

plt.show()
plt.close()
del(fig2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16442417331813, 'fraction_numerical': 0.05169137210186241, 'mean_word_length': 4.5023002927645335, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '162744', 'license': 'MIT', 'n_tokens_mistral': 6073, 'n_tokens_neox': 5732, 'n_words': 820, 'path': 'scripts/fig3.figS2.figS6.land.AR_freq_extent_change.ipynb', 'repo_name': 'lucas-uw/Chen-2020-GRL', 'size': '13155'}