Modelling Local Convective Extreme Events - Does Increasing the Resolution Help?

In [92]:
HTML(showcode)
Out[92]:

2.1 The diurnal cycle in different Monsoon stages

In [19]:
#Plot the diurnal Cycle
fig, ax = plt.subplots(1, figsize=(12,6))
ax.set_ylabel('Rain-rate mm/d', fontsize=22)
ax.set_xlabel('Local-time', fontsize=22)
ax.plot(idx, dc_burst.values, 'fuchsia', linewidth=2.0, label='Bursts')
ax.plot(idx, dc_break.values, 'g', linewidth=2.0, label='Breaks')
ax.plot(idx, (dc_break+dc_burst)/2, 'firebrick', linewidth=2.0, label='All')
ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%_H'))
ax.legend(loc=0)
ax.grid(linewidth=0)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)

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.

In [24]:
#Plot area avg TS
import matplotlib.dates as mdates
myFmt = mdates.DateFormatter('%d/%m %H LT')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(20,10), dpi=72)
ax = fig.add_subplot(111)
obs_time_1 = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
obs_time_2 = pd.DatetimeIndex(OBS.dataset[2].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
um044_time = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
um133_time = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
## First plot the area avg of the cpol data
#ax.plot(OBS.dataset[0].coords['t'], np.nanmean(OBS.dataset[0].variables['lsrain'].values[:,0,:]*mask.filled(np.nan),axis=(1,2)), 
#        label='CPOL', color='firebrick')
#ax.plot(obs_time_2, OBS.dataset[2].variables['lsrain'].mean(axis=(1,2)), label='Cmorph 8km',
#        color='fuchsia', linestyle='-')
ax.plot(obs_time_1, OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)), linestyle='-',
        color='lightseagreen', label='CPOL')
#############
ax.plot(um133_time, UM133ens.variables['lsrain'].mean(axis=(0,2,3,4)), color='fuchsia',
        label='UM 1.33km ensmean', linestyle='--')
#ax.fill_between(UM133ens.coords['t'].values, 
#                UM133ens.variables['lsrain'].mean(axis=(2,3,4)).min(axis=0).values,
#                UM133ens.variables['lsrain'].mean(axis=(2,3,4)).max(axis=0).values, color='fuchsia',
#                alpha=0.2)
############
ax.plot(um044_time, UM044ens.variables['lsrain'].mean(axis=(0,2,3,4)), color='fuchsia',
        label='UM 0.44km ensmean')
ax.fill_between(um044_time, 
                UM044ens.variables['lsrain'].mean(axis=(2,3,4)).min(axis=0).values,
                UM044ens.variables['lsrain'].mean(axis=(2,3,4)).max(axis=0).values, color='cornflowerblue',
                alpha=0.2)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
ax.set_ylim(0.0,0.42)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=18)
ax.xaxis.set_major_formatter(myFmt)
ax.grid(color='r', linestyle='-', linewidth=0)

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

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

Maps of Rainfall (0.44km)

In [27]:
#Create the avg. maps of rainfall 0.44km
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.35, top=0.95,left=0.05, hspace=0, wspace=0)
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)
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)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
             resolution='f', area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m[0].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 = ['Obs']
ax[-1].set_title('Obs', fontsize=24)
for i in range(1,9):
    ax.append(fig.add_subplot(3,3,i+1))
    m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
                     resolution='f', area_thresh=100, ax=ax[-1]))
    m[-1].drawcoastlines()
            #m[-1].shadedrelief()
    im.append(m[-1].pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
    tit.append(datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime('%e %b %H UTC'))
    ax[-1].set_title(tit[-1],fontsize=24)

Maps of Rainfall (1.33km)

In [28]:
#Create the avg maps of rainfall for 1.33km
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.35, top=0.95,left=0.05, hspace=0, wspace=0)
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)
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)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
             area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m[0].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 = ['Obs']
ax[-1].set_title('Obs', fontsize=24)
for i in range(1,9):
    ax.append(fig.add_subplot(3,3,i+1))
    m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
                     area_thresh=100, ax=ax[-1]))
    m[-1].drawcoastlines()
            #m[-1].shadedrelief()
    im.append(m[-1].pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
    tit.append(datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime('%e %b %H UTC'))
    ax[-1].set_title(tit[-1],fontsize=24)

Ensemble mean and std

In [29]:
#Plet the avg maps for comparison
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.01, top=0.6,left=0.01, hspace=0, wspace=0)
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 = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
             area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m[0].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.43, 5:0.23, 8:0.03, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=24)
cbar_ax = [fig.add_axes([0.95, height[2], 0.01, 0.15])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=24)
cbar[-1].set_label('Rain-rate [mm/d]',size=24)
for i in range(1,len(data)):
    if data[i] is not None:
        ax.append(fig.add_subplot(3,3,i+1))
        m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), 
                         resolution='f', area_thresh=100, ax=ax[-1]))
        m[-1].drawcoastlines()
            #m[-1].shadedrelief()
        im.append(m[-1].pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
        ax[-1].set_title(tit[i],fontsize=24)
        if i in (5, 8, 11):
            cbar_ax.append(fig.add_axes([0.95, height[i], 0.01, 0.15]))
            cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
            cbar[-1].ax.tick_params(labelsize=24)
            cbar[-1].set_label('Rain-rate [mm/d]',size=24)

2.1.2 The Simulated Diurnal Cycle

In [25]:
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(20,10), dpi=72)
ax = fig.add_subplot(111)
####################
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)
ax.plot(obs_S.index, obs_S.values, linestyle='-', color='lightseagreen', label='CPOL')
##################
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).mean(dim='ens')
um_2 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).min(dim='ens')
um_3 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).max(dim='ens')
um_044 = pd.DataFrame({'mean':um_1.variables['lsrain'].values,
                       'min' :um_2.variables['lsrain'].values,
                       'max' :um_3.variables['lsrain'].values},
                       index=(um_1.variables['hour'].values + 9) % 24).sort_index(inplace=False)
ax.plot(um_044.index, um_044['mean'].values, color='fuchsia', linestyle='-', label='UM 0.44km ensmean')
###################
um_4 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).mean(dim='ens')
um_5 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).min(dim='ens')
um_6 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).max(dim='ens')
um_133 = pd.DataFrame({'mean':um_4.variables['lsrain'].values,
                       'min' :um_5.variables['lsrain'].values,
                       'max' :um_6.variables['lsrain'].values},
                       index=(um_4.variables['hour'].values + 9) % 24).sort_index(inplace=False)

ax.plot(um_133['mean'].index, um_133['mean'].values, color='fuchsia', linestyle='--', label='UM 1.33km ensmean')
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)

#####################
ax.fill_between(um_044.index, um_044['min'], um_044['max'], color='cornflowerblue', alpha=0.2, label='UM 0.44km ens. spread')
ax.fill_between(um_133.index, um_133['min'], um_133['max'], color='coral', alpha=0.2, label='UM 1.33km ens. spread')
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=18)
ax.set_xlabel('Local time [h]', fontsize=18)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
ax.grid(color='r', linestyle='-', linewidth=0)
In [55]:
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Diurnal-2.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))
Out[55]:

Location and timing of daily rainfall peaks

In [48]:
#Plot the maps
plt.close()
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=(25,15))
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.67,left=0.01, hspace=0, wspace=0)
first = True
val_max = 3
colm = colmap2.Blues
colm.set_under('w', alpha=0)
m, ax, im = [], [], []

for i in range(24):
    break
    fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
    if first:
        for col in range(3*2):
            if col != 3:
                ax.append(fig.add_subplot(2,3,col+1))
                m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), 
                                 resolution='f', area_thresh=1, ax=ax[-1]))
                m[-1].drawcoastlines()
    
        ###MEAN
        ax[0].set_title('CPOL',fontsize=22)
        ax[1].set_title('UM 1.33km mean', fontsize=22)
        ax[2].set_title('UM 0.44km mean', fontsize=22)
        ax[-2].set_title('UM 1.33km std', fontsize=22)
        ax[-1].set_title('UM 0.44km std', fontsize=22)
        im.append(m[0].pcolormesh(lon1, lat1, obs['lsrain'].values[i],  shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
        im.append(m[1].pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
        im.append(m[2].pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
        ###STD
        im.append(m[-2].pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
        im.append(m[-1].pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
        first = False
        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=24)
        
    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)%24),size=24)
    fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
    #break

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 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9) % 24 / 2).round(0) * 2 
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9) % 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)

m, ax = [], []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km ens. mean', 'UM 0.44km ens. mean']
fig = plt.figure(figsize=(25,15))
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.67,left=0.01, hspace=0, wspace=0)
cbar_ax = fig.add_axes([0.12, 0.17, 0.74, 0.02])

for i in range(3):
    ax.append(fig.add_subplot(1,3,i+1))
    
    m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), 
            resolution='f', area_thresh=1, ax=ax[-1]))
    m[-1].drawcoastlines()
    try:
        im = m[-1].pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
    except TypeError:
        im = m[-1].pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
    ax[-1].set_title(names[i], fontsize=22)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Local Time [h]', size=24)

3.1 Pdf's of Rainfall

In [62]:
#Plot the pdfs
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.UM044['shape']*3.4], loc = np.log(LogNorm.UM044['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.UM044['popt']/4, PowerFit.UM044['pcov']
fig = plt.figure()
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um044vec.loc[um044vec>0], um133vec.loc[um133vec>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,1.4)
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 = ax2.hist(histdata, bins=nbins, color=facecolors, normed=1, alpha=0.25)
#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')
fig.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.legend(loc=4, fontsize=28)
ax.tick_params(labelsize=28)

ax.set_ylim(0.,.4)
ax.set_xlim(0.,40)
ax.set_ylabel('Density', fontsize=28)
ax.set_xlabel('Rain-rate mm/h', fontsize=28)
axin = inset_axes(ax, width = "60%", height = "60%", loc=1)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#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.plot(np.linspace(0.01,80,200), dist.pdf(np.linspace(0.01, 80,200)), lw=4,  color='firebrick', alpha=0.8,
          label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=24)
axin.set_xlim(0.01, 90)
axin.set_ylim(0.000025, 1)
axin.grid(color='r', linestyle='-', linewidth=0)
ann = axin.annotate("b)", (0.01,0.03), xycoords="axes fraction", fontsize=28)
ax.grid(color='r', linestyle='-', linewidth=0)
axin.plot(np.linspace(0.01, 50,200), fitpdf(np.linspace(0.01,50,200),popt[0]), lw=4, color='navy')
axin.set_ylabel('$\\log$(Density)', fontsize=24)
axin.set_xlabel('$\\log$(Rain-rate)', fontsize=24)
Out[62]:
<matplotlib.text.Text at 0x7f4c2f438ef0>

3.2 Percentiles

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

In [60]:
#Plot the percentages
fig=plt.figure()
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.plot(PERC.index, PERC['Obs'].values, color='lightseagreen', linestyle='-', label='CPOL',lw=3)
ax.plot(PERC.index, PERC['UM 1.33km'].values, color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], color='lightseagreen', linestyle='-', label='CPOL',lw=3)
axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=24)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])

axin.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
#axin.get_xaxis().get_major_formatter().labelOnlyBase = False
axin.grid(color='r', linestyle='-', linewidth=0)
axin.annotate("b)", (0.01,0.03), xycoords="axes fraction", fontsize=28)
axin.set_xticks(range(50,100,10), minor=False)
ax.grid(color='r', linestyle='-', linewidth=0)
axin.set_ylabel('Rain-rate[mm/h]', fontsize=24)
axin.set_xlabel('Percentile [%]', fontsize=24)
Out[60]:
<matplotlib.text.Text at 0x7f4c2fcfc748>