diff --git a/graphics.py b/graphics.py
index 85ef48eca512f3ef8852086102a637957aaa9122..aaefb23353f8a8ab3dcf6ce107aa9394019fc8b6 100644
--- a/graphics.py
+++ b/graphics.py
@@ -3,7 +3,7 @@
import numpy as np
from matplotlib import pyplot as p
import matplotlib.colors as colors
-from ENDA import experiment
+from ENDA import compute_distances, gaspari_cohn
is_sorted = lambda a: np.all(a[:-1] <= a[1:])
@@ -344,6 +344,7 @@ def plot_CBL_identifiability(cbl,sigma_o,ax=None):
C = deltap.T
D = ( (sigma2_p/sigma2_yb)**0.5 ).T
E = ( innov*sigma2_yb/(sigma2_yb+sigma_o**2) ).T
+ F = innov.T
# For plotting
contours = np.linspace(cbl.nr.theta_0,cbl.nr.theta_0+cbl.nr.gamma*zmax,ncont)
@@ -370,7 +371,8 @@ def plot_CBL_identifiability(cbl,sigma_o,ax=None):
elif type(ax) is list and len(ax)==2:
- c0 = ax[0].pcolormesh(times/3600.,zobs,D,cmap=new_cmap)
+ #c0 = ax[0].pcolormesh(times/3600.,zobs,D,cmap=new_cmap) # For old version of Fig 8
+ c0 = ax[0].pcolormesh(times/3600.,zobs,B,norm=colors.CenteredNorm(),cmap='RdBu_r')
ax[0].contour(times/3600.,zobs,theta.mean(axis=2).T,
contours,
colors='black',
@@ -378,7 +380,8 @@ def plot_CBL_identifiability(cbl,sigma_o,ax=None):
linewidths=0.75)
ax[0].set_ylim([0,zmax])
- c1 = ax[1].pcolormesh(times/3600.,zobs,E,norm=colors.CenteredNorm(),cmap='RdBu_r')
+ #c1 = ax[1].pcolormesh(times/3600.,zobs,E,norm=colors.CenteredNorm(),cmap='RdBu_r') # For old version of Fig 8
+ c1 = ax[1].pcolormesh(times/3600.,zobs,C,norm=colors.CenteredNorm(),cmap='RdBu_r')
ax[1].contour(times/3600.,zobs,theta.mean(axis=2).T,
contours,
colors='black',
@@ -484,13 +487,20 @@ def plot_free_ensemble(cbl,show='spread',ax=None,cmap='viridis'):
elif show=='mean':
bkgd = theta.mean(axis=0)
contour_color = 'white'
- if show=='cov':
+ elif show=='cov':
cov = np.zeros((nz,ntimes))+np.nan
for ii in range(nz):
for jj in range(ntimes):
cov[ii,jj] = np.cov(param,theta[:,ii,jj])[0,1]
bkgd = cov
contour_color = 'black'
+ elif show=='corr':
+ corr = np.zeros((nz,ntimes))+np.nan
+ for ii in range(nz):
+ for jj in range(ntimes):
+ corr[ii,jj] = np.corrcoef(param,theta[:,ii,jj])[0,1]
+ bkgd = corr
+ contour_color = 'black'
# Plot
zmax = 2000
@@ -498,7 +508,7 @@ def plot_free_ensemble(cbl,show='spread',ax=None,cmap='viridis'):
if ax is not None:
# Make plots
- c = ax.pcolormesh(times/3600.,zt,bkgd,cmap=cmap)
+ c = ax.pcolormesh(times/3600.,zt,bkgd,cmap=cmap,norm=colors.CenteredNorm())
ax.contour(times/3600.,zt,theta.mean(axis=0),
np.linspace(cbl.theta_0,cbl.theta_0+cbl.gamma*zmax,ncont),
colors=contour_color,
@@ -507,7 +517,7 @@ def plot_free_ensemble(cbl,show='spread',ax=None,cmap='viridis'):
ax.set_ylim([0,zmax])
return ax,c
-def plot_diagnostics(exp,show='errors',axes=None,label='',linecolor='blue'):
+def plot_diagnostics(exp,show='errors',axes=None,label='',linecolor='blue',zmax=None):
# See if the experiment is an OSSE
if (not hasattr(exp.da,'truths')):
@@ -517,9 +527,11 @@ def plot_diagnostics(exp,show='errors',axes=None,label='',linecolor='blue'):
# Get the data
sigma_o = exp.dg.varo**0.5
+ sigma_d = exp.dg.sigma_ominb
armsd = exp.dg.aRMSD_t
brmsd = exp.dg.bRMSD_t
- bsprd = exp.dg.bSprd_t
+ #bsprd = exp.dg.bSprd_t
+ bsprd = exp.dg.varb_ta**0.5
if truth_exists:
armse = exp.dg.aRMSE_t
brmse = exp.dg.bRMSE_t
@@ -537,15 +549,20 @@ def plot_diagnostics(exp,show='errors',axes=None,label='',linecolor='blue'):
elif exp.da.obs_coordinates.ndim == 2:
z = np.sort(exp.da.obs_coordinates[0,:])
+ if zmax is None:
+ idx = np.arange(z.size)
+ else:
+ idx = (z<zmax)
+
if axes is not None and len(axes)==4:
# Plot errors (wrt to true state) only if true state exists
if show=='errors' and truth_exists:
- axes[0].plot(armse,z,label=label,color=linecolor)#,color=linecolors[i])
- axes[1].plot(brmse,z,color=linecolor)
- axes[2].plot(brmse-armse,z,color=linecolor)
- axes[3].plot(np.sqrt(bsprd**2+sigma_o**2)/brmse,z,color=linecolor)
+ axes[0].plot(armse[idx],z[idx],label=label,color=linecolor)#,color=linecolors[i])
+ axes[1].plot(brmse[idx],z[idx],color=linecolor)
+ axes[2].plot(brmse[idx]-armse[idx],z[idx],color=linecolor)
+ axes[3].plot(np.sqrt(bsprd[idx]**2+sigma_o[idx]**2)/brmse[idx],z[idx],color=linecolor)
elif show=='deviations':
@@ -553,10 +570,10 @@ def plot_diagnostics(exp,show='errors',axes=None,label='',linecolor='blue'):
mean_bias_a_i1 = exp.dg.omina.mean(axis=0)
mean_bias_b_i1 = exp.dg.ominb.mean(axis=0)
# Plots
- axes[0].plot(armsd,z,label=label,color=linecolor)
- axes[1].plot(brmsd,z,color=linecolor)
- axes[2].plot(brmsd-armsd,z,color=linecolor)
- axes[3].plot(np.sqrt(bsprd**2+sigma_o**2)/brmsd,z,color=linecolor)
+ axes[0].plot(armsd[idx],z[idx],label=label,color=linecolor)
+ axes[1].plot(brmsd[idx],z[idx],color=linecolor)
+ axes[2].plot(brmsd[idx]-armsd[idx],z[idx],color=linecolor)
+ axes[3].plot(np.sqrt(bsprd[idx]**2+sigma_o[idx]**2)/brmsd[idx],z,color=linecolor)
return axes
@@ -593,7 +610,7 @@ def plot_linearity(cbl,ax=None):
return ax
-def plot_B_matrix(exp,timeid=-1):
+def plot_B_matrix(exp,timeid=-1,localize=True):
# Read the background state for the requested analysis time
xx = exp.da.backgrounds[timeid]
@@ -612,17 +629,30 @@ def plot_B_matrix(exp,timeid=-1):
Pb = np.matmul(X.T,X)/(nens-1)
v = np.abs(Pb[:-1,:-1]).max()
+ # Compute localized covariance matrix
+ xcoords = exp.nr.zt
+ dist = compute_distances(xcoords,xcoords)
+ loc = np.ones((nvar,nvar))
+ loc[:(nvar-1),:(nvar-1)] = gaspari_cohn(dist,cutoff=800)
+
# Plot
- fig = p.figure(0)
- fig.set_size_inches(4,3)
- sp1 = fig.add_subplot(1,1,1)
- c = sp1.pcolormesh(range(nvar),range(nvar),Pb,vmin=-v,vmax=v,cmap='RdBu_r',shading='nearest')
+ fig, [sp1,sp2] = p.subplots(1,2,constrained_layout=True)
+ fig.set_size_inches(6,3)
+ #
+ c1 = sp1.pcolormesh(range(nvar),range(nvar),Pb,vmin=-v,vmax=v,cmap='RdBu_r',shading='nearest')
sp1.set_aspect('equal','box')
sp1.set_xlim(-0.5,nvar-0.5)
sp1.set_ylim(-0.5,nvar-0.5)
- sp1.set_title('$1/(N-1)\cdot X^TX$, $t=%4.2f$ h'%(tt/3600))
- p.colorbar(c)
- fig.tight_layout()
+ sp1.set_title('a) $P^b$, $t=%4.2f$ h'%(tt/3600))
+ #
+ c2 = sp2.pcolormesh(range(nvar),range(nvar),loc*Pb,vmin=-v,vmax=v,cmap='RdBu_r',shading='nearest')
+ sp2.set_aspect('equal','box')
+ sp2.set_xlim(-0.5,nvar-0.5)
+ sp2.set_ylim(-0.5,nvar-0.5)
+ sp2.set_title('b) Localized $P^b$, $t=%4.2f$ h'%(tt/3600))
+ #
+ p.colorbar(c1,orientation='horizontal')
+ p.colorbar(c2,orientation='horizontal')
fig.savefig('Covariance_matrix_%05u_s.png'%tt,format='png',dpi=150)
p.close(fig)
@@ -636,18 +666,21 @@ def plot_desroziers(dgs,filename='desroziers_diagnostics.png',labels=None,logsca
fig.set_size_inches(4,3)
ax = fig.add_subplot(1,1,1)
- for i,dg in zip( range(len(dgs)) , dgs ):
+ colors = p.rcParams['axes.prop_cycle'].by_key()['color'][:len(dgs)]
+ for i,dg,color in zip( range(len(dgs)) , dgs , colors):
if labels is not None:
# Two formulations of the same thing; they give the same numbers!
#ax1.scatter( dg.dobdob_ta-dg.mean_ominb**2 , dg.varo+dg.varb_ta , s=5, label = labels[i])
- ax.scatter( dg.sigma_ominb**2 , dg.varo+dg.varb_ta , s=5, label = labels[i])
+ ax.scatter( np.mean(dg.sigma_ominb**2) , np.mean(dg.varo+dg.varb_ta) , s=13, label = labels[i],zorder=10, c=color, edgecolors='black')
+ ax.scatter( dg.sigma_ominb**2 , dg.varo+dg.varb_ta , s=5, alpha = 0.3, c=color, linewidths=0)
else:
# Two formulations of the same thing; they give the same numbers!
#ax1.scatter( dg.dobdob_ta-dg.mean_ominb**2 , dg.varo+dg.varb_ta , s=5)
- ax.scatter( dg.sigma_ominb**2 , dg.varo+dg.varb_ta , s=5) #
+ ax.scatter( np.mean(dg.sigma_ominb**2) , np.mean(dg.varo+dg.varb_ta) , s=13, zorder = 10)
+ ax.scatter( dg.sigma_ominb**2 , dg.varo+dg.varb_ta , s=5, alpha = 0.3) #
#ax.set_title("Consistency of innovations",fontsize=8)
- ax.set_xlabel(r'${\langle}\sigma_d^2{\rangle}_t$') # '${\langle}d_{ob}^2{\rangle}_t$'
- ax.set_ylabel(r'$\sigma_o^2+{\langle}\sigma_b^2{\rangle}_t$')
+ ax.set_xlabel(r'${\langle}\sigma_d^2{\rangle}$ (K)') # '${\langle}d_{ob}^2{\rangle}_t$'
+ ax.set_ylabel(r'$\sigma_{y_o}^2+{\langle}\sigma_{y_b}^2{\rangle}$ (K)')
if labels is not None:
ax.legend(frameon=False, fontsize=6, ncol=2)
@@ -685,21 +718,24 @@ def plot_overview(exp,label,filename,show_spread=False,show_transformed_paramete
ax2.set_title(r'c) Exp. %s, $\delta\overline{p}\prime\prime$'%label)
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
- ax3 = plot_CBL_PE(exp,None,ax=ax3,show_spread=show_spread,show_transformed_parameter=show_transformed_parameter,show_true_parameter=show_true_parameter)
- if show_transformed_parameter:
- ax3.set_title(r'd) Exp. %s, evolution of $p\prime\prime$'%label)
- else:
- ax3.set_title(r'd) Exp. %s, evolution of $p$'%label)
- ax3.set_yticks(np.arange(5+1))
- ax3.set_xlabel('Time (h)')
ax0.set_xlim([0,maxtime/3600])
ax1.set_xlim([0,maxtime/3600])
ax2.set_xlim([0,maxtime/3600])
- ax3.set_xlim([0,maxtime/3600])
ax0.set_xticks(np.arange((maxtime+1)/3600))
ax1.set_xticks(np.arange((maxtime+1)/3600))
ax2.set_xticks(np.arange((maxtime+1)/3600))
- ax3.set_xticks(np.arange((maxtime+1)/3600))
+ ax3 = plot_CBL_PE(exp,None,ax=ax3,show_spread=show_spread,show_transformed_parameter=show_transformed_parameter,show_true_parameter=show_true_parameter)
+ if ax3 is None:
+ pass
+ else:
+ if show_transformed_parameter:
+ ax3.set_title(r'd) Exp. %s, evolution of $p\prime\prime$'%label)
+ else:
+ ax3.set_title(r'd) Exp. %s, evolution of $p$'%label)
+ ax3.set_yticks(np.arange(5+1))
+ ax3.set_xlabel('Time (h)')
+ ax3.set_xlim([0,maxtime/3600])
+ ax3.set_xticks(np.arange((maxtime+1)/3600))
p.colorbar(c0,orientation='horizontal')
p.colorbar(c1,orientation='horizontal')
p.colorbar(c2,orientation='horizontal')