error bars in parameter estimation convergence plots are now based on percentiles

graphics.py

@@ -9,17 +9,15 @@ def plot_CBL_assimilation(exp,figtitle,which_cycle=1):
     ensmembers_color = 'lightskyblue'
     ensmean_color = 'royalblue'
     
-    def ensplot(ax,x,nr,x1,x2,z):
+    def ensplot(ax,x,x1,x2,z):
         nens = x.shape[1]
-        ax.plot(nr[x1:x2],zt,color=naturerun_color,zorder=2)
         ax.plot(x[x1:x2].mean(axis=1),z,color=ensmean_color,zorder=5)
         for i in range(nens): ax.plot(x[x1:x2,i],z,color=ensmembers_color,alpha=0.2,zorder=1)
 
     # Shorthand
-    nr = exp.nr
-    o = exp.observations[which_cycle-1]
-    b = exp.da.backgrounds[which_cycle-1,:,:]
-    a = exp.da.analyses[which_cycle-1,:,:]
+    o = exp.observations[which_cycle]
+    b = exp.da.backgrounds[which_cycle,:,:]
+    a = exp.da.analyses[which_cycle,:,:]
     ocoords = exp.obs_coordinates
     nz = exp.nr.nz
     zt = exp.nr.zt 
@@ -36,19 +34,18 @@ def plot_CBL_assimilation(exp,figtitle,which_cycle=1):
     # Potential temperature
     #
     ax1 = fig.add_subplot(ii,4,1)
-    ax1.plot(nr.x[:nz],zt,color=naturerun_color)
     ax1.scatter(o,ocoords,marker="o",s=20,color=naturerun_color,zorder=10)
     ax1.set_title(r'$\theta$ (K), nature run')
     ax1.set_ylabel(r'$z$ (m)')
     ax1.set_ylim([0,zmax])
     #
     ax2 = fig.add_subplot(ii,4,2,sharey=ax1)
-    ensplot(ax2,b,nr.x,0,nz,zt)
+    ensplot(ax2,b,0,nz,zt)
     ax2.scatter(o,ocoords,marker="o",s=20,color=naturerun_color,zorder=10)
     ax2.set_title(r'$\theta$ (K), prior')
     #
     ax3 = fig.add_subplot(ii,4,3,sharey=ax1)
-    ensplot(ax3,a,nr.x,0,nz,zt)
+    ensplot(ax3,a,0,nz,zt)
     ax3.scatter(o,ocoords,marker="o",s=20,color=naturerun_color,zorder=10)
     ax3.set_title(r'$\theta$ (K), posterior')
     #
@@ -66,17 +63,16 @@ def plot_CBL_assimilation(exp,figtitle,which_cycle=1):
         # u wind component
         #
         ax5 = fig.add_subplot(ii,4,5)
-        ax5.plot(nr.x[nz:nz*2],zt,color=naturerun_color)
         ax5.set_title(r'$u$ (m/s), nature run')
         ax5.set_ylabel(r'$z$ (m)')
         ax5.set_ylim([0,zmax])
         #
         ax6 = fig.add_subplot(ii,4,6,sharey=ax1)
-        ensplot(ax6,b,nr.x,nz,nz*2,zt)
+        ensplot(ax6,b,nz,nz*2,zt)
         ax6.set_title(r'$u$ (m/s), prior')
         #
         ax7 = fig.add_subplot(ii,4,7,sharey=ax1)
-        ensplot(ax7,a,nr.x,nz,nz*2,zt)
+        ensplot(ax7,a,nz,nz*2,zt)
         ax7.set_title(r'$u$ (m/s), posterior')
         #
         ax8 = fig.add_subplot(ii,4,8,sharey=ax1)
@@ -86,17 +82,16 @@ def plot_CBL_assimilation(exp,figtitle,which_cycle=1):
         # v wind component
         #
         ax9 = fig.add_subplot(ii,4,9)
-        ax9.plot(nr.x[nz*2:nz*3],zt,color=naturerun_color)
         ax9.set_title(r'$v$ (m/s), nature run')
         ax9.set_ylabel(r'$z$ (m)')
         ax9.set_ylim([0,zmax])
         #
         axa = fig.add_subplot(ii,4,10,sharey=ax1)
-        ensplot(axa,b,nr.x,nz*2,nz*3,zt)
+        ensplot(axa,b,nz*2,nz*3,zt)
         axa.set_title(r'$v$ (m/s), prior')
         #
         axb = fig.add_subplot(ii,4,11,sharey=ax1)
-        ensplot(axb,a,nr.x,nz*2,nz*3,zt)
+        ensplot(axb,a,nz*2,nz*3,zt)
         axb.set_title(r'$v$ (m/s), posterior')
         #
         axc = fig.add_subplot(ii,4,12,sharey=ax1)
@@ -126,7 +121,7 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False):
         parameter_number = exp.nr.parameter_number
         b = exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:]
         par_b = pt[parameter_id](exp.da.backgrounds[:,-parameter_number:,:][:,parameter_id,:],kind='inv')
-        t = exp.da.time
+        t = exp.da.time/3600.
 
         if plot_spread:
             ii=2
@@ -139,8 +134,8 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False):
         ax1 = fig.add_subplot(ii,1,1)
         ax1.step(t,par_b.mean(axis=1),color=mean_color)
         for i in range(t.size-1):
-            bmin=(par_b.mean(axis=1)-par_b.std(axis=1))[i+1]
-            bmax=(par_b.mean(axis=1)+par_b.std(axis=1))[i+1]
+            bmin = np.percentile(par_b, 10, axis=1)[i+1]
+            bmax = np.percentile(par_b, 90, axis=1)[i+1]
             ax1.fill_between(t[i:i+2],\
                         y1=[bmin,bmin],\
                         y2=[bmax,bmax],\
@@ -149,7 +144,8 @@ def plot_CBL_PE(exp,figtitle,parameter_id=0,plot_spread=False):
                         alpha=0.3)
         ax1.axhline(y=true_value,linestyle='--',color='black')
         ax1.set_title('Parameter evolution')
-        ax1.set_xlabel('Time (s)')
+        ax1.set_xlabel('Time (h)')
+        ax1.set_xlim([0,t.max()])
         #
         if plot_spread:
             ax2 = fig.add_subplot(ii,1,2)