Modeling Extreme Events

HTML(showcode)

2.1 The diurnal cycle in different Monsoon stages

Burst vs break

2.1 A Week of Hector

Now we compare the ensemble timeseries of the simulated rainfall for the 1.33km and 0.44km Simulations with the CPOL observations.

time, obs, rain = [], [], []
for e in range(len(ensembles)):
    r1 = list(UM133ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
    r2 = list(UM044ens.variables['lsrain'][e].mean(axis=(1,2,3)).values)
    obs+=len(r1)*['UM 1.33km'] + len(r2)*['UM 0.44km']
    t1 = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
    rain += r1 + r2
    t2 = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
    time += list(t1) + list(t2)

r3 = list(OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)).values)
t3 = list(pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime())
o3 = len(r3)*['CPOL']
rain += r3
obs += o3
time += t3
plot_data = pd.DataFrame({'rain-rate':rain, 'time':time, 'Type':obs})

#Plot area avg TS
import matplotlib.dates as mdates
mpld3.disable_notebook()
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
myFmt = mdates.DateFormatter('%d/%m')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
sns.lineplot(x='time', y='rain-rate', hue='Type', ax=ax, data=plot_data, lw=1.3)

plt.subplots_adjust(left=0.1, bottom=0.1, right=0.98, top=0.9)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
#ax.set_xlim(plot_data.time.iloc[0], plot_data.time.iloc[-1])
ax.set_ylim(-0.001,0.5)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
_ = ax.xaxis.set_major_formatter(myFmt)
_ = ax.set_xlabel('')
#_ = ax.set_title('Comparison of Area-Avg Rainfall', fontsize=fontsize)
fig.subplots_adjust(bottom=0.4, top=0.93, right=0.83, left=0.1)
mpld3.display(fig)
#nplot+=1; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

mpld3

open('../slides/figure.html','w').write(d3plot.data)

790506

2.1.1 Maps of Rainfall

Now lets look at some animations of UM rainfall occuring over the Tiwi-islands and compare it to the cpol observations

os.chdir('/home/unimelb.edu.au/mbergemann/Work/Programming/Extremes/python/')

#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Ens-3.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" loop="true" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))

Maps of Rainfall (0.44km)

#Create the avg. maps of rainfall 0.44km
#m1, m2 = None, None
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
rom = {1:'i', 2:'ii', 3:'iii', 4:'iv', 5:'v', 6:'vi', 7:'vii', 8:'viii'}
fig = plt.figure(figsize=(15,12))

UM = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues

colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m2 is None:
    m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
             resolution='f', area_thresh=1)
if m1 is None:
    m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1),
             resolution='f', area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
    if m2 is None:
        m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
                         resolution='f', area_thresh=1)
    ax.append(fig.add_subplot(3,3,i+1))
    m2.drawcoastlines(linewidth=0.8)
    im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
    tit.append('Ens. %s'%rom[i].upper())
    ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

Maps of Rainfall (1.33km)

#Create the avg maps of rainfall for 1.33km
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':True})
fig = plt.figure(figsize=(15,12))
UM = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
if m1 is None:
    m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1), resolution='f',
                 area_thresh=1)
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m2.drawcoastlines(linewidth=0.8)
cbar_ax = fig.add_axes([0.12, 0.12, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-Rate [mm/d]',size=24)
tit = ['CPOL']
ax[-1].set_title('CPOL', fontsize=fontsize)
for i in range(1,9):
    if m2 is None:
        m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
                     area_thresh=1)
    ax.append(fig.add_subplot(3,3,i+1))
   
    m2.drawcoastlines(linewidth=0.8)
            #m[-1].shadedrelief()
    im.append(m2.pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
    tit.append('Ens. %s'%rom[i].upper())
    ax[-1].set_title(tit[-1],fontsize=fontsize)
fig.subplots_adjust(right=0.99, bottom=0.15, top=0.95,left=0.01, hspace=0.2, wspace=0.0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

Ensemble mean and std

#Plet the avg maps for comparison
fig = plt.figure(figsize=(18,12))
#sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
UM_2 = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
UM_1 = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01)
data = [obs_data[1:-1]/6,
        np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
        np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
        None,
        np.ma.masked_less(UM_1.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0],
        np.ma.masked_less(UM_2.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0]
       ]

data.append(None)

data.append(data[0] - data[1])
data.append(data[0] - data[2])
UM_1min = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_2min = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).min(axis=0)[0]

UM_1max = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
UM_2max = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).max(axis=0)[0]

#data.append(None)
#data.append(UM_1max - UM_1min)
#data.append(UM_1max - UM_1min)


lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM_1.coords['lat'][:], UM_1.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m2.pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m2.drawcoastlines()
#cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
#cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
#cbar.ax.tick_params(labelsize=24)
#cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['CPOL avg.', 'UM 0.44km ens avg.', 'UM 1.33km ens avg.',
       None, 'UM 0.44 ens std.', 'UM 1.33km ens std.',
       None, 'CPOL - UM 0.44 ens avg.' ,'CPOL - UM 1.33km ens avg.',
       None, 'UM 0.44 ens range', 'UM 1.33km ens range']
colors = {1:colm, 2:colm, 4:colm, 5:colm,
          7:colmap.GMT_polar_r, 8:colmap.GMT_polar_r,
         10:colm, 11:colm}
Range = {1:(0.0, 25/6), 2:(0.0, 25/6), 4:(0,8/6), 5:(0,8/6), 7:(-20/6,20/6), 8:(-20/6,20/6), 10:(0,30/6), 11:(0,30/6)}
height = {2:0.682, 5:0.37, 8:0.044, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=fontsize-2)
cbar_ax = [fig.add_axes([0.91, height[2], 0.01, 0.23])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=fontsize-4)
cbar[-1].set_label('Rain-Rate [mm/d]',size=fontsize-4)
for i in range(1,len(data)):
    if data[i] is not None:
        ax.append(fig.add_subplot(3,3,i+1))
            #m[-1].shadedrelief()
        im.append(m2.pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
        m2.drawcoastlines()
        ax[-1].set_title(tit[i],fontsize=fontsize-2)
        if i in (5, 8, 11):
            cbar_ax.append(fig.add_axes([0.91, height[i], 0.01, 0.23]))
            cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
            cbar[-1].ax.tick_params(labelsize=fontsize-4)
            cbar[-1].set_label('Rain-rate [mm/d]',size=fontsize-2)
fig.subplots_adjust(right=0.9, bottom=0.01, top=0.95,left=0.01, hspace=0.0, wspace=0)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%i03.png'%nplot, bbox_inches='tight', format='png', dpi=300)

2.1.2 The Simulated Diurnal Cycle

data, time, model = [], [], []
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
for i in range(len(ensembles)):
    d1 = UM133.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
    d2 = UM044.dataset[i].groupby('t.hour').mean(dim='t').mean(dim=('surface','lat','lon'))
    t1 = pd.Series(d1.variables['lsrain'].values, index=(d1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
    t2 = pd.Series(d2.variables['lsrain'].values, index=(d2.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
    time += list(t1.index) + list(t2.index)
    data += list(t1.values) + list(t2.values)
    model += len(t1)*['UM 1.33km'] + len(t2)*['UM 0.44km']
cpol_1 = OBS.dataset[1].groupby('t.hour').mean(dim='t').mean(dim=('y','x'))
obs_S = pd.Series(cpol_1.variables['lsrain'].values, index=(cpol_1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
data += list(obs_S.values)
time += list(obs_S.index)
model += len(obs_S)*['Cpol']
plot_data = pd.DataFrame(dict(data=data, time=time, Type=model))
#nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#mpld3.enable_notebook()

Location and timing of daily rainfall peaks

#Plot the maps
plt.close()
mpld3.disable_notebook()
sns.set_style("whitegrid", {'axes.grid':True, 'ticks':False})
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
um_2 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
obs  = OBS.dataset[1].groupby('t.hour').mean(dim='t')

time = (um_1.hour + 9) % 24

data = [obs, um_2, um_1]

fig = plt.figure(figsize=(15,7))
lat1, lon1 = OBS.dataset[1].variables['latitude'].values[:,0], OBS.dataset[1].variables['longitude'].values[0,:]
lat2, lon2 = um_1.coords['lat'][:], um_1.coords['lon'][:]
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.98,left=0.01, hspace=0, wspace=0)
first = True
val_max = 2
colm = colmap2.Blues
colm.set_under('w', alpha=0)
colm.set_bad('w', alpha=0)

#ax = [fig.add_subplot(2,3,1)]
mm, im = [], []
col = 2
for i in range(24):
    #break
    fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
    if first:
        for col in range(1, 7):
            if col != 4:
                ax1 = fig.add_subplot(2,3,col)
                m2.drawcoastlines(linewidth=0.8)
            if col == 1:
                ax1.set_title('CPOL',fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon1, lat1, obs['lsrain'].values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 2:
                ax1.set_title('UM 1.33km mean', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 3:
                ax1.set_title('UM 0.44km mean', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 5:
                ax1.set_title('UM 1.33km std', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
            elif col == 6:
                ax1.set_title('UM 0.44km std', fontsize=fontsize)
                _ = im.append(m2.pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))      
         
                cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.02])
                cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
                _ = cbar.ax.tick_params(labelsize=fontsize-4)
        first = False
    else:
        ###MEAN
        im[0].set_array(obs['lsrain'].values[i].ravel())
        im[1].set_array(um_2['lsrain'].mean(dim='ens').values[i].ravel())
        im[2].set_array(um_1['lsrain'].mean(dim='ens').values[i].ravel())
        ###STD
        im[-2].set_array(um_2['lsrain'].std(dim='ens').values[i].ravel())
        im[-1].set_array(um_1['lsrain'].std(dim='ens').values[i].ravel())
    cbar.set_label('Rain-Rate at %2i LT [mm/h]'%((i+9.5)%24),size=fontsize)
    fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
    #break
dest_dir = os.path.abspath('Vid')
make_mp4_from_frames(outdir, dest_dir, 'WeekOfHector-Diurnal-2.mp4', 4, glob='DiurnalCycle_??')
fig.clf()
plt.close()
day = (time[1:5] - 9) % 24
night = (time[5:5+4] - 9) % 24
obs_diff = obs['lsrain'][day].sum(dim='hour') - obs['lsrain'][night].sum('hour')
um_1_diff = um_1['lsrain'][day].sum(dim='hour') - um_1['lsrain'][night].sum('hour')
um_2_diff = um_2['lsrain'][day].sum(dim='hour') - um_2['lsrain'][night].sum('hour')

mask = np.ma.masked_invalid(obs['lsrain'].values) * 0 + 1
obs_max = (((np.argmax(obs['lsrain'].values * mask, axis=0)*mask[0] + 9.5 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2 
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9.5) % 24 / 2).round(0) * 2
colm = CubeYF_20.get_mpl_colormap(N=24, gamma=2.0)
#colm = qualitative.Paired[12].get_mpl_colormap(N=24, gamma=2.0)
colm.set_bad('w')
#m, ax = [], []
ax = []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km', 'UM 0.44km']
gs = gridspec.GridSpec(30, 3)

fig = plt.figure(figsize=(15,8))
fig.subplots_adjust(right=0.98, bottom=0.03, top=0.98,left=0.01, hspace=0, wspace=0)
cbar_ax = plt.subplot(gs[13, :])

for i in range(3):
    ax.append(plt.subplot(gs[:13, i]))
    if m2 is None:
        m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
                     area_thresh=1)
    
    m2.drawcoastlines()
    try:
        im = m2.pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
    except TypeError:
        im = m2.pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
    ax[-1].set_title(names[i], fontsize=fontsize)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=fontsize-4)
#cbar.set_label('Local Time [h]', size=fontsize)
mpld3.disable_notebook()
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
#fig = plt.figure(figsize=(15,8), dpi=72)
#ax = fig.add_subplot(111)
sns.set_style("darkgrid", {'axes.grid':True, 'ticks':True})
ax = plt.subplot(gs[18:, :])

####################
sns.lineplot(x=plot_data.time, y=plot_data.data, hue=plot_data.Type, ax=ax, lw=3, data=plot_data)
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-Rate [mm/h]', fontsize=fontsize)
ax.set_xlabel('Local Time [h]', fontsize=fontsize)
ax.legend(loc=0, fontsize=fontsize-4)
ax.tick_params(labelsize=fontsize-4)
fig.subplots_adjust(bottom=0.145, top=0.99, right=0.99, left=0.1)
#ax.grid(color='r', linestyle='-', linewidth=0)
mpld3.disable_notebook()
nplot = 12; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

3.1 Pdf's of Rainfall

#Plot the pdfs
sns.set_style('darkgrid')
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.CPOL['shape']*3.4], loc = np.log(LogNorm.CPOL['scale']*1.9))#, shape  # mu, sigma
med = np.median(um044vec.loc[um044vec>0].values)
X = np.linspace(0,65,200)
#b = 2.5
nbins = 100
def fitpdf(x, b):
    a=med*100
    return ((b/a) * (x/a)**(b-1)) / (1+(x/a)**b)**2
#popt, pcov = curve_fit(powerlaw, um044cnt['Rain-rate'].values[1:],
#                       um044cnt['count'].values[1:])
popt, pcov = PowerFit.CPOL['popt']/4, PowerFit.CPOL['pcov']
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um133vec.loc[um133vec>0], um044vec.loc[um044vec>0], obs_vec.loc[obs_vec>0])
#facecolors = ['fuchsia', 'firebrick', 'lightseagreen']
#n1_1, bins1_1, patches1_1 = ax2.hist(, bins=nbins, normed=1, facecolor='firebrick', alpha=0.25)
#n2, bins2, patches2 = ax2.hist(, bins=nbins, normed=1, facecolor='lightseagreen', alpha=0.25)
ax2.set_ylim(0,0.9)
ax2.set_yticks([])
#ax.plot(np.linspace(0.01, 100,200), fitpdf(np.linspace(0.01,100,200),popt[0]), lw=4, color='navy',
#        label='Log-logistic $\\beta = %02.2e$'%popt[0])
#ax.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
#ax.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#ax.scatter(um044cnt['Rain-rate'].values, um044cnt['count'].values, c='red', s=125, alpha=0.24)
#ax.scatter(obscnt['Rain-rate'].values, obscnt['count'].values, c='purple', s=125, alpha=0.24)
#ax.plot(np.linspace(0, 50,200), um044pdf, lw=4, label='UM 0.44km ens', color='lightseagreen')
#ax.plot(np.linspace(0, 50,200), obspdf, lw=4, label='CPOL', color='fuchsia')
#ax.plot(np.linspace(0.0,50,200), dist.pdf(np.linspace(0.0, 50,200)), color='firebrick', alpha=0.8,
#        lw=4, label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
n1, bins1, patches1 = ax.hist(histdata, bins=nbins, normed=1, alpha=0.55)
#fitted_pl.plot_pdf(ax=ax,color='k',lw=4)
#ax.plot(np.linspace(0.1, 100), um044dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 0.44km area-avg')
#ax.plot(np.linspace(0.1, 100), um133dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 1.33km area-avg')


ax.tick_params(labelsize=fontsize-4)

#ax.set_ylim(0.,.25)
ax.set_xlim(0.,20)
ax.set_ylabel('Density', fontsize=fontsize)
ax.set_xlabel('Rain-Rate [mm/h]', fontsize=fontsize)
axin = inset_axes(ax, width = "80%", height = "80%", loc=1)

axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'-', label='UM 1.33 km ens',lw=2)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=2)
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=2)
#axin.plot(np.linspace(0.1,80,200), dist.pdf(np.linspace(0.1, 80,200)), lw=1, alpha=0.8,
#          label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
#axin.set_xlabel('$\\log$(Rain-rate)', fontsize=fontsize-2)
#axin.plot(X, um044pdf,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#axin.plot(X, obspdf,'-', label='CPOL',lw=4, color='lightseagreen')

axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(0.1, 100)
axin.set_ylim(0.000023, 1)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
axin.legend(loc=3, fontsize=fontsize-4)
#n.set_ylabel('$\\log$(Density)', fontsize=fontsize-2)

fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)
#mpld3.enable_notebook()

3.2 Percentiles

Now get the percentiles for the ensemble simulation and compare it against the observations

#Plot the percentages
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)

ax.plot(PERC.index, PERC['UM 1.33km'].values, label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, label='UM 0.44km', lw=3)
ax.plot(PERC.index, PERC['Obs'].values, label='CPOL',lw=2)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile []', fontsize=fontsize)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=fontsize)

ax.tick_params(labelsize=fontsize-4)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)

axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], label='UM 0.44km', lw=3)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], label='CPOL',lw=2)

#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=fontsize-4)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])
axin.legend(loc=0, fontsize=fontsize-4)
#axin.get_xaxis().set_major_formatter(mpl.ticker.ScalarFormatter(useMathText=False))
#axin.get_xaxis().get_major_formatter().labelOnlyBase = True
#axin.set_xticklabels([], rotation=0)
ann = axin.annotate("(b)", (0.95,0.93), xycoords="axes fraction", fontsize=fontsize)
#axin.set_xticks(range(50,100,10), minor=False)
axin.set_yscale('log')
axin.set_xscale('log')
fig.subplots_adjust(bottom=0.13, right=0.98, top=0.97)
nplot+=1 ; _ = fig.savefig('Vid/Fig_%03i.png'%nplot, bbox_inches='tight', format='png', dpi=300)

nplot

15

mpl.ticker.ScalarFormatter?

2.1 Distribution of Percentiles during break and bursts

Now we compare the distributions of percentiles during break and burst periods. We also add

Created with Jupyter, delivered by Fastly, rendered by Rackspace.