[Dart-dev] [7331] DART/trunk/diagnostics/matlab: These scripts subtract the number of observations with DART QCs of 5 and 6
nancy at ucar.edu
nancy at ucar.edu
Wed Dec 24 14:22:45 MST 2014
Revision: 7331
Author: thoar
Date: 2014-12-24 14:22:45 -0700 (Wed, 24 Dec 2014)
Log Message:
-----------
These scripts subtract the number of observations with DART QCs of 5 and 6
from the number of observations possible. QCs of 5,6 are namelist settings
that FORCE the observation to not be assimilated, so this is a more meaningful
representation of the possible observations.
They also plot only 1 region per figure instead of trying to get clever
with 2 or 3 or 4 or ???
The line types and symbols are now consistent across the routines.
The number of observations used is now represented by an asterisk
instead of a plus symbol to avoid confusion with the plus symbol
occasionally used as part of some lines.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
DART/trunk/diagnostics/matlab/plot_evolution.m
DART/trunk/diagnostics/matlab/plot_obs_netcdf.m
DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m
DART/trunk/diagnostics/matlab/plot_profile.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
DART/trunk/diagnostics/matlab/read_obs_netcdf.m
DART/trunk/diagnostics/matlab/two_experiments_evolution.m
DART/trunk/diagnostics/matlab/two_experiments_profile.m
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2014-12-24 20:54:42 UTC (rev 7330)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2014-12-24 21:22:45 UTC (rev 7331)
@@ -24,11 +24,17 @@
%
% DART $Id$
+%%--------------------------------------------------------------------
+% Decode,Parse,Check the input
+%---------------------------------------------------------------------
+
if (exist(fname,'file') ~= 2)
error('file/fname <%s> does not exist',fname)
end
+%%--------------------------------------------------------------------
% Harvest plotting info/metadata from netcdf file.
+%---------------------------------------------------------------------
plotdat.fname = fname;
plotdat.copystring = copystring;
@@ -71,58 +77,51 @@
skip_seconds = timefloats(4)*3600 + timefloats(5)*60 + timefloats(6);
iskip = timefloats(3) + skip_seconds/86400.0;
+% Set up a structure to use for plotting
+
plotdat.bincenters = plotdat.bincenters + timeorigin;
plotdat.binedges = plotdat.binedges + timeorigin;
plotdat.Nbins = length(plotdat.bincenters);
plotdat.toff = plotdat.binedges(1) + iskip;
-
plotdat.timespan = sprintf('%s through %s', datestr(plotdat.toff), ...
- datestr(max(plotdat.binedges(:))));
+ datestr(max(plotdat.binedges(:))));
+plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copystring);
-% set up a structure with all static plotting components
-
-plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copystring);
-plotdat.linewidth = 2.0;
-
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindVerticalVars(plotdat);
-plotdat.nvars = length(plotdat.varnames);
+plotdat.nvars = length(plotdat.varnames);
+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');
+plotdat.NQC5index = get_copy_index(fname,'N_DARTqc_5');
+plotdat.NQC6index = get_copy_index(fname,'N_DARTqc_6');
-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');
+figuredata = setfigure();
-%----------------------------------------------------------------------
+%%---------------------------------------------------------------------
% Loop around (copy-level-region) observation types
%----------------------------------------------------------------------
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);
- 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);
-
+
% 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')
@@ -132,18 +131,18 @@
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
@@ -161,21 +160,21 @@
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);
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 = NaN*ones(varinfo.Size(1),varinfo.Size(2),1);
+ ted = NaN*ones(varinfo.Size(1),varinfo.Size(2),1);
bob(:,:,1) = guess;
ted(:,:,1) = analy;
guess = bob; clear bob
@@ -188,58 +187,65 @@
guess = bob; clear bob
analy = ted; clear ted
end
-
+
% check to see if there is anything to plot
+ % The number possible is decreased by the number of observations
+ % rejected by namelist control.
- nposs = sum(guess(plotdat.Npossindex,:,:));
-
+ fprintf('%d %s observations had DART QC of 5 (all regions).\n', ...
+ sum(sum(guess(plotdat.NQC5index, :,:))),plotdat.myvarname)
+ fprintf('%d %s observations had DART QC of 6 (all regions).\n', ...
+ sum(sum(guess(plotdat.NQC6index, :,:))),plotdat.myvarname)
+
+ nposs = sum(guess(plotdat.Npossindex,:,:)) - ...
+ sum(guess(plotdat.NQC5index ,:,:)) - ...
+ sum(guess(plotdat.NQC6index ,:,:));
+
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_Nqc5 = guess(plotdat.NQC5index, :, :);
+ plotdat.anl_Nqc5 = analy(plotdat.NQC5index, :, :);
+ plotdat.ges_Nqc6 = guess(plotdat.NQC6index, :, :);
+ plotdat.anl_Nqc6 = analy(plotdat.NQC6index, :, :);
plotdat.ges_Nused = guess(plotdat.Nusedindex, :, :);
plotdat.anl_Nused = guess(plotdat.Nusedindex, :, :);
+ plotdat.ges_Nposs = guess(plotdat.Npossindex, :, :) - ...
+ plotdat.ges_Nqc5 - plotdat.ges_Nqc6;
+ plotdat.anl_Nposs = analy(plotdat.Npossindex, :, :) - ...
+ plotdat.anl_Nqc5 - plotdat.anl_Nqc6;
plotdat.Xrange = FindRange(plotdat);
-
- % plot by region
-
- clf; orient tall; wysiwyg
-
+
+ % plot by region - each in its own figure.
+
for iregion = 1:plotdat.nregions
- plotdat.region = iregion;
+ figure(iregion); clf; orient(figuredata.orientation); wysiwyg
+ plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
-
- myplot(plotdat);
+ myplot(plotdat, figuredata);
+ BottomAnnotation(fname)
+
+ psfname = sprintf('%s_bias_%s_profile_region%d', ...
+ plotdat.varnames{ivar}, plotdat.copystring, iregion);
+ print(gcf,'-dpdf',psfname);
end
-
- if (plotdat.nregions > 2)
- CenterAnnotation(plotdat.myvarname)
- end
- BottomAnnotation(plotdat.timespan, fname)
-
- % create a postscript file
- print(gcf,'-dpsc','-append',psfname);
-
+
end
-
%=====================================================================
% 'Helper' functions
%=====================================================================
+function myplot(plotdat,figdata)
-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.
@@ -251,6 +257,7 @@
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;
@@ -284,33 +291,24 @@
% 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
- % 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
+ax1 = subplot('position',figdata.position);
% add type of vertical coordinate info for adjusting axes to accomodate legend
-Stripes(plotdat.Xrange, plotdat.level_edges, plotdat.level_units, plotdat.nregions);
+Stripes(plotdat.Xrange, plotdat.level_edges, plotdat.level_units);
set(ax1, 'YDir', plotdat.YDir, 'YTick', plotdat.level, 'Layer', 'top')
-ylabel(plotdat.level_units)
+set(ax1,'YAxisLocation','left','FontSize',figdata.fontsize)
% draw the result of the experiment
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);
+h1 = plot(MG,plotdat.level,'k+-',MA,plotdat.level,'k+--', ...
+ CG,plotdat.level,'ro-',CA,plotdat.level,'ro--');
+set(h1,'LineWidth',figdata.linewidth);
hold off;
-zeroline = line([0 0],plotdat.Yrange,'Color','k','Parent',ax1);
-set(zeroline,'LineWidth',1.0,'LineStyle','-.')
+zeroline = line([0 0],plotdat.Yrange,'Color',[0 100 0]/255,'Parent',ax1);
+set(zeroline,'LineWidth',2.5,'LineStyle','-')
h = legend(h1, str_bias_pr, str_bias_po, str_other_pr, str_other_po, 'Location', 'NorthWest');
set(h,'Interpreter','none','Box','off')
@@ -319,56 +317,41 @@
ax2 = axes('position',get(ax1,'Position'), ...
'XAxisLocation','top', ...
- 'YAxisLocation','right',...
- 'Color','none',...
- 'XColor','b','YColor','b',...
+ 'YAxisLocation','right', ...
+ 'Color','none', ...
+ 'XColor','b', ...
+ 'YColor',get(ax1,'YColor'), ...
'YLim',get(ax1,'YLim'), ...
- 'YDir',get(ax1,'YDir'));
+ 'YDir',get(ax1,'YDir'), ...
+ 'FontSize',get(ax1,'FontSize'));
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(h3,'LineStyle','none','Marker','*');
-% use same Y ticks
-set(ax2,'YTick',get(ax1,'YTick'), ...
- 'YTicklabel',get(ax1,'YTicklabel'));
+% use same Y ticks - but no labels.
+set(ax2,'YTick',get(ax1,'YTick'), 'YTicklabel',[]);
% 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')
+set(get(ax1,'Ylabel'),'String',plotdat.level_units, ...
+ 'Interpreter','none','FontSize',figdata.fontsize)
+set(get(ax1,'Xlabel'),'String',{plotdat.xlabel, plotdat.timespan}, ...
+ 'Interpreter','none','FontSize',figdata.fontsize)
+set(get(ax2,'Xlabel'),'String', ...
+ ['# of obs (o=poss, \ast=used) x' int2str(uint32(xscale))],'FontSize',figdata.fontsize)
-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
+title({plotdat.myregion, plotdat.myvarname}, ...
+ 'Interpreter', 'none', 'FontSize', figdata.fontsize, 'FontWeight', 'bold')
%=====================================================================
-function CenterAnnotation(main)
-%% 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')
-
-
-%=====================================================================
-
-
-function BottomAnnotation(timespan,main)
-%% annotates the timespan and the filename containing the data being plotted
+function BottomAnnotation(main)
+%% annotates the filename containing the data being plotted
subplot('position',[0.10 0.01 0.8 0.04])
axis off
@@ -379,12 +362,6 @@
string1 = sprintf('data file: %s/%s',mydir,main);
end
-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', ...
@@ -412,7 +389,7 @@
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
@@ -506,30 +483,30 @@
glommed = bob(inds);
ymin = min(glommed);
ymax = max(glommed);
-
+
if ( ymax > 1.0 )
ymin = floor(min(glommed));
ymax = ceil(max(glommed));
end
-
+
if (ymin == 0 && ymax == 0)
ymax = 1;
end
-
+
if (ymin == 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
-
+
x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
end
@@ -537,7 +514,7 @@
%=====================================================================
-function h = Stripes(x,edges,units,nregions)
+function h = Stripes(x,edges,units)
%% plot the subtle background stripes that indicate the vertical
% extent of what was averaged.
%
@@ -556,18 +533,17 @@
h = plot([min(x) max(x)],[min(edges) max(edges)]);
axlims = axis;
-region_factor = floor((nregions+1)/2.0);
-legend_fraction = 0.1375;
+legend_fraction = 0.22;
% 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));
+ axlims(3) = min(edges) - legend_fraction*(axlims(4)-min(edges));
case 'm'
axlims(3) = min(edges);
- axlims(4) = max(edges) + region_factor*legend_fraction*(max(edges)-axlims(3));
+ axlims(4) = max(edges) + legend_fraction*(max(edges)-axlims(3));
otherwise
end
axis(axlims)
@@ -587,6 +563,27 @@
set(h,'Visible','off') % make the dots invisible
+%=====================================================================
+
+
+function figdata = setfigure()
+%%
+% figure out a page layout
+% extra space at the bottom for the date/file annotation
+% extra space at the top because the titles have multiple lines
+
+orientation = 'tall';
+fontsize = 16;
+position = [0.15 0.12 0.7 0.75];
+linewidth = 2.0;
+
+figdata = struct('expcolors', {{'k','r','g','m','b','c','y'}}, ...
+ 'expsymbols', {{'o','s','d','p','h','s','*'}}, ...
+ 'prpolines', {{'-','--'}}, 'position', position, ...
+ 'fontsize',fontsize, 'orientation',orientation, ...
+ 'linewidth',linewidth);
+
+
% <next few lines under version control, do not edit>
% $URL$
% $Revision$
Modified: DART/trunk/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_evolution.m 2014-12-24 20:54:42 UTC (rev 7330)
+++ DART/trunk/diagnostics/matlab/plot_evolution.m 2014-12-24 21:22:45 UTC (rev 7331)
@@ -62,7 +62,9 @@
error('file/fname <%s> does not exist',fname)
end
+%%--------------------------------------------------------------------
% Harvest plotting info/metadata from netcdf file.
+%---------------------------------------------------------------------
plotdat.fname = fname;
plotdat.copystring = copystring;
@@ -108,15 +110,13 @@
error('time_to_skip variable has unusual length. Should be either 0 or 6.')
end
+% set up a structure with all static plotting components
+
plotdat.bincenters = plotdat.bincenters + timeorigin;
plotdat.binedges = plotdat.binedges + timeorigin;
plotdat.Nbins = length(plotdat.bincenters);
plotdat.toff = plotdat.bincenters(1) + iskip;
-% set up a structure with all static plotting components
-
-plotdat.linewidth = 2.0;
-
if (nvars == 0)
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindTemporalVars(plotdat);
@@ -134,7 +134,9 @@
plotdat.NQC6index = get_copy_index(fname,'N_DARTqc_6');
plotdat.NQC7index = get_copy_index(fname,'N_DARTqc_7');
-%----------------------------------------------------------------------
+figuredata = setfigure();
+
+%%---------------------------------------------------------------------
% Loop around (time-copy-level-region) observation types
%----------------------------------------------------------------------
@@ -149,9 +151,12 @@
% remove any existing postscript file - will simply append each
% level as another 'page' in the .ps file.
- psfname = sprintf('%s_%s_evolution.ps',plotdat.varnames{ivar},plotdat.copystring);
- fprintf('Removing %s from the current directory.\n',psfname)
- system(sprintf('rm %s',psfname));
+ for iregion = 1:plotdat.nregions
+ psfname{iregion} = sprintf('%s_%s_evolution_region%d.ps', ...
+ plotdat.varnames{ivar}, plotdat.copystring, iregion);
+ fprintf('Removing %s from the current directory.\n',psfname{iregion})
+ system(sprintf('rm %s',psfname{iregion}));
+ end
% remove any existing log file -
@@ -193,9 +198,21 @@
plotdat.nlevels, plotdat.nregions);
% check to see if there is anything to plot
-
- nposs = sum(guess(:,plotdat.Npossindex,:,:));
+ % The number possible is decreased by the number of observations
+ % rejected by namelist control.
+ nqc5 = guess(:,plotdat.NQC5index,:,:);
+ nqc6 = guess(:,plotdat.NQC6index,:,:);
+
+ fprintf('%d %s observations had DART QC of 5 (all levels, all regions).\n', ...
+ sum(nqc5(:)),plotdat.myvarname)
+ fprintf('%d %s observations had DART QC of 6 (all levels, all regions).\n', ...
+ sum(nqc6(:)),plotdat.myvarname)
+
+ nposs = sum(guess(:,plotdat.Npossindex,:,:)) - ...
+ sum(guess(:,plotdat.NQC5index ,:,:)) - ...
+ sum(guess(:,plotdat.NQC6index ,:,:));
+
if ( sum(nposs(:)) < 1 )
fprintf('%s no obs for %s... skipping\n', plotdat.varnames{ivar})
continue
@@ -224,8 +241,10 @@
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,:);
+ plotdat.ges_Nposs = guess(:,plotdat.Npossindex, ilevel,:) - ...
+ plotdat.ges_Nqc5 - plotdat.ges_Nqc6;
+ plotdat.anl_Nposs = analy(:,plotdat.Npossindex, ilevel,:) - ...
+ plotdat.anl_Nqc5 - plotdat.anl_Nqc6;
fprintf(logfid,'# obs poss, prior/post %d %d\n',sum(plotdat.ges_Nposs(:)), ...
sum(plotdat.anl_Nposs(:)));
@@ -239,15 +258,10 @@
plotdat.Yrange = FindRange(plotdat);
- % plot by region
+ % plot each region, each level to a separate figure
- if (plotdat.nregions > 2)
- clf; orient tall; wysiwyg
- else
- clf; orient landscape; wysiwyg
- end
-
for iregion = 1:plotdat.nregions
+ figure(iregion); clf; orient(figuredata.orientation); wysiwyg
plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
@@ -260,16 +274,16 @@
plotdat.level_units);
end
- myplot(plotdat);
- end
+ myplot(plotdat,figuredata);
- % create a postscript file
- print(gcf,'-dpsc','-append',psfname);
+ % create/append to the postscript file
+ print(gcf,'-dpsc','-append',psfname{iregion});
- % block to go slow and look at each one ...
- % disp('Pausing, hit any key to continue ...')
- % pause
+ % block to go slow and look at each one ...
+ % disp('Pausing, hit any key to continue ...')
+ % pause
+ end
end
end
@@ -278,7 +292,7 @@
%=====================================================================
-function myplot(plotdat)
+function myplot(plotdat,figdata)
%% The prior and posterior are plotted as separate items.
% By this point, the middle two dimensions are singletons.
@@ -309,7 +323,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);
+plotdat.subtitle = sprintf('%s %s',string_guess, string_analy);
% Plot the requested quantity on the left axis.
% The observation count will use the axis on the right.
@@ -317,11 +331,12 @@
% so we manually set some values that normally
% don't need to be set.
-ax1 = subplot(plotdat.nregions,1,plotdat.region);
+ax1 = subplot('position',figdata.position);
+set(ax1,'YAxisLocation','left','FontSize',figdata.fontsize)
-h1 = plot(tg,cg,'k+-',ta,ca,'ro-','LineWidth',plotdat.linewidth);
+h1 = plot(tg,cg,'k+-',ta,ca,'ro-','LineWidth',figdata.linewidth);
h = legend('forecast', 'analysis');
-legend(h,'boxoff','Interpreter','none')
+set(h,'Interpreter','none','Box','off')
% get the range of the existing axis and replace with
% replace y axis values
@@ -331,15 +346,14 @@
axlims = [axlims(1:2) plotdat.Yrange];
axis(axlims)
-plotdat.ylabel{1} = plotdat.myregion;
switch lower(plotdat.copystring)
case 'bias'
% plot a zero-bias line
- 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);
+ zeroline = line(axlims(1:2),[0 0], 'Color',[0 100 0]/255,'Parent',ax1);
+ set(zeroline,'LineWidth',2.5,'LineStyle','-')
+ plotdat.ylabel = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
otherwise
- plotdat.ylabel{2} = sprintf('%s',plotdat.copystring);
+ plotdat.ylabel = sprintf('%s',plotdat.copystring);
end
% hokey effort to decide to plot months/days vs. daynum vs.
@@ -356,44 +370,37 @@
else
xlabelstring = 'days';
end
+set(get(ax1,'Xlabel'),'String',xlabelstring, ...
+ 'Interpreter','none','FontSize',figdata.fontsize)
-% only put x axis on last/bottom plot
-if (plotdat.region == plotdat.nregions)
- xlabel(xlabelstring)
-end
+title({plotdat.myregion, plotdat.title, plotdat.subtitle}, ...
+ 'Interpreter', 'none', 'Fontsize', figdata.fontsize, 'FontWeight', 'bold')
+BottomAnnotation(plotdat.fname)
-% 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
-
% create a separate scale for the number of observations
ax2 = axes('position',get(ax1,'Position'), ...
'XAxisLocation','top', ...
- 'YAxisLocation','right',...
- 'Color','none',...
- 'XColor','b','YColor','b');
+ 'YAxisLocation','right', ...
+ 'Color','none', ...
+ 'XColor',get(ax1,'Xcolor'), ...
+ 'YColor','b', ...
+ 'FontSize',get(ax1,'FontSize'));
+
h2 = line(t,nobs_poss,'Color','b','Parent',ax2);
h3 = line(t,nobs_used,'Color','b','Parent',ax2);
set(h2,'LineStyle','none','Marker','o');
-set(h3,'LineStyle','none','Marker','+');
+set(h3,'LineStyle','none','Marker','*');
-% use same X ticks - with no labels
-set(ax2,'XTick',get(ax1,'XTick'),'XTickLabel',[])
-
+% use same X ticks
% use the same Y ticks, but find the right label values
-% Make sure that values at ticks are whole numbers.
+set(ax2,'XTick', get(ax1,'XTick'), 'XTicklabel', []);
matchingYticks(ax1,ax2);
-set(get(ax2,'Ylabel'),'String','# of obs : o=poss, +=used')
-set(get(ax1,'Ylabel'),'String',plotdat.ylabel,'Interpreter','none')
+set(get(ax1,'Ylabel'), 'String', plotdat.ylabel, ...
+ 'Interpreter','none','FontSize',figdata.fontsize)
+set(get(ax2,'Ylabel'),'String','# of obs : o=poss, \ast=used', ...
+ 'FontSize',figdata.fontsize)
-
%=====================================================================
@@ -516,6 +523,27 @@
%=====================================================================
+function figdata = setfigure()
+%%
+% figure out a page layout
+% extra space at the bottom for the date/file annotation
+% extra space at the top because the titles have multiple lines
+
+orientation = 'landscape';
+fontsize = 16;
+position = [0.10 0.15 0.8 0.7];
+linewidth = 2.0;
+
+figdata = struct('expcolors', {{'k','r','g','m','b','c','y'}}, ...
+ 'expsymbols', {{'o','s','d','p','h','s','*'}}, ...
+ 'prpolines', {{'-','--'}}, 'position', position, ...
+ 'fontsize',fontsize, 'orientation',orientation, ...
+ 'linewidth',linewidth);
+
+
+%=====================================================================
+
+
function value = local_nc_varget(fname,varname)
%% If the variable exists in the file, return the contents of the variable.
% if the variable does not exist, return empty value instead of error-ing
Modified: DART/trunk/diagnostics/matlab/plot_obs_netcdf.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_obs_netcdf.m 2014-12-24 20:54:42 UTC (rev 7330)
+++ DART/trunk/diagnostics/matlab/plot_obs_netcdf.m 2014-12-24 21:22:45 UTC (rev 7331)
@@ -5,6 +5,9 @@
% plotted on a separate figure ... color-coded to its QC value, not the
% observation value.
%
+% If you only want to plot the locations, it is easiest to simply use
+% read_obs_netcdf (not plot_obs_netcdf)
+%
%--------------------------------------------------
% EXAMPLE 1: plotting just one type of observation
%--------------------------------------------------
@@ -20,7 +23,7 @@
% bob = plot_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, maxgoodQC, verbose, twoup);
%
%--------------------------------------------------
-% EXAMPLE 2: plotting all the observation types
+% EXAMPLE 2: plotting all the observation types
%--------------------------------------------------
% fname = 'obs_epoch_001.nc';
% ObsTypeString = 'ALL';
@@ -32,6 +35,21 @@
% twoup = 1; % anything > 0 == 'true'
%
% bob = plot_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, maxgoodQC, verbose, twoup);
+%
+%--------------------------------------------------
+% EXAMPLE 3: just plotting the locations without regard to value.
+% NOTE: this uses READ_obs_netcdf
+%--------------------------------------------------
+% fname = 'obs_epoch_001.nc';
+% ObsTypeString = 'ALL';
+% region = [0 360 -90 90 -Inf Inf];
+% CopyString = 'observation';
+% QCString = 'DART quality control ';
+% verbose = 1; % anything > 0 == 'true'
+%
+% bob = read_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, verbose);
+% plot(bob.lons,bob.lats,'*')
+% continents('light');
%% 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
@@ -43,7 +61,7 @@
error('%s does not exist.',fname)
end
-if ( twoup > 0 )
+if ( twoup > 0 )
clf; orient tall
positions = [0.1, 0.55, 0.8, 0.35 ; ...
0.1, 0.10, 0.8, 0.35 ; ...
@@ -77,7 +95,7 @@
flaggedobs.qc = obsstruct.qc( inds);
end
- fprintf('Removing %d obs with a %s value greater than %f\n', ...
+ fprintf('%d obs have a %s value greater than %f\n', ...
length(inds), QCString, maxQC)
inds = find(obsstruct.qc <= maxQC);
@@ -91,6 +109,50 @@
end
+%% Separate out all the NaN values
+% If there are NaN's that make it through to the plotting phase,
+% the entire plot gets wiped out when the continents get plotted.
+% The locations with NaN values will be plotted with separate symbols.
+
+ntotal = length(obsstruct.obs);
+ngood = sum(isfinite(obsstruct.obs));
+
+if ntotal == ngood
+
+ nanobs.numflagged = 0;
+
+else
+
+ % Creating a structure with the unplottable values.
+
+ inds = find(isfinite(obsstruct.obs) == 0);
+
+ nanobs.numflagged = length(inds);
+ nanobs.string = sprintf('%d obs have unplottable values.\n', nanobs.numflagged);
+
+ if (~isempty(inds))
+ nanobs.lons = obsstruct.lons(inds);
+ nanobs.lats = obsstruct.lats(inds);
+ nanobs.Ztyp = obsstruct.Ztyp(inds);
+ nanobs.z = obsstruct.z( inds);
+ nanobs.obs = 1;
+ nanobs.qc = obsstruct.qc( inds);
+ end
+
+ % Removing the unplottable values from the plottable ones.
+ fprintf('%d obs have unplottable values.\n', nanobs.numflagged)
+
+ inds = isfinite(obsstruct.obs);
+
+ bob = obsstruct.lons(inds); obsstruct.lons = bob;
+ bob = obsstruct.lats(inds); obsstruct.lats = bob;
+ bob = obsstruct.Ztyp(inds); obsstruct.Ztyp = bob;
+ bob = obsstruct.z( inds); obsstruct.z = bob;
+ bob = obsstruct.obs( inds); obsstruct.obs = bob;
+ bob = obsstruct.qc( inds); obsstruct.qc = bob;
+
+end
+
%% Create graphic with area-weighted symbols for the good observations.
% It has happened that there have been zero good observations in a file.
@@ -103,42 +165,55 @@
pstruct.colorbarstring = obsstruct.ObsTypeString;
pstruct.region = region;
+pstruct.str1 = sprintf('%s',obsstruct.ObsTypeString);
pstruct.str3 = sprintf('%s - %s',obsstruct.timestring(1,:),obsstruct.timestring(2,:));
-if ( length(obsstruct.obs) < 1 )
+subplot('position',positions(1,:))
+
+if ( length(obsstruct.obs) < 1 )
fprintf('There are no ''good'' observations to plot\n')
- return
+ % may still want to plot the obs with bad qc values.
+ str1 = sprintf('There are no ''good'' observations to plot\n');
+ text(0.5,0.67,str1,'HorizontalAlignment','center')
+ text(0.5,0.33,nanobs.string,'HorizontalAlignment','center')
+ title( {pstruct.str1, obsstruct.CopyString, pstruct.str3}, 'Interpreter','none','FontSize',14);
else
- subplot('position',positions(1,:))
-
% choose a symbol size based on the number of obs to plot.
- if (length(obsstruct.obs) < 1000)
+ if (length(obsstruct.obs) < 1000)
pstruct.scalearray = scaleme(obsstruct.obs, 30);
else
- pstruct.scalearray = 128.0 * ones(size(obsstruct.obs));
+ pstruct.scalearray = 30.0 * ones(size(obsstruct.obs));
end
pstruct.clim = [min(obsstruct.obs) max(obsstruct.obs)];
pstruct.str2 = sprintf('%s (%d locations)',obsstruct.CopyString,length(obsstruct.obs));
% If all the observations live on the same level ... make a 2D plot.
- if ( zmin ~= zmax )
+ if ( zmin == zmax )
+ pstruct.axis = [xmin xmax ymin ymax];
+
+ plot_2D(obsstruct, pstruct);
+
+ else
+
pstruct.axis = [xmin xmax ymin ymax zmin zmax];
pstruct.str1 = sprintf('%s level (%.2f - %.2f)',obsstruct.ObsTypeString,zmin,zmax);
plot_3D(obsstruct, pstruct);
- else
+ end
- pstruct.axis = [xmin xmax ymin ymax ];
- pstruct.str1 = sprintf('%s',obsstruct.ObsTypeString);
+ % If there are unplottable values ... report that.
- plot_2D(obsstruct, pstruct);
+ if nanobs.numflagged > 0
+ subplot('position',positions(3,:))
+ axis off
+ text(0.5,0.5,nanobs.string,'HorizontalAlignment','center')
+ end
- end
end
%% Create graphic of spatial distribution of 'flagged' observations & their QC value.
@@ -169,21 +244,21 @@
if (twoup <= 0)
figure(gcf+1); clf
end
-
+
subplot('position',positions(2,:))
-
+
zmin = min(flaggedobs.z);
zmax = max(flaggedobs.z);
prej = 100.0 * length(flaggedobs.obs) / ...
(length(flaggedobs.obs) + length(obsstruct.obs));
- pstruct.scalearray = 128 * ones(size(flaggedobs.obs));
+ pstruct.scalearray = 30 * ones(size(flaggedobs.obs));
pstruct.colorbarstring = QCString;
pstruct.clim = [min(flaggedobs.qc) max(flaggedobs.qc)];
pstruct.str1 = sprintf('%s level (%.2f - %.2f)',obsstruct.ObsTypeString,zmin,zmax);
pstruct.str2 = sprintf('%s (%d ''good'', %d ''flagged'' -- %.2f %%)', obsstruct.CopyString, ...
length(obsstruct.obs), length(flaggedobs.obs), prej);
-
+
flaggedobs.obs = flaggedobs.qc; % plot QC values, not obs values
if ( zmin ~= zmax )
@@ -198,7 +273,7 @@
plot_2D(flaggedobs, pstruct);
end
-
+
subplot('position',positions(3,:))
axis off
@@ -207,7 +282,7 @@
switch lower(strtrim(QCString))
case 'dart quality control',
- qcvals = unique(flaggedobs.qc);
+ qcvals = unique(flaggedobs.qc);
qccount = zeros(size(qcvals));
s = cell(length(qcvals));
for i = 1:length(qcvals)
@@ -215,7 +290,7 @@
s{i} = sprintf('%d obs with qc == %d %s',qccount(i),qcvals(i), ...
dartqc_strings{qcvals(i)});
end
-
+
dy = 0.8*1.0/length(s);
for i = 1:length(s)
text(0.0, (i-1)*dy ,s{i})
@@ -280,7 +355,7 @@
orgholdstate = ishold;
hold on;
-if ( length(ax) > 4)
+if ( length(ax) > 4)
switch get(gca,'ZDir')
case 'reverse'
zlevel = max(ax(5:6));
@@ -292,7 +367,7 @@
fcolor = [0.7 0.7 0.7]; % light grey
myclim = get(gca,'CLim');
-[c,h] = contourf(x,y,elev,[0.0 0.0],'k-');
+[~,h] = contourf(x,y,elev,[0.0 0.0],'k-');
set(gca,'CLim',myclim)
h_patch = get(h, 'Children');
@@ -315,7 +390,7 @@
function s = scaleme(x,minsize)
% scaleme returns a uniformly scaled array the same size as the input
-% array where the maximum is 10 times the minimum
+% array where the maximum is 10 times the minimum
maxsize = 10*minsize;
minx = min(x);
maxx = max(x);
@@ -330,7 +405,7 @@
if (pstruct.clim(1) == pstruct.clim(2))
% If all the observations have the same value, setting the
- % colorbar limits is a real pain. Fundamentally, I am
+ % colorbar limits is a real pain. Fundamentally, I am
% forcing the plot symbols to be the lowest color of the
% colormap and setting the colorbar to have some more
% colors 'on top' - that are never used.
@@ -378,7 +453,7 @@
if (pstruct.clim(1) == pstruct.clim(2))
% If all the observations have the same value, setting the
- % colorbar limits is a real pain. Fundamentally, I am
+ % colorbar limits is a real pain. Fundamentally, I am
% forcing the plot symbols to be the lowest color of the
% colormap and setting the colorbar to have some more
% colors 'on top' - that are never used.
@@ -389,11 +464,9 @@
set(gca,'XGrid','on','YGrid','on')
else
-
- scatter3(obsstruct.lons, obsstruct.lats, zeros(size(obsstruct.obs)), ...
+ scatter3(obsstruct.lons, obsstruct.lats, 0.0, ...
pstruct.scalearray, obsstruct.obs, 'd', 'filled');
end
-
h1 = gca;
clim = get(h1,'CLim');
Modified: DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m 2014-12-24 20:54:42 UTC (rev 7330)
+++ DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m 2014-12-24 21:22:45 UTC (rev 7331)
@@ -7,34 +7,34 @@
%
% fname the name of the netCDF file (from obs_seq_to_netcdf)
%
-% ObsTypeString the variable of interest (from ObsTypesMetaData variable)
+% ObsTypeString the variable of interest (from ObsTypesMetaData variable)
%
% region region of interest [lonmin lonmax latmin latmax zmin zmax]
%
% CopyString1 the difference is taken 'CopyString2 - CopyString1'
% CopyString2
%
-% QCString There are multiple QC copies
+% QCString There are multiple QC copies
% maxQC The highest QC value of interest. Anything more than this
% will not be differenced. The locations will be plotted on
% a separate axis.
-% verbose logical flag ... if 'true', a table listing the possible
+% verbose logical flag ... if 'true', a table listing the possible
% observation types and observation counts is displayed.
%
-% twoup logical flag indicating that both the plot of the rejected
-% observations and the plot of the differences is
+% twoup logical flag indicating that both the plot of the rejected
+% observations and the plot of the differences is
% created on the same figure window.
%
-% The 'copies' are recorded in the netCDF 'CopyMetaData' variable -
-% the observation types are recorded in the 'ObsTypesMetaData' variable,
+% The 'copies' are recorded in the netCDF 'CopyMetaData' variable -
+% the observation types are recorded in the 'ObsTypesMetaData' variable,
% and the QC strings of interest are recorded in QCMetaData - so
% ncdump -v CopyMetaData,ObsTypesMetaData,QCMetaData obs_epoch_001.nc
% is a useful endeavor.
%
%--------------------------------------------------
-% EXAMPLE : plot the difference between the ensemble mean of the prior
+% EXAMPLE : plot the difference between the ensemble mean of the prior
% and the actual observation value - while rejecting any obs that had
-% a DART QC greater than 3 ( prior forward operator failed ... or worse)
+% a DART QC greater than 3 ( prior forward operator failed ... or worse)
%--------------------------------------------------
% fname = 'obs_epoch_001.nc';
% ObsTypeString = 'RADIOSONDE_U_WIND_COMPONENT';
@@ -59,7 +59,7 @@
error('%s does not exist.',fname)
end
-if ( twoup > 0 )
+if ( twoup > 0 )
clf; orient tall
positions = [0.1, 0.55, 0.8, 0.35 ; ...
0.1, 0.10, 0.8, 0.35 ; ...
@@ -100,8 +100,8 @@
flaggedobs.qc = obsstruct.qc( inds);
end
- fprintf('Removing %d obs with a %s value greater than %f\n', ...
- length(inds),QCString,maxQC)
+ fprintf('%d obs have a %s value greater than %f\n', ...
+ length(inds), QCString, maxQC)
inds = find(obsstruct.qc <= maxQC);
@@ -114,6 +114,50 @@
end
+%% Separate out all the NaN values
+% If there are NaN's that make it through to the plotting phase,
+% the entire plot gets wiped out when the continents get plotted.
+% The locations with NaN values will be plotted with separate symbols.
+
+ntotal = length(obsstruct.obs);
+ngood = sum(isfinite(obsstruct.obs));
+
+if ntotal == ngood
+
+ nanobs.numflagged = 0;
+
+else
+
+ % Creating a structure with the unplottable values.
+
+ inds = find(isfinite(obsstruct.obs) == 0);
+
+ nanobs.numflagged = length(inds);
+ nanobs.string = sprintf('%d obs have unplottable values.\n', nanobs.numflagged);
+
+ if (~isempty(inds))
+ nanobs.lons = obsstruct.lons(inds);
+ nanobs.lats = obsstruct.lats(inds);
+ nanobs.Ztyp = obsstruct.Ztyp(inds);
+ nanobs.z = obsstruct.z( inds);
+ nanobs.obs = 1;
+ nanobs.qc = obsstruct.qc( inds);
+ end
+
+ % Removing the unplottable values from the plottable ones.
+ fprintf('%d obs have unplottable values.\n', nanobs.numflagged)
+
+ inds = isfinite(obsstruct.obs);
+
+ bob = obsstruct.lons(inds); obsstruct.lons = bob;
+ bob = obsstruct.lats(inds); obsstruct.lats = bob;
+ bob = obsstruct.Ztyp(inds); obsstruct.Ztyp = bob;
+ bob = obsstruct.z( inds); obsstruct.z = bob;
+ bob = obsstruct.obs( inds); obsstruct.obs = bob;
+ bob = obsstruct.qc( inds); obsstruct.qc = bob;
+
+end
+
%% Create graphic with area-weighted symbols for the good observations.
% It has happened that there have been zero good observations in a file.
@@ -126,41 +170,55 @@
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list