[Dart-dev] DART/branches Revision: 12976

[dart at ucar.edu]

thoar at ucar.edu
2019-02-08 16:51:18 -0700 (Fri, 08 Feb 2019)
263
First attempt at fully supporting obs_diag output with no posterior information.
Attempts to plot number of observations used for the prior _and_ the posterior, separately.
Checks to see if the observation was 'assimilated' or 'evaluated' and annotates as such.




Modified: DART/branches/recam/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/branches/recam/diagnostics/matlab/plot_rmse_xxx_evolution.m	2019-02-08 19:52:39 UTC (rev 12975)
+++ DART/branches/recam/diagnostics/matlab/plot_rmse_xxx_evolution.m	2019-02-08 23:51:18 UTC (rev 12976)
@@ -124,11 +124,11 @@
 plotdat.copystring    = copy;
 plotdat.bincenters    = ncread(fname,'time');
 plotdat.binedges      = ncread(fname,'time_bounds');
-plotdat.mlevel        = local_nc_varget(fname,'mlevel');
-plotdat.plevel        = local_nc_varget(fname,'plevel');
-plotdat.plevel_edges  = local_nc_varget(fname,'plevel_edges');
-plotdat.hlevel        = local_nc_varget(fname,'hlevel');
-plotdat.hlevel_edges  = local_nc_varget(fname,'hlevel_edges');
+plotdat.mlevel        = local_ncread(fname,'mlevel');
+plotdat.plevel        = local_ncread(fname,'plevel');
+plotdat.plevel_edges  = local_ncread(fname,'plevel_edges');
+plotdat.hlevel        = local_ncread(fname,'hlevel');
+plotdat.hlevel_edges  = local_ncread(fname,'hlevel_edges');
 [plotdat.ncopies,~]   = nc_dim_info(fname,'copy');
 [plotdat.nregions,~]  = nc_dim_info(fname,'region');
 plotdat.region_names  = strtrim(ncread(fname,'region_names')');
@@ -178,6 +178,10 @@
 plotdat.rmseindex   = get_copy_index(fname,'rmse');
 plotdat.Npossindex  = get_copy_index(fname,'Nposs');
 plotdat.Nusedindex  = get_copy_index(fname,'Nused');
+plotdat.NQC0index   = get_copy_index(fname,'N_DARTqc_0');
+plotdat.NQC1index   = get_copy_index(fname,'N_DARTqc_1');
+plotdat.NQC2index   = get_copy_index(fname,'N_DARTqc_2');
+plotdat.NQC3index   = get_copy_index(fname,'N_DARTqc_3');
 plotdat.NQC4index   = get_copy_index(fname,'N_DARTqc_4');
 plotdat.NQC5index   = get_copy_index(fname,'N_DARTqc_5');
 plotdat.NQC6index   = get_copy_index(fname,'N_DARTqc_6');
@@ -185,6 +189,14 @@
 
 figuredata = setfigure();
 
+prior_green = [  0/255 128/255   0/255];
+poste_blue  = [  0/255   0/255 255/255];
+obs_red     = [215/255  10/255  83/255];
+
+plotdat.ges_color = prior_green;
+plotdat.anl_color = poste_blue;
+plotdat.obs_color = obs_red;
+
 %%---------------------------------------------------------------------
 % Loop around (time-copy-level-region) observation types
 %----------------------------------------------------------------------
@@ -246,10 +258,17 @@
     guess = reshape(guess_raw, plotdat.Nbins,   plotdat.ncopies, ...
         plotdat.nlevels, plotdat.nregions);
     
-    analy_raw = ncread(fname, plotdat.analyvar);
-    analy_raw = permute(analy_raw,length(size(analy_raw)):-1:1);
-    analy = reshape(analy_raw, plotdat.Nbins,   plotdat.ncopies, ...
-        plotdat.nlevels, plotdat.nregions);
+    analy_raw = local_ncread(fname, plotdat.analyvar);
+    if ( isempty(analy_raw) )
+        % force analysis to be the same shape as the guess,
+        % and full of NaNs instead of 0s.
+        analy = guess;
+        analy(:) = NaN;
+    else
+        analy_raw = permute(analy_raw,length(size(analy_raw)):-1:1);
+        analy = reshape(analy_raw, plotdat.Nbins,   plotdat.ncopies, ...
+            plotdat.nlevels, plotdat.nregions);
+    end
     
     % check to see if there is anything to plot
     % The number possible is decreased by the number of observations
@@ -281,6 +300,26 @@
     for ilevel = wantedlevels
         
         fprintf(logfid,'\nlevel %d %f %s\n',ilevel,plotdat.level(ilevel),plotdat.level_units);
+        plotdat.ges_Nqc0  = guess(:,plotdat.NQC0index  ,ilevel,:);
+        plotdat.anl_Nqc0  = analy(:,plotdat.NQC0index  ,ilevel,:);
+        fprintf(logfid,'DART QC == 0, prior/post %d %d\n',sum(plotdat.ges_Nqc0(:)), ...
+            sum(plotdat.anl_Nqc0(:)));
+
+        plotdat.ges_Nqc1  = guess(:,plotdat.NQC1index  ,ilevel,:);
+        plotdat.anl_Nqc1  = analy(:,plotdat.NQC1index  ,ilevel,:);
+        fprintf(logfid,'DART QC == 1, prior/post %d %d\n',sum(plotdat.ges_Nqc1(:)), ...
+            sum(plotdat.anl_Nqc1(:)));
+
+        plotdat.ges_Nqc2  = guess(:,plotdat.NQC2index  ,ilevel,:);
+        plotdat.anl_Nqc2  = analy(:,plotdat.NQC2index  ,ilevel,:);
+        fprintf(logfid,'DART QC == 2, prior/post %d %d\n',sum(plotdat.ges_Nqc2(:)), ...
+            sum(plotdat.anl_Nqc2(:)));
+
+        plotdat.ges_Nqc3  = guess(:,plotdat.NQC3index  ,ilevel,:);
+        plotdat.anl_Nqc3  = analy(:,plotdat.NQC3index  ,ilevel,:);
+        fprintf(logfid,'DART QC == 3, prior/post %d %d\n',sum(plotdat.ges_Nqc3(:)), ...
+            sum(plotdat.anl_Nqc3(:)));
+
         plotdat.ges_Nqc4  = guess(:,plotdat.NQC4index  ,ilevel,:);
         plotdat.anl_Nqc4  = analy(:,plotdat.NQC4index  ,ilevel,:);
         fprintf(logfid,'DART QC == 4, prior/post %d %d\n',sum(plotdat.ges_Nqc4(:)), ...
@@ -360,47 +399,51 @@
 
 function myplot(plotdat,figdata)
 
-% Interlace the [ges,anl] to make a sawtooth plot.