[Dart-dev] [6324] DART/trunk/diagnostics/matlab: Implemented Kevin' s algorithm for computing nice ticks on the dual-axis plots.
nancy at ucar.edu
nancy at ucar.edu
Mon Jul 29 10:49:33 MDT 2013
Revision: 6324
Author: thoar
Date: 2013-07-29 10:49:33 -0600 (Mon, 29 Jul 2013)
Log Message:
-----------
Implemented Kevin's algorithm for computing nice ticks on the dual-axis plots.
Ticklabels are now whole numbers and are more robust in the face of axis rescaling.
Some of the labels have also moved around to be consistent across plot types.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
DART/trunk/diagnostics/matlab/plot_evolution.m
DART/trunk/diagnostics/matlab/plot_profile.m
DART/trunk/diagnostics/matlab/plot_rank_histogram.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
DART/trunk/diagnostics/matlab/two_experiments_evolution.m
DART/trunk/diagnostics/matlab/two_experiments_profile.m
Added Paths:
-----------
DART/trunk/diagnostics/matlab/matchingXticks.m
DART/trunk/diagnostics/matlab/matchingYticks.m
-------------- next part --------------
Added: DART/trunk/diagnostics/matlab/matchingXticks.m
===================================================================
--- DART/trunk/diagnostics/matlab/matchingXticks.m (rev 0)
+++ DART/trunk/diagnostics/matlab/matchingXticks.m 2013-07-29 16:49:33 UTC (rev 6324)
@@ -0,0 +1,94 @@
+function xscale = matchingXticks(ax1, ax2)
+%% This takes the existing X ticks from ax1 (presumed nice)
+% and determines the nice, matching labels for ax2.
+%
+% In our implementation (for _profile.m applications):
+% the data plotted on ax1 are the quantities of interest.
+% the data plotted on ax2 are the observation counts.
+%
+% xscale = matchingXticks(ax1, ax2)
+%
+% xscale the scaling factor between ticks and labels ...
+%
+% USAGE - given two coincident axes (see plot_profile.m for an example):
+%
+% xscale = matchingXticks(ax1,ax2);
+% set(get(ax2,'Xlabel'),'String',['# of obs (o=poss, +=used) x' int2str(uint32(xscale))])
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+Dlimits = get(ax1,'XLim');
+DXticks = get(ax1,'XTick');
+nXticks = length(DXticks);
+xlimits = get(ax2,'XLim');
+
+% make sure the #_obs ticks have whole numbers.
+% The log of the x range is useful because it separates the magnitude from the number.
+% The xlimits sometimes exaggerate actual range of obs values too much.
+
+xrange_log = log10(xlimits(2) - xlimits(1));
+
+% Separate the log of the range (magnitude.span) into the 2 parts.
+% E.g. xrange = 1500 means that xrange_log = 3.176,
+% so the size is 10^3 and the span is 10^.176 = 1.5
+
+xrange_magnitude = floor(xrange_log);
+xrange_span = 10^(xrange_log - xrange_magnitude);
+
+% Match top axis ticks to bottom axis ticks
+% Generate a nicely sized spacing between ticks, rounded up from the original minimum size.
+% In matlab 10^[log10(x)] = x + delta where delta can be > eps (2.e-16)
+
+delta = 1.0e-14;
+
+% Rescale the magnitude and span if the span is to small to be useful.
+
+if ((nXticks -1 - xrange_span) > delta )
+ xrange_span = xrange_span * 10 ;
+ xrange_magnitude = xrange_magnitude - 1;
+end
+
+% pass the scaling factor back for the axis label
+
+xscale = 10.^xrange_magnitude;
+
+% Here's the distance between obs ticks, in units of SCALED # of observations.
+% It's used to label the ticks and determine the new XLim.
+% Prevent fractional tick spans (when there are only a few obs, and relatively many ticks).
+
+xrange_tick_size = ceil((xrange_span - delta) / (nXticks - 1)) * xscale ;
+xrange_tick_size = max([xrange_tick_size 1.0]);
+
+% Set the tick values for the obs number axis
+
+xticks0 = floor(xlimits(1)/xrange_tick_size - 1) * xrange_tick_size;
+iticks = 1:nXticks;
+xticks = xticks0 + xrange_tick_size*iticks;
+
+% New axis left and right for obs number axis.
+% These are returned to the calling program.
+
+D_tick_size = DXticks(2) - DXticks(1);
+xlimits(1) = xticks( 1 ) - xrange_tick_size*((DXticks(1) - Dlimits( 1 ))/D_tick_size);
+xlimits(2) = xticks(nXticks) + xrange_tick_size*((Dlimits(2) - DXticks(nXticks))/D_tick_size);
+
+newticklabels = num2str(round((10/xscale)*xticks') /10);
+
+% use the new ticks and labels.
+
+set(ax2,'XTick', xticks, 'XTicklabel', newticklabels, 'XLim', xlimits)
+
+% sometimes the ax2 position changes from the new ticks. Make sure ax1 follows.
+
+set(ax1,'Position',get(ax2,'Position'))
+grid
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Property changes on: DART/trunk/diagnostics/matlab/matchingXticks.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/trunk/diagnostics/matlab/matchingYticks.m
===================================================================
--- DART/trunk/diagnostics/matlab/matchingYticks.m (rev 0)
+++ DART/trunk/diagnostics/matlab/matchingYticks.m 2013-07-29 16:49:33 UTC (rev 6324)
@@ -0,0 +1,85 @@
+function matchingYticks(ax1, ax2)
+%% This takes the existing Y ticks from ax1 (presumed nice)
+% and determines nice-looking matching labels for ax2.
+%
+% In our implementation (for _evolution.m applications):
+% the data plotted on ax1 are the quantities of interest.
+% the data plotted on ax2 are the observation counts.
+%
+% USAGE - given two coincident axes (see plot_evolution.m for an example):
+%
+% matchingYticks(ax1,ax2);
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+Dlimits = get(ax1,'YLim');
+DYticks = get(ax1,'YTick');
+nYticks = length(DYticks);
+ylimits = get(ax2,'YLim');
+
+% make sure the #_obs ticks have whole numbers.
+% The log of the y range is useful because it separates the magnitude from the number.
+% The ylimits sometimes exaggerate actual range of obs values too much.
+
+yrange_log = log10(ylimits(2) - ylimits(1));
+
+% Separate the log of the range (magnitude.span) into the 2 parts.
+% E.g. yrange = 1500 means that yrange_log = 3.176,
+% so the size is 10^3 and the span is 10^.176 = 1.5
+
+yrange_magnitude = floor(yrange_log);
+yrange_span = 10^(yrange_log - yrange_magnitude);
+
+% Match right axis ticks to left axis ticks
+% Generate a nicely sized spacing between ticks, rounded up from the original minimum size.
+% In matlab 10^[log10(x)] = x + delta where delta can be > eps (2.e-16)
+
+delta = 1.0e-14;
+
+% Rescale the magnitude and span if the span is to small to be useful.
+
+if ((nYticks -1 - yrange_span) > delta )
+ yrange_span = yrange_span * 10;
+ yrange_magnitude = yrange_magnitude - 1;
+end
+
+% Here's the distance between obs ticks, in units of # of observations.
+% It's used to label the ticks and determine the new YLim.
+% Prevent fractional tick spans (when there are only a few obs, and relatively many ticks).
+
+yrange_tick_size = ceil((yrange_span - delta) / (nYticks -1)) * (10^yrange_magnitude);
+yrange_tick_size = max([yrange_tick_size 1.0]);
+
+% Set the tick values for the obs number axis
+
+yticks0 = floor(ylimits(1)/yrange_tick_size -1) * yrange_tick_size;
+iticks = 1:nYticks;
+yticks = yticks0 + yrange_tick_size*iticks;
+
+% New axis bottom and top for obs number axis.
+% These are returned to the calling program.
+
+D_tick_size = DYticks(2) - DYticks(1);
+ylimits(1) = yticks( 1 ) - yrange_tick_size*((DYticks(1) - Dlimits( 1 ))/D_tick_size);
+ylimits(2) = yticks(nYticks) + yrange_tick_size*((Dlimits(2) - DYticks(nYticks))/D_tick_size);
+
+newticklabels = num2str(round(10*yticks')/10);
+
+% use the new ticks and labels.
+
+set(ax2,'YTick', yticks, 'YTickLabel', newticklabels,'YLim', ylimits)
+
+% sometimes the ax2 position changes from the new ticks. Make sure ax1 follows.
+
+set(ax1,'Position',get(ax2,'Position'))
+grid
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Property changes on: DART/trunk/diagnostics/matlab/matchingYticks.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2013-07-26 19:32:37 UTC (rev 6323)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2013-07-29 16:49:33 UTC (rev 6324)
@@ -8,7 +8,7 @@
%
% fname : netcdf file produced by 'obs_diag'
% copystring : 'copy' string == quantity of interest. These
-% can be any of the ones available in the netcdf
+% can be any of the ones available in the netcdf
% file 'CopyMetaData' variable.
% (ncdump -v CopyMetaData obs_diag_output.nc)
%
@@ -52,9 +52,14 @@
diminfo = nc_getdiminfo(fname,'region');
plotdat.nregions = diminfo.Length;
-region_names = nc_varget(fname,'region_names');
-plotdat.region_names = deblank(region_names);
+plotdat.region_names = nc_varget(fname,'region_names');
+% Matlab wants character matrices to be Nx1 instead of 1xN.
+
+if (plotdat.nregions == 1 && (size(plotdat.region_names,2) == 1) )
+ plotdat.region_names = deblank(plotdat.region_names');
+end
+
% Coordinate between time types and dates
calendar = nc_attget(fname,'time','calendar');
@@ -72,11 +77,11 @@
plotdat.toff = plotdat.binedges(1) + iskip;
plotdat.timespan = sprintf('%s through %s', datestr(plotdat.toff), ...
- datestr(max(plotdat.binedges(:))));
+ datestr(max(plotdat.binedges(:))));
% set up a structure with all static plotting components
-plotdat.xlabel = {sprintf('bias (%s) and %s',plotdat.biasconv,copystring), plotdat.timespan};
+plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copystring);
plotdat.linewidth = 2.0;
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
@@ -84,7 +89,7 @@
plotdat.nvars = length(plotdat.varnames);
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.biasindex = get_copy_index(fname,'bias');
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
@@ -94,267 +99,307 @@
%----------------------------------------------------------------------
for ivar = 1:plotdat.nvars
-
+
% create the variable names of interest.
-
+
plotdat.myvarname = plotdat.varnames{ivar};
plotdat.guessvar = sprintf('%s_VPguess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_VPanaly',plotdat.varnames{ivar});
-
+
% remove any existing postscript file - will simply append each
% level as another 'page' in the .ps file.
psfname = sprintf('%s_bias_%s_profile.ps',plotdat.varnames{ivar},plotdat.copystring);
- disp(sprintf('Removing %s from the current directory.',psfname))
+ fprintf('Removing %s from the current directory.\n',psfname)
system(sprintf('rm %s',psfname));
-
+
% get appropriate vertical coordinate variable
-
+
guessdims = nc_var_dims( fname, plotdat.guessvar);
analydims = nc_var_dims( fname, plotdat.analyvar);
varinfo = nc_getvarinfo(fname, plotdat.analyvar);
-
- if ( findstr('surface',guessdims{2}) > 0 )
+
+ % this is a superfluous check ... FindVerticalVars already weeds out
+ % variables only present on surface or undef because obs_diag
+ % does not time-average statistics for these.
+
+ if (~ isempty(strfind(guessdims{2},'surface')))
fprintf('%s is a surface field.\n',plotdat.guessvar)
fprintf('Cannot display a surface field this way.\n')
- elseif ( findstr('undef',guessdims{2}) > 0 )
+ continue
+ elseif (~ isempty(strfind(guessdims{2},'undef')))
fprintf('%s has no vertical definition.\n',plotdat.guessvar)
fprintf('Cannot display this field this way.\n')
+ continue
end
-
+
[level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
plotdat.level_org = level_org;
plotdat.level_units = level_units;
plotdat.nlevels = nlevels;
plotdat.level_edges = level_edges;
plotdat.Yrange = Yrange;
-
+
% Matlab likes strictly ASCENDING order for things to be plotted,
% then you can impose the direction. The data is stored in the original
% order, so the sort indices are saved to reorder the data.
-
+
if (plotdat.level_org(1) > plotdat.level_org(plotdat.nlevels))
plotdat.YDir = 'reverse';
else
plotdat.YDir = 'normal';
end
+
+ % Add error-checking for output from older versions of obs_diag.
+
[levels, indices] = sort(plotdat.level_org);
- plotdat.level = levels;
+ plotdat.level = unique(levels);
+ if (length(plotdat.level) ~= length(levels))
+ error('There is a duplicated value in the array specifying the levels - must change your input.nml and rerun obs_diag')
+ end
+
plotdat.indices = indices;
level_edges = sort(plotdat.level_edges);
plotdat.level_edges = level_edges;
- guess = nc_varget(fname, plotdat.guessvar);
- analy = nc_varget(fname, plotdat.analyvar);
+ guess = nc_varget(fname, plotdat.guessvar);
+ analy = nc_varget(fname, plotdat.analyvar);
n = size(analy);
-
+
% singleton dimensions are auto-squeezed - which is unfortunate.
% We want these things to be 3D. [copy-level-region]
% Sometimes there is one region, sometimes one level, ...
% To complicate matters, the stupid 'ones' function does not allow
% the last dimension to be unity ... so you have double the size
% of the array ...
-
+
if ( plotdat.nregions == 1 )
- bob = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
- ted = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
- bob(:,:,1) = guess;
- ted(:,:,1) = analy;
- guess = bob; clear bob
- analy = ted; clear ted
+ bob = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
+ ted = NaN*ones(varinfo.Size(1),varinfo.Size(2),2);
+ bob(:,:,1) = guess;
+ ted(:,:,1) = analy;
+ guess = bob; clear bob
+ analy = ted; clear ted
elseif ( plotdat.nlevels == 1 )
- bob = NaN*ones(varinfo.Size);
- ted = NaN*ones(varinfo.Size);
- bob(:,1,:) = guess;
- ted(:,1,:) = analy;
- guess = bob; clear bob
- analy = ted; clear ted
+ bob = NaN*ones(varinfo.Size);
+ ted = NaN*ones(varinfo.Size);
+ bob(:,1,:) = guess;
+ ted(:,1,:) = analy;
+ guess = bob; clear bob
+ analy = ted; clear ted
end
% check to see if there is anything to plot
+
nposs = sum(guess(plotdat.Npossindex,:,:));
-
+
if ( sum(nposs(:)) < 1 )
fprintf('No obs for %s... skipping\n', plotdat.varnames{ivar})
continue
end
-
+
plotdat.ges_copy = guess(plotdat.copyindex, :, :);
plotdat.anl_copy = analy(plotdat.copyindex, :, :);
plotdat.ges_bias = guess(plotdat.biasindex, :, :);
plotdat.anl_bias = analy(plotdat.biasindex, :, :);
-
+
plotdat.ges_Nposs = guess(plotdat.Npossindex, :, :);
plotdat.anl_Nposs = analy(plotdat.Npossindex, :, :);
plotdat.ges_Nused = guess(plotdat.Nusedindex, :, :);
plotdat.anl_Nused = guess(plotdat.Nusedindex, :, :);
plotdat.Xrange = FindRange(plotdat);
-
+
% plot by region
-
- clf; orient tall
-
+
+ clf; orient tall; wysiwyg
+
for iregion = 1:plotdat.nregions
- plotdat.region = iregion;
+ plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
-
+
myplot(plotdat);
end
-
+
if (plotdat.nregions > 2)
CenterAnnotation(plotdat.myvarname)
end
- BottomAnnotation(fname)
-
+ BottomAnnotation(plotdat.timespan, fname)
+
% create a postscript file
print(gcf,'-dpsc','-append',psfname);
-
+
end
-%----------------------------------------------------------------------
+
+
+%=====================================================================
% 'Helper' functions
-%----------------------------------------------------------------------
+%=====================================================================
+
+
function myplot(plotdat)
- % Interlace the [ges,anl] to make a sawtooth plot.
- % By this point, the middle two dimensions are singletons.
- % The data must be sorted to match the order of the levels.
- cg = plotdat.ges_copy(:,:,plotdat.region); CG = cg(plotdat.indices);
- ca = plotdat.anl_copy(:,:,plotdat.region); CA = ca(plotdat.indices);
+%% Interlace the [ges,anl] to make a sawtooth plot.
+% By this point, the middle two dimensions are singletons.
+% The data must be sorted to match the order of the levels.
+cg = plotdat.ges_copy(:,:,plotdat.region); CG = cg(plotdat.indices);
+ca = plotdat.anl_copy(:,:,plotdat.region); CA = ca(plotdat.indices);
- mg = plotdat.ges_bias(:,:,plotdat.region); MG = mg(plotdat.indices);
- ma = plotdat.anl_bias(:,:,plotdat.region); MA = ma(plotdat.indices);
+mg = plotdat.ges_bias(:,:,plotdat.region); MG = mg(plotdat.indices);
+ma = plotdat.anl_bias(:,:,plotdat.region); MA = ma(plotdat.indices);
- g = plotdat.ges_Nposs(:,:,plotdat.region); G = g(plotdat.indices);
- a = plotdat.anl_Nposs(:,:,plotdat.region); A = a(plotdat.indices);
- nobs_poss = G;
- nposs_delta = G - A;
+g = plotdat.ges_Nposs(:,:,plotdat.region); G = g(plotdat.indices);
+a = plotdat.anl_Nposs(:,:,plotdat.region); A = a(plotdat.indices);
+nobs_poss = G;
+nposs_delta = G - A;
- g = plotdat.ges_Nused(:,:,plotdat.region); G = g(plotdat.indices);
- a = plotdat.anl_Nused(:,:,plotdat.region); A = a(plotdat.indices);
- nobs_used = G;
- nused_delta = G - A;
+g = plotdat.ges_Nused(:,:,plotdat.region); G = g(plotdat.indices);
+a = plotdat.anl_Nused(:,:,plotdat.region); A = a(plotdat.indices);
+nobs_used = G;
+nused_delta = G - A;
- % Determine some quantities for the legend
- nobs = sum(nobs_used);
- if ( nobs > 1 )
- bias_guess = mean(MG(isfinite(MG)));
- bias_analy = mean(MA(isfinite(MA)));
- other_guess = mean(CG(isfinite(CG)));
- other_analy = mean(CA(isfinite(CA)));
- else
- bias_guess = NaN;
- bias_analy = NaN;
- other_guess = NaN;
- other_analy = NaN;
- end
+% Determine some quantities for the legend
+nobs = sum(nobs_used);
+if ( nobs > 1 )
+ bias_guess = mean(MG(isfinite(MG)));
+ bias_analy = mean(MA(isfinite(MA)));
+ other_guess = mean(CG(isfinite(CG)));
+ other_analy = mean(CA(isfinite(CA)));
+else
+ bias_guess = NaN;
+ bias_analy = NaN;
+ other_guess = NaN;
+ other_analy = NaN;
+end
- str_bias_pr = sprintf('%s pr=%.5g','bias',bias_guess);
- str_bias_po = sprintf('%s po=%.5g','bias',bias_analy);
- str_other_pr = sprintf('%s pr=%.5g',plotdat.copystring,other_guess);
- str_other_po = sprintf('%s po=%.5g',plotdat.copystring,other_analy);
+str_bias_pr = sprintf('%s pr=%.5g','bias',bias_guess);
+str_bias_po = sprintf('%s po=%.5g','bias',bias_analy);
+str_other_pr = sprintf('%s pr=%.5g',plotdat.copystring,other_guess);
+str_other_po = sprintf('%s po=%.5g',plotdat.copystring,other_analy);
- % Plot the bias and 'xxx' on the same (bottom) axis.
- % The observation count will use the axis on the top.
- % Ultimately, 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.
-
- % if more then 4 regions, this will not work (well) ...
- if ( plotdat.nregions > 2 )
- ax1 = subplot(2,2,plotdat.region);
- else
- ax1 = subplot(1,plotdat.nregions,plotdat.region);
- axpos = get(ax1,'Position');
- axpos(4) = 0.925*axpos(4);
- set(ax1,'Position',axpos);
- end
+% Plot the bias and 'xxx' on the same (bottom) axis.
+% The observation count will use the axis on the top.
+% Ultimately, 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.
- Stripes(plotdat.Xrange, plotdat.level_edges);
- set(ax1,'YDir', plotdat.YDir,'YTick',plotdat.level,'Layer','top')
- ylabel(plotdat.level_units)
+% if more then 4 regions, this will not work (well) ...
+if ( plotdat.nregions > 2 )
+ ax1 = subplot(2,2,plotdat.region);
+else
+ % subplot(1,N,N) will plot 1 row of N subplots
+ ax1 = subplot(1,plotdat.nregions,plotdat.region);
+ axpos = get(ax1,'Position');
+ axpos(4) = 0.925*axpos(4); % make room for title
+ set(ax1,'Position',axpos);
+end
- %% draw the result of the experiment
+% add type of vertical coordinate info for adjusting axes to accomodate legend
- hold on;
- h1 = plot(MG,plotdat.level,'k+-',MA,plotdat.level,'k+:', ...
- CG,plotdat.level,'ro-',CA,plotdat.level,'ro:');
- set(h1,'LineWidth',plotdat.linewidth);
+Stripes(plotdat.Xrange, plotdat.level_edges, plotdat.level_units, plotdat.nregions);
+set(ax1, 'YDir', plotdat.YDir, 'YTick', plotdat.level, 'Layer', 'top')
+ylabel(plotdat.level_units)
- biasline = line([0 0],plotdat.Yrange,'Color','k','Parent',ax1);
- set(biasline,'LineWidth',2.0,'LineStyle','-.')
+% draw the result of the experiment
- hold off;
- h = legend(h1, str_bias_pr, str_bias_po, str_other_pr, str_other_po);
- set(h,'Interpreter','none','Box','off')
+hold on;
+h1 = plot(MG,plotdat.level,'k+-',MA,plotdat.level,'k+:', ...
+ CG,plotdat.level,'ro-',CA,plotdat.level,'ro:');
+set(h1,'LineWidth',plotdat.linewidth);
+hold off;
- %% Create another axes to use for plotting the observation counts
+zeroline = line([0 0],plotdat.Yrange,'Color','k','Parent',ax1);
+set(zeroline,'LineWidth',1.0,'LineStyle','-.')
- ax2 = axes('position',get(ax1,'Position'), ...
- 'XAxisLocation','top', ...
- 'YAxisLocation','right',...
- 'Color','none',...
- 'XColor','b','YColor','b',...
- 'YLim',get(ax1,'YLim'), ...
- 'YDir',get(ax1,'YDir'));
+h = legend(h1, str_bias_pr, str_bias_po, str_other_pr, str_other_po, 'Location', 'NorthWest');
+set(h,'Interpreter','none','Box','off')
- h2 = line(nobs_poss,plotdat.level,'Color','b','Parent',ax2);
- h3 = line(nobs_used,plotdat.level,'Color','b','Parent',ax2);
- set(h2,'LineStyle','none','Marker','o');
- set(h3,'LineStyle','none','Marker','+');
+% Create another axes to use for plotting the observation counts
- % use same Y ticks
- set(ax2,'YTick', get(ax1,'YTick'), ...
- 'YTicklabel',get(ax1,'YTicklabel'));
+ax2 = axes('position',get(ax1,'Position'), ...
+ 'XAxisLocation','top', ...
+ 'YAxisLocation','right',...
+ 'Color','none',...
+ 'XColor','b','YColor','b',...
+ 'YLim',get(ax1,'YLim'), ...
+ 'YDir',get(ax1,'YDir'));
- % use the same X ticks, but find the right label values
- [xticks, newticklabels] = matchingXticks(ax1,ax2);
- set(ax2,'XTick', xticks, 'XTicklabel', newticklabels)
+h2 = line(nobs_poss,plotdat.level,'Color','b','Parent',ax2);
+h3 = line(nobs_used,plotdat.level,'Color','b','Parent',ax2);
+set(h2,'LineStyle','none','Marker','o');
+set(h3,'LineStyle','none','Marker','+');
- set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
- set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
+% use same Y ticks
+set(ax2,'YTick',get(ax1,'YTick'), ...
+ 'YTicklabel',get(ax1,'YTicklabel'));
- if (plotdat.nregions <=2 )
- title({plotdat.myvarname, plotdat.myregion}, ...
- 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
- else
- title(plotdat.myregion, ...
- 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
- end
+% use the same X ticks, but find the right label values
+xscale = matchingXticks(ax1,ax2);
+set(get(ax2,'Xlabel'),'String',['# of obs (o=poss, +=used) x' int2str(uint32(xscale))])
+set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
+if (plotdat.nregions <=2 )
+ title({plotdat.myvarname, plotdat.myregion}, ...
+ 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
+else
+ title(plotdat.myregion, ...
+ 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
+end
+
+%=====================================================================
+
+
function CenterAnnotation(main)
-subplot('position',[0.48 0.48 0.04 0.04])
+%% annotates the variable being plotted
+subplot('position',[0.48 0.49 0.04 0.04])
axis off
h = text(0.5,0.5,main);
set(h,'HorizontalAlignment','center', ...
- 'VerticalAlignment','bottom', ...
- 'Interpreter','none', ...
- 'FontSize',12, ...
- 'FontWeight','bold')
+ 'VerticalAlignment','bottom', ...
+ 'Interpreter','none', ...
+ 'FontSize',12, ...
+ 'FontWeight','bold')
+%=====================================================================
-function BottomAnnotation(main)
-% annotates the filename containing the data being plotted
+
+function BottomAnnotation(timespan,main)
+%% annotates the timespan and the filename containing the data being plotted
subplot('position',[0.10 0.01 0.8 0.04])
axis off
-if ( main(1) == '/' ) % must be a absolute pathname
+
+if ( main(1) == '/' ) % must be an absolute pathname
string1 = sprintf('data file: %s',main);
else
mydir = pwd;
string1 = sprintf('data file: %s/%s',mydir,main);
end
-h = text(0.5, 0.5, string1);
-set(h, 'Interpreter', 'none', 'FontSize', 8);
-set(h, 'HorizontalAlignment','center');
+h = text(0.5, 0.67, timespan);
+set(h, 'HorizontalAlignment', 'center', ...
+ 'VerticalAlignment','middle', ...
+ 'Interpreter', 'none', ...
+ 'FontSize', 8);
+h = text(0.5, 0.33, string1);
+set(h, 'HorizontalAlignment', 'center', ...
+ 'VerticalAlignment','middle', ...
+ 'Interpreter', 'none', ...
+ 'FontSize', 8);
+
+%=====================================================================
+
+
function [y,ydims] = FindVerticalVars(x)
-% Returns UNIQUE (i.e. base) vertical variable names
+%% Returns UNIQUE (i.e. base) vertical variable names
+% In this context, if the variable has a 'time' dimension
+% it cannot be a variable of interest.
+
if ( ~(isfield(x,'allvarnames') && isfield(x,'allvardims')))
error('Doh! no ''allvarnames'' and ''allvardims'' components')
end
@@ -364,29 +409,31 @@
basedims = struct([]);
for i = 1:length(x.allvarnames)
- indx = findstr('time',x.allvardims{i});
- if (isempty(indx))
+ dimnames = lower(x.allvardims{i});
+ if (isempty(strfind(dimnames,'time')))
j = j + 1;
-
+
basenames{j} = ReturnBase(x.allvarnames{i});
basedims{j} = x.allvardims{i};
end
end
-[b,i,j] = unique(basenames);
+[~,i,j] = unique(basenames);
y = struct([]);
ydims = struct([]);
for k = 1:length(i)
fprintf('%2d is %s\n',k,basenames{i(k)})
- y{k} = basenames{i(k)};
-ydims{k} = basedims{i(k)};
+ y{k} = basenames{i(k)};
+ ydims{k} = basedims{i(k)};
end
+%=====================================================================
+
function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname)
-% Find the vertical dimension and harvest some info
+%% Find the vertical dimension and harvest some info
varinfo = nc_getvarinfo(fname,varname);
leveldim = [];
@@ -408,39 +455,49 @@
Yrange = [min(level_edges) max(level_edges)];
+%=====================================================================
+
+
function s = ReturnBase(string1)
-inds = findstr('_guess',string1);
+%% Pick off the variable name.
+s = [];
+inds = strfind(string1,'_guess');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-inds = findstr('_analy',string1);
+inds = strfind(string1,'_analy');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-inds = findstr('_VPguess',string1);
+inds = strfind(string1,'_VPguess');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-inds = findstr('_VPanaly',string1);
+inds = strfind(string1,'_VPanaly');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
+%=====================================================================
+
function x = FindRange(y)
-% Trying to pick 'nice' limits for plotting.
+%% Trying to pick 'nice' limits for plotting.
% Completely ad hoc ... and not well posed.
%
% In this scope, y is bounded from below by 0.0
%
-% If the numbers are very small ...
+% If the numbers are very small ...
-bob = [y.ges_copy(:) ; y.ges_bias(:); ...
- y.anl_copy(:) ; y.anl_bias(:)];
+bob = [y.ges_copy(:) ; y.ges_bias(:); y.anl_copy(:) ; y.anl_bias(:)];
inds = find(isfinite(bob));
if ( isempty(inds) )
@@ -449,86 +506,89 @@
glommed = bob(inds);
ymin = min(glommed);
ymax = max(glommed);
-
- if ( ymax > 1.0 )
+
+ if ( ymax > 1.0 )
ymin = floor(min(glommed));
ymax = ceil(max(glommed));
end
-
+
if (ymin == 0 && ymax == 0)
- ymax = 1;
+ ymax = 1;
end
if (ymin == ymax)
- ymin = ymin - 0.1*ymin;
- ymax = ymax + 0.1*ymax;
+ ymin = ymin - 0.1*ymin;
+ ymax = ymax + 0.1*ymax;
end
-
+
Yrange = [ymin ymax];
-
+
% Make sure a zero bias is visible on plot
- if ymax < 0
- Yrange = [ ymin 0.0 ];
- elseif ymin > 0
- Yrange = [ 0.0 ymax ];
- end
+ if ymax < 0
+ Yrange = [ ymin 0.0 ];
+ elseif ymin > 0
+ Yrange = [ 0.0 ymax ];
+ end
+
x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
end
+%=====================================================================
-function h = Stripes(x,edges)
+
+function h = Stripes(x,edges,units,nregions)
%% plot the subtle background stripes that indicate the vertical
% extent of what was averaged.
%
% FIXME:
% This really should be modified to add a percentage of the data
% range to provide space for the legend. Right now it is hardwired
-% to assume that we are plotting hPa, on a 'reverse' axis.
+% to assume that we are plotting hPa, on a 'reverse' axis.
+% kdr axlims(3) should be conditional on the observation vertical coordinate:
+% values for pressure coordinates are inappropriate for height coord.
+% It also assumes 4 plots/page, but 2 works better for plotting all levels of CAM5.
+% That requires a smaller % of vertical range for the legend.
-hold on;
-
% plot two little dots at the corners to make Matlab
% choose some plot limits. Given those nice limits and
% tick labels ... KEEP THEM. Later, make the dots invisible.
h = plot([min(x) max(x)],[min(edges) max(edges)]);
-axlims = axis;
-axlims(4) = max(edges);
-axlims(3) = -100; % This gives extra space for the legend.
+axlims = axis;
+region_factor = floor((nregions+1)/2.0);
+legend_fraction = 0.1375;
+
+% partial fix to legend space; add in option for vert coord = height.
+
+switch lower(units)
+ case 'hpa'
+ axlims(4) = max(edges);
+ axlims(3) = min(edges) - region_factor*legend_fraction*(axlims(4)-min(edges));
+ case 'm'
+ axlims(3) = min(edges);
+ axlims(4) = max(edges) + region_factor*legend_fraction*(max(edges)-axlims(3));
+ otherwise
+end
axis(axlims)
+% set up list of x,y values defining corner of every other stripe.
+
xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ];
+hold on;
for i = 1:2:(length(edges)-1)
- yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
- hf = fill(xc,yc,[0.8 0.8 0.8],'EdgeColor','none');
+ yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
+ hf = fill(xc,yc,[0.8 0.8 0.8],'EdgeColor','none');
end
-set(gca,'XGrid','on')
hold off;
-set(h,'Visible','off')
+set(gca,'XGrid','on')
+set(h,'Visible','off') % make the dots invisible
-function [xticks newticklabels] = matchingXticks(ax1, ax2)
-%% This takes the existing X ticks from ax1 (presumed nice)
-% and determines the matching labels for ax2 so we can keep
-% at least one of the axes looking nice.
-
-Dlimits = get(ax1,'XLim');
-DXticks = get(ax1,'XTick');
-nXticks = length(DXticks);
-xlimits = get(ax2,'XLim');
-
-slope = (xlimits(2) - xlimits(1))/(Dlimits(2) - Dlimits(1));
-xtrcpt = xlimits(2) -slope*Dlimits(2);
-
-xticks = slope*DXticks + xtrcpt;
-newticklabels = num2str(round(10*xticks')/10);
-
-
% <next few lines under version control, do not edit>
% $URL$
+% $Revision$
% $Date$
-% $Revision$
Modified: DART/trunk/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_evolution.m 2013-07-26 19:32:37 UTC (rev 6323)
+++ DART/trunk/diagnostics/matlab/plot_evolution.m 2013-07-29 16:49:33 UTC (rev 6324)
@@ -27,6 +27,7 @@
% about how many observations were assimilated, how many were available, etc.
% Both of these filenames contain the observation type as part of the name.
%
+%
% 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'
@@ -77,19 +78,20 @@
plotdat.region_names = nc_varget(fname,'region_names');
% Matlab wants character matrices to be Nx1 instead of 1xN.
+
if (plotdat.nregions == 1 && (size(plotdat.region_names,2) == 1) )
plotdat.region_names = deblank(plotdat.region_names');
end
-dimensionality = nc_attget(fname, nc_global, 'LocationRank');
-plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
-plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
-time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
-plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
-plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
-plotdat.latlim1 = local_nc_attget(fname, nc_global, 'latlim1');
-plotdat.latlim2 = local_nc_attget(fname, nc_global, 'latlim2');
+dimensionality = local_nc_attget(fname, nc_global, 'LocationRank');
+plotdat.binseparation = local_nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = local_nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = local_nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.lonlim1 = local_nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = local_nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = local_nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = local_nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = local_nc_attget(fname, nc_global, 'bias_convention');
% Coordinate between time types and dates
@@ -97,8 +99,14 @@
timeunits = nc_attget(fname,'time','units');
timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
timeorigin = datenum(timebase(1),timebase(2),timebase(3));
-skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
-iskip = time_to_skip(3) + skip_seconds/86400;
+if ( isempty(time_to_skip) == 1)
+ iskip = 0;
+elseif ( numel(time_to_skip) == 6)
+ skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
+ iskip = time_to_skip(3) + skip_seconds/86400;
+else
+ error('time_to_skip variable has unusual length. Should be either 0 or 6.')
+end
plotdat.bincenters = plotdat.bincenters + timeorigin;
plotdat.binedges = plotdat.binedges + timeorigin;
@@ -118,7 +126,6 @@
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');
@@ -186,6 +193,7 @@
plotdat.nlevels, plotdat.nregions);
% check to see if there is anything to plot
+
nposs = sum(guess(:,plotdat.Npossindex,:,:));
if ( sum(nposs(:)) < 1 )
@@ -234,9 +242,9 @@
% plot by region
if (plotdat.nregions > 2)
- clf; orient tall
+ clf; orient tall; wysiwyg
else
- clf; orient landscape
+ clf; orient landscape; wysiwyg
end
for iregion = 1:plotdat.nregions
@@ -265,13 +273,14 @@
end
end
-%----------------------------------------------------------------------
+%=====================================================================
% 'Helper' functions
-%----------------------------------------------------------------------
+%=====================================================================
+
function myplot(plotdat)
-% Interlace the [ges,anl] to make a sawtooth plot.
+%% The prior and posterior are plotted as separate items.
% By this point, the middle two dimensions are singletons.
cg = plotdat.ges_copy(:,:,:,plotdat.region);
ca = plotdat.anl_copy(:,:,:,plotdat.region);
@@ -311,9 +320,13 @@
ax1 = subplot(plotdat.nregions,1,plotdat.region);
h1 = plot(tg,cg,'k+-',ta,ca,'ro-','LineWidth',plotdat.linewidth);
-h = legend('forecast', 'analysis');
+h = legend('forecast', 'analysis');
legend(h,'boxoff','Interpreter','none')
+% get the range of the existing axis and replace with
+% replace y axis values
+% reset the axes limits
+
axlims = axis;
axlims = [axlims(1:2) plotdat.Yrange];
axis(axlims)
@@ -322,8 +335,8 @@
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',':')
+ zeroline = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
+ set(zeroline,'LineWidth',1.5,'LineSTyle',':')
plotdat.ylabel{2} = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
otherwise
plotdat.ylabel{2} = sprintf('%s',plotdat.copystring);
@@ -350,14 +363,13 @@
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')
+ title(plotdat.subtitle, ...
+ 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
end
% create a separate scale for the number of observations
@@ -375,19 +387,18 @@
set(ax2,'XTick',get(ax1,'XTick'),'XTickLabel',[])
% use the same Y ticks, but find the right label values
-[yticks, newticklabels] = matchingYticks(ax1,ax2);
-set(ax2,'YTick', yticks, 'YTickLabel', newticklabels)
+% Make sure that values at ticks are whole numbers.
+matchingYticks(ax1,ax2);
set(get(ax2,'Ylabel'),'String','# of obs : o=poss, +=used')
set(get(ax1,'Ylabel'),'String',plotdat.ylabel,'Interpreter','none')
-set(ax1,'Position',get(ax2,'Position'))
-grid
+%=====================================================================
function BottomAnnotation(main)
-% annotates the directory containing the data being plotted
+%% annotates the full path of the file being plotted
subplot('position',[0.48 0.01 0.04 0.04])
axis off
fullname = which(main); % Could be in MatlabPath
@@ -406,12 +417,14 @@
set(h,'HorizontalAlignment','center', ...
'VerticalAlignment','middle',...
'Interpreter','none',...
- 'FontSize',10)
+ 'FontSize',8)
+%=====================================================================
+
function [y,ydims] = FindTemporalVars(x)
-% Returns UNIQUE (i.e. base) temporal variable names
+%% Returns UNIQUE (i.e. base) temporal variable names
if ( ~(isfield(x,'allvarnames') && isfield(x,'allvardims')))
error('Doh! no ''allvarnames'' and ''allvardims'' components')
end
@@ -424,7 +437,7 @@
j = j + 1;
basenames{j} = ReturnBase(x.allvarnames{i});
- basedims{j} = x.allvardims{i};
+ basedims{ j} = x.allvardims{i};
end
end
@@ -433,37 +446,43 @@
ydims = cell(length(i),1);
for k = 1:length(i)
fprintf('%2d is %s\n',k,basenames{i(k)})
- y{k} = basenames{i(k)};
- ydims{k} = basedims{i(k)};
+ y{k} = basenames{i(k)};
+ ydims{k} = basedims{ i(k)};
end
+%=====================================================================
+
function s = ReturnBase(string1)
+%% Pick off the variable name.
inds = strfind(string1,'_guess');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-
inds = strfind(string1,'_analy');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-
inds = strfind(string1,'_VPguess');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-
inds = strfind(string1,'_VPanaly');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
+%=====================================================================
+
function x = FindRange(y)
-% Trying to pick 'nice' limits for plotting.
+%% Trying to pick 'nice' limits for plotting.
% Completely ad hoc ... and not well posed.
%
% In this scope, y is bounded from below by 0.0
@@ -471,7 +490,7 @@
% If the numbers are very small ...
bob = [y.ges_copy(:) ; ...
- y.anl_copy(:)];
+ y.anl_copy(:)];
inds = find(isfinite(bob));
if ( isempty(inds) )
@@ -494,25 +513,9 @@
end
+%=====================================================================
-function [yticks newticklabels] = matchingYticks(ax1, ax2)
-%% This takes the existing Y ticks from ax1 (presumed nice)
-% and determines the matching labels for ax2 so we can keep
-% at least one of the axes looking nice.
-Dlimits = get(ax1,'YLim');
-DYticks = get(ax1,'YTick');
-nYticks = length(DYticks);
-ylimits = get(ax2,'YLim');
-
-slope = (ylimits(2) - ylimits(1))/(Dlimits(2) - Dlimits(1));
-xtrcpt = ylimits(2) - slope*Dlimits(2);
-
-yticks = slope*DYticks + xtrcpt;
-newticklabels = num2str(round(10*yticks')/10);
-
-
-
function value = local_nc_varget(fname,varname)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list