[Dart-dev] [4322] DART/trunk/diagnostics/matlab: Changed the legend format. Moved the

nancy at ucar.edu nancy at ucar.edu
Thu Mar 18 13:14:12 MDT 2010


Revision: 4322
Author:   thoar
Date:     2010-03-18 13:14:11 -0600 (Thu, 18 Mar 2010)
Log Message:
-----------
Changed the legend format. Moved the summary statistics to be part of the
title - removed some of the axis labelling clutter. Also prints the full path
of the input file at the bottom of the page. Added feature to be able to
plot just one observation type instead of ALL the observation types in the file.
Also prints out a separate file for each observation type with summary 
observation counts by QC flag, etc - for each level. 

Modified Paths:
--------------
    DART/trunk/diagnostics/matlab/plot_evolution.m
    DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m

-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_evolution.m	2010-03-18 15:50:22 UTC (rev 4321)
+++ DART/trunk/diagnostics/matlab/plot_evolution.m	2010-03-18 19:14:11 UTC (rev 4322)
@@ -1,22 +1,46 @@
-function plotdat = plot_evolution(fname,copystring)
+function plotdat = plot_evolution(fname, copystring, varargin)
 %% plot_evolution plots the temporal evolution of the observation-space quantities for all possible levels, all possible variables.
 % Part of the observation-space diagnostics routines.
 %
 % 'obs_diag' produces a netcdf file containing the diagnostics.
+% obs_diag condenses the obs_seq.final information into summaries for a few specified
+% regions - on a level-by-level basis.
 %
-% USAGE: plotdat = plot_evolution(fname,copystring);
+% USAGE: plotdat = plot_evolution(fname, copystring [,varargin]);
 %
-% fname  :  netcdf file produced by 'obs_diag'
-% copystring :  'copy' string == quantity of interest. These
-%            can be any of the ones available in the netcdf 
-%            file 'CopyMetaData' variable.
-%            (ncdump -v CopyMetaData obs_diag_output.nc)
+% fname         : netcdf file produced by 'obs_diag'
+% copystring    : 'copy' string == quantity of interest. These
+%                 can be any of the ones available in the netcdf 
+%                 file 'CopyMetaData' variable.
+%                 (ncdump -v CopyMetaData obs_diag_output.nc)
+% obstypestring : 'observation type' string. Optional.
+%                 Must match something in the netcdf 
+%                 file 'ObservationTypes' variable.
+%                 (ncdump -v ObservationTypes obs_diag_output.nc)
+%                 If specified, only this observation type will be plotted.
+%                 If not specified, all observation types incluced in the netCDF file
+%                 will be plotted.
 %
-% EXAMPLE:
+% OUTPUT: two files will result for each observation type plotted. One is a 
+%         postscript file containing a page for each level - all regions.
+%         The other file is a simple text file containing summary information
+%         about how many observations were assimilated, how many were available, etc.
+%         Both of these filenames contain the observation type as part of the name.
 %
-% fname = 'obs_diag_output.nc';   % netcdf file produced by 'obs_diag'
-% copystring = 'bias';   % 'copy' string == quantity of interest
-% plotdat = plot_evolution(fname,copystring);
+% EXAMPLE 1 - plot the evolution of the bias for all observation types, all levels
+%
+% fname      = 'obs_diag_output.nc';   % netcdf file produced by 'obs_diag'
+% copystring = 'bias';                 % 'copy' string == quantity of interest
+% plotdat    = plot_evolution(fname, copystring);
+%
+%
+% EXAMPLE 2 - plot the evolution of the rmse for just the radiosonde temperature observations
+%             This requires that the 'RADIOSONDE_TEMPERATURE' is one of the known observation
+%             types in the netCDF file.
+%
+% fname      = 'obs_diag_output.nc';   % netcdf file produced by 'obs_diag'
+% copystring = 'rmse';                 % 'copy' string == quantity of interest
+% plotdat    = plot_evolution(fname, copystring, 'RADIOSONDE_TEMPERATURE');
 
 %% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
 % provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -28,6 +52,16 @@
 % $Revision$
 % $Date$
 
+if nargin == 2
+   nvars = 0;
+elseif nargin == 3
+   varname = varargin{1};
+   nvars = 1;
+else
+   error('wrong number of arguments ... ')
+end
+
+
 if (exist(fname,'file') ~= 2)
    error('file/fname <%s> does not exist',fname)
 end
@@ -47,7 +81,7 @@
 plotdat.nregions      = length(nc_varget(fname,'region'));
 plotdat.region_names  = nc_varget(fname,'region_names');
 
-if (plotdat.nregions == 1)
+if (plotdat.nregions == 1 && (size(plotdat.region_names,2) == 1) )
    plotdat.region_names = deblank(plotdat.region_names');
 end
 
@@ -80,13 +114,23 @@
 
 plotdat.linewidth = 2.0;
 
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
-[plotdat.varnames,    plotdat.vardims]    = FindTemporalVars(plotdat);
-plotdat.nvars       = length(plotdat.varnames);
+if (nvars == 0)
+   [plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
+   [plotdat.varnames,    plotdat.vardims]    = FindTemporalVars(plotdat);
+   plotdat.nvars       = length(plotdat.varnames);
+else
+   plotdat.varnames{1} = varname;
+   plotdat.nvars       = nvars;
+end
 
+
 plotdat.copyindex   = get_copy_index(fname,copystring); 
 plotdat.Npossindex  = get_copy_index(fname,'Nposs');
 plotdat.Nusedindex  = get_copy_index(fname,'Nused');
+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');
+plotdat.NQC7index   = get_copy_index(fname,'N_DARTqc_7');
 
 %----------------------------------------------------------------------
 % Loop around (time-copy-level-region) observation types
@@ -107,6 +151,14 @@
    disp(sprintf('Removing %s from the current directory.',psfname))
    system(sprintf('rm %s',psfname));
 
+   % remove any existing log file - 
+   
+   lgfname = sprintf('%s_obscount.txt',plotdat.varnames{ivar});
+   disp(sprintf('Removing %s from the current directory.',lgfname))
+   system(sprintf('rm %s',lgfname));
+   logfid = fopen(lgfname,'wt');
+   fprintf(logfid,'%s\n',lgfname);
+
    % get appropriate vertical coordinate variable
 
    guessdims = nc_var_dims(fname, plotdat.guessvar);
@@ -145,15 +197,42 @@
    
    for ilevel = 1:plotdat.nlevels
 
-      plotdat.ges_copy   = guess(:,plotdat.copyindex,  ilevel,:);
-      plotdat.anl_copy   = analy(:,plotdat.copyindex,  ilevel,:);
+      fprintf(logfid,'\nlevel %d %f %s\n',ilevel,plotdat.level(ilevel),plotdat.level_units);
+      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(:)), ...
+                                                 sum(plotdat.anl_Nqc4(:)));
 
-      plotdat.ges_Nposs  = guess(:,plotdat.Npossindex, ilevel,:);
-      plotdat.anl_Nposs  = analy(:,plotdat.Npossindex, ilevel,:);
-      plotdat.ges_Nused  = guess(:,plotdat.Nusedindex, ilevel,:);
-      plotdat.anl_Nused  = analy(:,plotdat.Nusedindex, ilevel,:);
-      plotdat.Yrange     = FindRange(plotdat);
-      
+      plotdat.ges_Nqc5  = guess(:,plotdat.NQC5index  ,ilevel,:);
+      plotdat.anl_Nqc5  = analy(:,plotdat.NQC5index  ,ilevel,:);
+      fprintf(logfid,'DART QC == 5, prior/post %d %d\n',sum(plotdat.ges_Nqc5(:)), ...
+                                                 sum(plotdat.anl_Nqc5(:)));
+
+      plotdat.ges_Nqc6  = guess(:,plotdat.NQC6index  ,ilevel,:);
+      plotdat.anl_Nqc6  = analy(:,plotdat.NQC6index  ,ilevel,:);
+      fprintf(logfid,'DART QC == 6, prior/post %d %d\n',sum(plotdat.ges_Nqc6(:)), ...
+                                                 sum(plotdat.anl_Nqc6(:)));
+
+      plotdat.ges_Nqc7  = guess(:,plotdat.NQC7index  ,ilevel,:);
+      plotdat.anl_Nqc7  = analy(:,plotdat.NQC7index  ,ilevel,:);
+      fprintf(logfid,'DART QC == 7, prior/post %d %d\n',sum(plotdat.ges_Nqc7(:)), ...
+                                                 sum(plotdat.anl_Nqc7(:)));
+
+      plotdat.ges_Nposs = guess(:,plotdat.Npossindex, ilevel,:);
+      plotdat.anl_Nposs = analy(:,plotdat.Npossindex, ilevel,:);
+      fprintf(logfid,'# obs poss,   prior/post %d %d\n',sum(plotdat.ges_Nposs(:)), ...
+                                                        sum(plotdat.anl_Nposs(:)));
+
+      plotdat.ges_Nused = guess(:,plotdat.Nusedindex, ilevel,:);
+      plotdat.anl_Nused = analy(:,plotdat.Nusedindex, ilevel,:);
+      fprintf(logfid,'# obs used,   prior/post %d %d\n',sum(plotdat.ges_Nused(:)), ...
+                                                        sum(plotdat.anl_Nused(:)));
+
+      plotdat.ges_copy  = guess(:,plotdat.copyindex,  ilevel,:);
+      plotdat.anl_copy  = analy(:,plotdat.copyindex,  ilevel,:);
+
+      plotdat.Yrange    = FindRange(plotdat);
+
       % plot by region
 
       if (plotdat.nregions > 2)
@@ -219,6 +298,7 @@
 
    string_guess = sprintf('forecast: mean=%.5g', mean_prior);
    string_analy = sprintf('analysis: mean=%.5g', mean_post);
+   plotdat.subtitle = sprintf('%s     %s     %s',plotdat.myregion, string_guess, string_analy);
 
    % Plot the requested quantity on the left axis.
    % The observation count will use the axis on the right.
@@ -228,22 +308,23 @@
    
    ax1 = subplot(plotdat.nregions,1,plotdat.region);
 
-   h1 = plot(ta,ca,'ro-',tg,cg,'k+-','LineWidth',plotdat.linewidth);
-   h = legend(string_analy, string_guess);
+   h1 = plot(tg,cg,'k+-',ta,ca,'ro-','LineWidth',plotdat.linewidth);
+   h = legend('forecast', 'analysis');
    legend(h,'boxoff')
 
    axlims = axis;
    axlims = [axlims(1:2) plotdat.Yrange];
    axis(axlims)
 
+   plotdat.ylabel{1} = plotdat.myregion;
    switch lower(plotdat.copystring)
       case 'bias'
          % plot a zero-bias line
          h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
          set(h4,'LineWidth',1.5,'LineSTyle',':')
-         plotdat.ylabel = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
+         plotdat.ylabel{2} = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
       otherwise
-         plotdat.ylabel = sprintf('%s',plotdat.copystring);
+         plotdat.ylabel{2} = sprintf('%s',plotdat.copystring);
    end
    
    % hokey effort to decide to plot months/days vs. daynum vs.
@@ -252,14 +333,30 @@
    if ((plotdat.bincenters(1) > 1000) && (ttot > 5))
       datetick('x',6,'keeplimits','keepticks');
       monstr = datestr(plotdat.bincenters(1),21);
-      xlabel(sprintf('month/day - %s start',monstr))
+      xlabelstring = sprintf('month/day - %s start',monstr);
    elseif (plotdat.bincenters(1) > 1000)
       datetick('x',15,'keeplimits','keepticks')
       monstr = datestr(plotdat.bincenters(1),21);
-      xlabel(sprintf('%s start',monstr))
+      xlabelstring = sprintf('%s start',monstr);
    else
-      xlabel('days')
+      xlabelstring = 'days';
    end
+
+   % only put x axis on last/bottom plot
+   if (plotdat.region == plotdat.nregions)
+      xlabel(xlabelstring)
+   end
+
+   % more annotation ...
+
+   if (plotdat.region == 1)
+      title({plotdat.title, plotdat.subtitle}, ...
+         'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
+      BottomAnnotation(plotdat.fname)
+   else
+      title(plotdat.subtitle, 'Interpreter', 'none', ...
+         'Fontsize', 12, 'FontWeight', 'bold')
+   end
    
    % use same X,Y limits for all plots in this region
    nYticks = length(get(ax1,'YTick'));
@@ -285,7 +382,7 @@
    yinc   = (ylimits(2)-ylimits(1))/(nYticks-1);
    yticks = ylimits(1):yinc:ylimits(2);
    niceyticks = round(10*yticks')/10;
-   set(ax2,'XTick',get(ax1,'XTick'),'XTicklabel',get(ax1,'XTicklabel'), ...
+   set(ax2,'XTick',get(ax1,'XTick'),'XTicklabel',[], ...
            'YTick',          yticks,'YTicklabel',num2str(niceyticks))
        
    set(get(ax2,'Ylabel'),'String','# of obs : o=poss, +=used')
@@ -293,13 +390,6 @@
    set(ax1,'Position',get(ax2,'Position'))
    grid
 
-   if (plotdat.region == 1)
-   title({plotdat.title, plotdat.myregion}, ...
-         'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
-   else
-   title(plotdat.myregion, 'Interpreter', 'none', ...
-         'Fontsize', 12, 'FontWeight', 'bold')
-   end
 
 
 
@@ -307,11 +397,18 @@
 % annotates the directory containing the data being plotted
 subplot('position',[0.48 0.01 0.04 0.04])
 axis off
-h = text(0.0,0.5,main);
+str = which(main);
+if ( isempty(str) )   % must be a relative filename 
+   mydir = pwd;
+   string1 = sprintf('data file: %s/%s',mydir,main);
+else
+   string1 = sprintf('data file: %s',str);
+end
+h = text(0.0, 0.5, string1);
 set(h,'HorizontalAlignment','center', ...
       'VerticalAlignment','middle',...
       'Interpreter','none',...
-      'FontSize',8)
+      'FontSize',10)
 
 
 

Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m	2010-03-18 15:50:22 UTC (rev 4321)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m	2010-03-18 19:14:11 UTC (rev 4322)
@@ -1,22 +1,46 @@
-function plotdat = plot_rmse_xxx_evolution(fname,copystring)
+function plotdat = plot_rmse_xxx_evolution(fname, copystring, varargin)
 %% plot_rmse_xxx_evolution plots the temporal evolution of the observation-space quantity RMSE and any other for all possible levels, all possible variables.
 % Part of the observation-space diagnostics routines.
 %
 % 'obs_diag' produces a netcdf file containing the diagnostics.
+% obs_diag condenses the obs_seq.final information into summaries for a few specified
+% regions - on a level-by-level basis.
 %
-% USAGE: plotdat = plot_rmse_xxx_evolution(fname,copystring);
+% USAGE: plotdat = plot_rmse_xxx_evolution(fname, copystring [,varargin]);
 %
-% fname  :  netcdf file produced by 'obs_diag'
-% copystring :  'copy' string == quantity of interest. These
-%            can be any of the ones available in the netcdf 
-%            file 'CopyMetaData' variable.
-%            (ncdump -v CopyMetaData obs_diag_output.nc)
+% fname         : netcdf file produced by 'obs_diag'
+% copystring    : 'copy' string == quantity of interest. These
+%                 can be any of the ones available in the netcdf 
+%                 file 'CopyMetaData' variable.
+%                 (ncdump -v CopyMetaData obs_diag_output.nc)
+% obstypestring : 'observation type' string. Optional.
+%                 Must match something in the netcdf 
+%                 file 'ObservationTypes' variable.
+%                 (ncdump -v ObservationTypes obs_diag_output.nc)
+%                 If specified, only this observation type will be plotted.
+%                 If not specified, all observation types incluced in the netCDF file
+%                 will be plotted.
 %
-% EXAMPLE: plot the RMSE and totalspread on the same axis.
+% OUTPUT: two files will result for each observation type plotted. One is a 
+%         postscript file containing a page for each level - all regions.
+%         The other file is a simple text file containing summary information
+%         about how many observations were assimilated, how many were available, etc.
+%         Both of these filenames contain the observation type as part of the name.
 %
-% fname = 'obs_diag_output.nc';   % netcdf file produced by 'obs_diag'
-% copystring = 'totalspread';   % 'copy' string == quantity of interest
-% plotdat = plot_rmse_xxx_evolution(fname,copystring);
+%
+% EXAMPLE 1 - plot the RMSE and totalspread on the same axis.
+%
+% fname      = 'obs_diag_output.nc';   % netcdf file produced by 'obs_diag'
+% copystring = 'totalspread';          % 'copy' string == quantity of interest
+% plotdat    = plot_rmse_xxx_evolution(fname,copystring);
+%
+%
+% EXAMPLE 2 - plot the RMSE and spread on the same axis - for just one observation type.
+%
+% fname      = 'obs_diag_output.nc';
+% copystring = 'totalspread';
+% obstype    = 'RADIOSONDE_TEMPERATURE';
+% plotdat    = plot_rmse_xxx_evolution(fname, copystring, obstype);
 
 %% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
 % provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -29,6 +53,15 @@
 % $Revision$
 % $Date$
 
+if nargin == 2
+   nvars = 0;
+elseif nargin == 3
+   varname = varargin{1};
+   nvars = 1;
+else
+   error('wrong number of arguments ... ')
+end
+
 if (exist(fname,'file') ~= 2)
    error('file/fname <%s> does not exist',fname)
 end
@@ -48,7 +81,7 @@
 plotdat.nregions      = length(nc_varget(fname,'region'));
 plotdat.region_names  = nc_varget(fname,'region_names');
 
-if (plotdat.nregions == 1)
+if (plotdat.nregions == 1 && (size(plotdat.region_names,2) == 1) )
    plotdat.region_names = deblank(plotdat.region_names');
 end
 
@@ -79,17 +112,26 @@
 
 % set up a structure with all static plotting components
 
-plotdat.ylabel    = sprintf('rmse and %s',copystring);
 plotdat.linewidth = 2.0;
 
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
-[plotdat.varnames,    plotdat.vardims]    = FindTemporalVars(plotdat);
-plotdat.nvars       = length(plotdat.varnames);
+if (nvars == 0)
+   [plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
+   [plotdat.varnames,    plotdat.vardims]    = FindTemporalVars(plotdat);
+   plotdat.nvars       = length(plotdat.varnames);
+else
+   plotdat.varnames{1} = varname;
+   plotdat.nvars       = nvars;
+end
 
+
 plotdat.copyindex   = get_copy_index(fname,copystring); 
 plotdat.rmseindex   = get_copy_index(fname,'rmse');
 plotdat.Npossindex  = get_copy_index(fname,'Nposs');
 plotdat.Nusedindex  = get_copy_index(fname,'Nused');
+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');
+plotdat.NQC7index   = get_copy_index(fname,'N_DARTqc_7');
 
 %----------------------------------------------------------------------
 % Loop around (time-copy-level-region) observation types
@@ -110,6 +152,14 @@
    disp(sprintf('Removing %s from the current directory.',psfname))
    system(sprintf('rm %s',psfname));
 
+   % remove any existing log file - 
+   
+   lgfname = sprintf('%s_obscount.txt',plotdat.varnames{ivar});
+   disp(sprintf('Removing %s from the current directory.',lgfname))
+   system(sprintf('rm %s',lgfname));
+   logfid = fopen(lgfname,'wt');
+   fprintf(logfid,'%s\n',lgfname);
+
    % get appropriate vertical coordinate variable
 
    guessdims = nc_var_dims(fname, plotdat.guessvar);
@@ -148,17 +198,44 @@
    
    for ilevel = 1:plotdat.nlevels
 
-      plotdat.ges_copy   = guess(:,plotdat.copyindex,  ilevel,:);
-      plotdat.anl_copy   = analy(:,plotdat.copyindex,  ilevel,:);
-      plotdat.ges_rmse   = guess(:,plotdat.rmseindex,  ilevel,:);
-      plotdat.anl_rmse   = analy(:,plotdat.rmseindex,  ilevel,:);
+      fprintf(logfid,'\nlevel %d %f %s\n',ilevel,plotdat.level(ilevel),plotdat.level_units);
+      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(:)), ...
+                                                 sum(plotdat.anl_Nqc4(:)));
 
-      plotdat.ges_Nposs  = guess(:,plotdat.Npossindex, ilevel,:);
-      plotdat.anl_Nposs  = analy(:,plotdat.Npossindex, ilevel,:);
-      plotdat.ges_Nused  = guess(:,plotdat.Nusedindex, ilevel,:);
-      plotdat.anl_Nused  = analy(:,plotdat.Nusedindex, ilevel,:);
-      plotdat.Yrange     = FindRange(plotdat);
-      
+      plotdat.ges_Nqc5  = guess(:,plotdat.NQC5index  ,ilevel,:);
+      plotdat.anl_Nqc5  = analy(:,plotdat.NQC5index  ,ilevel,:);
+      fprintf(logfid,'DART QC == 5, prior/post %d %d\n',sum(plotdat.ges_Nqc5(:)), ...
+                                                 sum(plotdat.anl_Nqc5(:)));
+
+      plotdat.ges_Nqc6  = guess(:,plotdat.NQC6index  ,ilevel,:);
+      plotdat.anl_Nqc6  = analy(:,plotdat.NQC6index  ,ilevel,:);
+      fprintf(logfid,'DART QC == 6, prior/post %d %d\n',sum(plotdat.ges_Nqc6(:)), ...
+                                                 sum(plotdat.anl_Nqc6(:)));
+
+      plotdat.ges_Nqc7  = guess(:,plotdat.NQC7index  ,ilevel,:);
+      plotdat.anl_Nqc7  = analy(:,plotdat.NQC7index  ,ilevel,:);
+      fprintf(logfid,'DART QC == 7, prior/post %d %d\n',sum(plotdat.ges_Nqc7(:)), ...
+                                                 sum(plotdat.anl_Nqc7(:)));
+
+      plotdat.ges_Nposs = guess(:,plotdat.Npossindex, ilevel,:);
+      plotdat.anl_Nposs = analy(:,plotdat.Npossindex, ilevel,:);
+      fprintf(logfid,'# obs poss,   prior/post %d %d\n',sum(plotdat.ges_Nposs(:)), ...
+                                                        sum(plotdat.anl_Nposs(:)));
+
+      plotdat.ges_Nused = guess(:,plotdat.Nusedindex, ilevel,:);
+      plotdat.anl_Nused = analy(:,plotdat.Nusedindex, ilevel,:);
+      fprintf(logfid,'# obs used,   prior/post %d %d\n',sum(plotdat.ges_Nused(:)), ...
+                                                        sum(plotdat.anl_Nused(:)));
+
+      plotdat.ges_copy  = guess(:,plotdat.copyindex,  ilevel,:);
+      plotdat.anl_copy  = analy(:,plotdat.copyindex,  ilevel,:);
+      plotdat.ges_rmse  = guess(:,plotdat.rmseindex,  ilevel,:);
+      plotdat.anl_rmse  = analy(:,plotdat.rmseindex,  ilevel,:);
+
+      plotdat.Yrange    = FindRange(plotdat);
+
       % plot by region
 
       if (plotdat.nregions > 2)
@@ -234,32 +311,33 @@
    string_rmse  = sprintf('%s pr=%.5g, po=%.5g','rmse', mean_pr_rmse, mean_po_rmse);
    string_other = sprintf('%s pr=%.5g, po=%.5g', plotdat.copystring, ...
                           mean_pr_other, mean_po_other);
+   plotdat.subtitle = sprintf('%s     %s     %s',plotdat.myregion,string_rmse, string_other);
 
-   % Plot the bias and 'xxx' on the same (left) axis.
+   % Plot the rmse and 'xxx' on the same (left) axis.
    % The observation count will use the axis on the right.
-   % Ultimately, we want to suppress the 'auto' feature of the
-   % axis labelling, so we manually set some values that normally
+   % We want to suppress the 'auto' feature of the axis labelling, 
+   % so we manually set some values that normally
    % don't need to be set.
    
    ax1 = subplot(plotdat.nregions,1,plotdat.region);
 
-   h1 = plot(t,rmse,'k+-',t,other,'ro-');
-   set(h1,'LineWidth',plotdat.linewidth);
-   h = legend(string_rmse, string_other);
+   h1 = plot(t,rmse,'k+-',t,other,'ro-','LineWidth',plotdat.linewidth);
+   h = legend(h1,'rmse', plotdat.copystring);
    legend(h,'boxoff')
 
    axlims = axis;
    axlims = [axlims(1:2) plotdat.Yrange];
    axis(axlims)
 
+   plotdat.ylabel{1} = plotdat.myregion;
    switch lower(plotdat.copystring)
       case 'bias'
          % plot a zero-bias line
          h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
          set(h4,'LineWidth',1.5,'LineSTyle',':')
-         plotdat.ylabel = sprintf('%s (%s) and rmse',plotdat.copystring,plotdat.biasconv);
+         plotdat.ylabel{2} = sprintf('rmse and %s (%s)',plotdat.copystring,plotdat.biasconv);
       otherwise
-         plotdat.ylabel = sprintf('%s and rmse',plotdat.copystring);
+         plotdat.ylabel{2} = sprintf('rmse and %s',plotdat.copystring);
    end
    
    % hokey effort to decide to plot months/days vs. daynum vs.
@@ -268,14 +346,30 @@
    if ((plotdat.bincenters(1) > 1000) && (ttot > 5))
       datetick('x',6,'keeplimits','keepticks');
       monstr = datestr(plotdat.bincenters(1),21);
-      xlabel(sprintf('month/day - %s start',monstr))
+      xlabelstring = sprintf('month/day - %s start',monstr);
    elseif (plotdat.bincenters(1) > 1000)
       datetick('x',15,'keeplimits','keepticks')
       monstr = datestr(plotdat.bincenters(1),21);
-      xlabel(sprintf('%s start',monstr))
+      xlabelstring = sprintf('%s start',monstr);
    else
-      xlabel('days')
+      xlabelstring = 'days';
    end
+
+   % only put x axis on last/bottom plot
+   if (plotdat.region == plotdat.nregions)
+      xlabel(xlabelstring)
+   end
+
+   % more annotation ...
+
+   if (plotdat.region == 1)
+      title({plotdat.title, plotdat.subtitle}, ...
+         'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
+      BottomAnnotation(plotdat.fname)
+   else
+      title(plotdat.subtitle, 'Interpreter', 'none', ...
+         'Fontsize', 12, 'FontWeight', 'bold')
+   end
    
    % use same X,Y limits for all plots in this region
    nYticks = length(get(ax1,'YTick'));
@@ -301,7 +395,7 @@
    yinc   = (ylimits(2)-ylimits(1))/(nYticks-1);
    yticks = ylimits(1):yinc:ylimits(2);
    niceyticks = round(10*yticks')/10;
-   set(ax2,'XTick',get(ax1,'XTick'),'XTicklabel',get(ax1,'XTicklabel'), ...
+   set(ax2,'XTick',get(ax1,'XTick'),'XTicklabel',[], ...
            'YTick',          yticks,'YTicklabel',num2str(niceyticks))
        
    set(get(ax2,'Ylabel'),'String','# of obs : o=poss, +=used')
@@ -309,13 +403,6 @@
    set(ax1,'Position',get(ax2,'Position'))
    grid
 
-   if (plotdat.region == 1)
-   title({plotdat.title, plotdat.myregion}, ...
-         'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
-   else
-   title(plotdat.myregion, 'Interpreter', 'none', ...
-         'Fontsize', 12, 'FontWeight', 'bold')
-   end
 
 
 
@@ -323,13 +410,18 @@
 % annotates the directory containing the data being plotted
 subplot('position',[0.48 0.01 0.04 0.04])
 axis off
-bob = which(main);
-[pathstr,name,ext,versn] = fileparts(bob);
-h = text(0.0,0.5,pathstr);
+str = which(main);
+if ( isempty(str) )   % must be a relative filename 
+   mydir = pwd;
+   string1 = sprintf('data file: %s/%s',mydir,main);
+else
+   string1 = sprintf('data file: %s',str);
+end
+h = text(0.0, 0.5, string1);
 set(h,'HorizontalAlignment','center', ...
       'VerticalAlignment','middle',...
       'Interpreter','none',...
-      'FontSize',8)
+      'FontSize',10)
 
 
 


More information about the Dart-dev mailing list