Skip to content
Snippets Groups Projects
Commit 66c2f0bf authored by Stefano Serafin's avatar Stefano Serafin
Browse files

more diagnostic plots

parent 0b45bdd8
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 123 additions and 4 deletions
PE_CBL.py
+ 123
− 4
View file @ 66c2f0bf
......@@ -56,7 +56,7 @@ if __name__ == '__main__':
fig05 = False # CHECKED OK # Experiments with systematic model errors
fig06 = False # CHECKED OK # Analysis increments in model error experiments
fig07 = False # CHECKED OK # Experiment with error growth
fig08 = False # CHECKED OK # Perfect-obs experiments
fig08 = True # CHECKED OK # Perfect-obs experiments
fig09 = False # CHECKED OK # DA diagnostics (profiles) for perfect-obs experiments
fig10 = False # CHECKED OK # Comparison EAKF vs LETKF
......@@ -69,11 +69,13 @@ if __name__ == '__main__':
opt06 = False # CHECKED OK # Check impact of RTPS variance inflation factor for parameters
opt07 = False # CHECKED OK # Check impact of RTPS variance inflation factor for state vector
opt08 = False # CHECKED OK # Check sensitivity to observation sorting
opt09 = False # CHECKED OK # Plot covariance matrix
opt10 = False # CHECKED OK # Consistency checks
opt09 = False # CHECKED OK # Background error covariance matrix
opt10 = False # CHECKED OK # Profiles of consistency diagnostics
opt11 = False # CHECKED OK # An experiment without any variance inflation
opt12 = False # CHECKED OK # Old version of Figure 08
opt13 = False # CHECKED OK # More profiles of consistency diagnostics
opt14 = False # CHECKED OK # More plots of consistency checks
opt15 = False # CHECKED OK # O-B for experiments C1 and C2
# A test to check if results are exactly reproducible
check_reproducibility = False
......@@ -361,7 +363,7 @@ if __name__ == '__main__':
ax2.set_title(r'a) Exp. B$_4$, $\delta\overline{y}$ (K)')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
ax3.set_title(r'b) Exp. B$_4$, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
ax3.set_title(r'b) Exp. B$_4$, $\delta\overline{p}\prime\prime$')
ax3.set_xlabel('Time (h)')
ax3.set_ylabel('Height (m)')
maxtime = exp_B4.trun
......@@ -776,6 +778,123 @@ if __name__ == '__main__':
fig.savefig('diagnostics_profiles_errors_vs_deviations.png',format='png',dpi=300)
p.close(fig)
if opt14:
workdir = "./runs/seed=181612_enssize=20_EAKF_6hrs/"
[exp_I1,exp_I2,exp_I3,exp_I4,
exp_I5,exp_I6,exp_I7,exp_I8,] = map(get_experiment,
[workdir+'/da_exp_I1.json',
workdir+'/da_exp_I2.json',
workdir+'/da_exp_A1.json',
workdir+'/da_exp_I4.json',
workdir+'/da_exp_I5.json',
workdir+'/da_exp_A1.json',
workdir+'/da_exp_I7.json',
workdir+'/da_exp_I8.json'])
exp_I3.label = 'I3'
exp_I6.label = 'I6'
# Plots: Scatter plots of Desroziers statistics
experiments1 = [exp_I1, exp_I2, exp_I3, exp_I4]
experiments2 = [exp_I5, exp_I6, exp_I7, exp_I8]
labels1 = [exp.label for exp in experiments1]
labels2 = [exp.label for exp in experiments2]
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
ax1 = plot_desroziers([i.dg for i in experiments1],labels=labels1,ax=ax1)
ax2 = plot_desroziers([i.dg for i in experiments2],labels=labels2,ax=ax2)
ax1.set_xlim([0,0.018])
ax1.set_ylim([0,0.018])
ax2.set_xlim([0,0.018])
ax2.set_ylim([0,0.018])
ax1.set_title(r"a) Variable $\alpha_{RTPS}$ for parameter")
ax2.set_title(r"b) Variable $\alpha_{RTPS}$ for state")
fig.savefig('Diag_desroziers_1.png',format='png',dpi=150)
p.close(fig)
[exp_L1,exp_L2,exp_L3,exp_L4] = map(get_experiment,
[workdir+'/da_exp_L1.json',
workdir+'/da_exp_A1.json',
workdir+'/da_exp_L3.json',
workdir+'/da_exp_L4.json'])
exp_L2.label = 'L2'
experiments3 = [exp_L1, exp_L2, exp_L3, exp_L4]
labels3 = [exp.label for exp in experiments3]
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
ax1 = plot_desroziers([i.dg for i in experiments3],labels=labels3,ax=ax1)
ax2 = plot_desroziers([i.dg for i in experiments2],labels=labels2,ax=ax2)
ax1.set_xlim([0,0.018])
ax1.set_ylim([0,0.018])
ax2.set_xlim([0,0.018])
ax2.set_ylim([0,0.018])
ax1.set_title(r"a) Variable localicazion scale")
ax2.set_title(r"b) Variable $\alpha_{RTPS}$ for state")
fig.savefig('Diag_desroziers_2.png',format='png',dpi=150)
p.close(fig)
if opt15:
[exp_C1,exp_C2] = map(get_experiment, [workdir+'/da_exp_C1.json',
workdir+'/da_exp_C2.json'])
ntimes,nobs = exp_C1.da.observations.shape
# Expand times array to 2D (matching the number of observations)
times = exp_C1.da.time/3600
times = times[:,None]+np.zeros((ntimes,nobs))
# Sort coordinate and innovation arrays in regular order
ominb_C1 = np.zeros((ntimes,nobs))+np.nan
ominb_C2 = np.zeros((ntimes,nobs))+np.nan
z_C1 = np.zeros((ntimes,nobs))+np.nan
z_C2 = np.zeros((ntimes,nobs))+np.nan
for [z_sorted,innov_sorted,z,innov] in zip([z_C1,z_C2],[ominb_C1,ominb_C2],[exp_C1.da.obs_coordinates,exp_C2.da.obs_coordinates],[exp_C1.da.innov,exp_C2.da.innov]):
for k in range(ntimes):
if is_sorted(z[k,:]):
indices = np.arange(nobs)
else:
indices = np.argsort(z[k,:])
innov_sorted[k,:] = innov[k,indices]
z_sorted[k,:] = z[k,indices]
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
c1 = ax1.pcolormesh(times,z_C1,ominb_C1,vmin=-0.5,vmax=0.5,cmap='RdBu_r')
c2 = ax2.pcolormesh(times,z_C2,ominb_C2,vmin=-0.5,vmax=0.5,cmap='RdBu_r')
ax1.text(0.1,0.2,r'$\mu_d$ = %05.4f'%ominb_C1.mean(), transform = ax1.transAxes)
ax1.text(0.1,0.1,r'$\sigma_d$ = %05.4f'%ominb_C1.std(),transform = ax1.transAxes)
ax2.text(0.1,0.2,r'$\mu_d$ = %05.4f'%ominb_C2.mean(), transform = ax2.transAxes)
ax2.text(0.1,0.1,r'$\sigma_d$ = %05.4f'%ominb_C2.std(),transform = ax2.transAxes)
#
ax1.set_ylim([0,2000])
ax2.set_ylim([0,2000])
ax1.set_title(r'a) Exp. C$_1$, $\overline{\theta}^o-\overline{\theta}^b$')
ax2.set_title(r'b) Exp. C$_2$, $\overline{\theta}^o-\overline{\theta}^b$')
ax1.set_xlabel('Time (h)')
ax2.set_xlabel('Time (h)')
ax1.set_ylabel('Height (m)')
ax2.sharey(ax1)
maxtime = exp_C1.trun
ax1.set_xlim([0,maxtime/3600])
ax2.set_xlim([0,maxtime/3600])
ax1.set_xticks(np.arange((maxtime+1)/3600))
ax2.set_xticks(np.arange((maxtime+1)/3600))
p.setp(ax2.get_yticklabels(), visible=False)
#
p.colorbar(c1,orientation='horizontal')
p.colorbar(c2,orientation='horizontal')
#
fig.savefig('ominb.png',format='png',dpi=300)
p.close(fig)
if check_reproducibility:
def get_obs(f):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment