[Dart-dev] [4341] DART/trunk/diagnostics/matlab: All functions that create vertical profiles now use the same methodology
nancy at ucar.edu
nancy at ucar.edu
Mon Apr 5 15:51:11 MDT 2010
Revision: 4341
Author: thoar
Date: 2010-04-05 15:51:11 -0600 (Mon, 05 Apr 2010)
Log Message:
-----------
All functions that create vertical profiles now use the same methodology
for creating the axes. All functions also annotate the files used in
the plotting as teeny tiny text at the bottom of the figure.
The range of time that was used in the averaging is also annotated
in the xlabel - and now properly reflects the 'time_to_skip' attribute
that was used when obs_diag was run.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
DART/trunk/diagnostics/matlab/plot_profile.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.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 2010-04-02 22:43:12 UTC (rev 4340)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2010-04-05 21:51:11 UTC (rev 4341)
@@ -35,50 +35,54 @@
% Harvest plotting info/metadata from netcdf file.
-plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.fname = fname;
+plotdat.copystring = copystring;
-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.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
-plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
-plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
-plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
-plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
-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.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.bincenters = nc_varget(fname, 'time');
-plotdat.binedges = nc_varget(fname, 'time_bounds');
-plotdat.mlevel = nc_varget(fname, 'mlevel');
-plotdat.plevel = nc_varget(fname, 'plevel');
-plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
-plotdat.hlevel = nc_varget(fname, 'hlevel');
-plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
+plotdat.bincenters = nc_varget(fname, 'time');
+plotdat.binedges = nc_varget(fname, 'time_bounds');
+plotdat.mlevel = nc_varget(fname, 'mlevel');
+plotdat.plevel = nc_varget(fname, 'plevel');
+plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
+plotdat.hlevel = nc_varget(fname, 'hlevel');
+plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
-diminfo = nc_getdiminfo(fname,'region');
-plotdat.nregions = diminfo.Length;
-region_names = nc_varget(fname,'region_names');
-plotdat.region_names = deblank(region_names);
+diminfo = nc_getdiminfo(fname,'region');
+plotdat.nregions = diminfo.Length;
+region_names = nc_varget(fname,'region_names');
+plotdat.region_names = deblank(region_names);
% Coordinate between time types and dates
-calendar = nc_attget(fname,'time','calendar');
-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;
+calendar = nc_attget(fname,'time','calendar');
+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));
+timefloats = zeros(size(time_to_skip)); % stupid int32 type conversion
+timefloats(:) = time_to_skip(:);
+skip_seconds = timefloats(4)*3600 + timefloats(5)*60 + timefloats(6);
+iskip = timefloats(3) + skip_seconds/86400.0;
-plotdat.bincenters = plotdat.bincenters + timeorigin;
-plotdat.binedges = plotdat.binedges + timeorigin;
-plotdat.Nbins = length(plotdat.bincenters);
-plotdat.toff = plotdat.bincenters(1) + iskip;
+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(:))));
+
% set up a structure with all static plotting components
-plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copystring);
+plotdat.xlabel = {sprintf('bias (%s) and %s',plotdat.biasconv,copystring), plotdat.timespan};
plotdat.linewidth = 2.0;
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
@@ -206,7 +210,7 @@
if (plotdat.nregions > 2)
CenterAnnotation(plotdat.myvarname)
end
- % BottomAnnotation(ges)
+ BottomAnnotation(fname)
% create a postscript file
print(gcf,'-dpsc','-append',psfname);
@@ -274,58 +278,48 @@
end
Stripes(plotdat.Xrange, plotdat.level_edges);
+ set(ax1,'YDir', plotdat.YDir,'YTick',plotdat.level)
+ ylabel(plotdat.level_units)
+
+ %% 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:');
- hold off;
set(h1,'LineWidth',plotdat.linewidth);
- h = legend(h1, str_bias_pr, str_bias_po, str_other_pr, str_other_po, ...
- 'Location','East');
- legend(h,'boxoff')
- set(h,'Interpreter','none')
- axlims = [plotdat.Xrange plotdat.Yrange];
- axis(axlims)
+ biasline = line([0 0],plotdat.Yrange,'Color','k','Parent',ax1);
+ set(biasline,'LineWidth',2.0,'LineStyle','-.')
- set(gca,'YDir', plotdat.YDir)
- hold on; plot([0 0],plotdat.Yrange,'k-')
+ hold off;
+ h = legend(h1, str_bias_pr, str_bias_po, str_other_pr, str_other_po);
+ set(h,'Interpreter','none','Box','off')
- set(gca,'YTick',plotdat.level,'Ylim',plotdat.Yrange)
- ylabel(plotdat.level_units)
-
- % use same X,Y limits for all plots in this region
- nXticks = length(get(ax1,'XTick'));
- xlimits = plotdat.Xrange;
- xinc = (xlimits(2)-xlimits(1))/(nXticks-1);
- xticks = xlimits(1):xinc:xlimits(2);
- set(ax1,'XTick',xticks,'Xlim',xlimits)
-
- % create a separate scale for the number of observations
+ %% Create another axes to use for plotting the observation counts
+
ax2 = axes('position',get(ax1,'Position'), ...
'XAxisLocation','top', ...
'YAxisLocation','right',...
'Color','none',...
'XColor','b','YColor','b',...
- 'YLim',plotdat.Yrange, ...
- 'YDir',plotdat.YDir);
+ 'YLim',get(ax1,'YLim'), ...
+ 'YDir',get(ax1,'YDir'));
+
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','+');
- % use same number of X ticks and the same Y ticks
-
- xlimits = get(ax2,'XLim');
- xinc = (xlimits(2)-xlimits(1))/(nXticks-1);
- xticks = xlimits(1):xinc:xlimits(2);
- nicexticks = round(10*xticks')/10;
- set(ax2,'YTick',get(ax1,'YTick'),'YTicklabel',[], ...
- 'XTick', xticks,'XTicklabel',num2str(nicexticks))
-
+ % use same Y ticks
+ set(ax2,'YTick', get(ax1,'YTick'), ...
+ 'YTicklabel',get(ax1,'YTicklabel'));
+
+ % use the same X ticks, but find the right label values
+ [xticks, newticklabels] = matchingXticks(ax1,ax2);
+ set(ax2,'XTick', xticks, 'XTicklabel', newticklabels)
+
set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
- set(ax1,'Position',get(ax2,'Position'))
- grid
if (plotdat.nregions <=2 )
title({plotdat.myvarname, plotdat.myregion}, ...
@@ -350,16 +344,18 @@
function BottomAnnotation(main)
-% annotates the directory containing the data being plotted
-subplot('position',[0.48 0.01 0.04 0.04])
+% annotates the filename containing the data being plotted
+subplot('position',[0.10 0.01 0.8 0.04])
axis off
-bob = which(main);
-[pathstr,name,ext,versn] = fileparts(bob);
-h = text(0.0,0.5,pathstr);
-set(h,'HorizontalAlignment','center', ...
- 'VerticalAlignment','middle',...
- 'Interpreter','none',...
- 'FontSize',8)
+if ( main(1) == '/' ) % must be a 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');
@@ -370,6 +366,8 @@
end
j = 0;
+basenames = struct([]);
+basedims = struct([]);
for i = 1:length(x.allvarnames)
indx = findstr('time',x.allvardims{i});
@@ -382,9 +380,11 @@
end
[b,i,j] = unique(basenames);
+y = struct([]);
+ydims = struct([]);
for j = 1:length(i)
- disp(sprintf('%2d is %s',j,basenames{j}))
+ fprintf('%2d is %s\n',j,basenames{j})
y{j} = basenames{j};
ydims{j} = basedims{j};
end
@@ -472,17 +472,65 @@
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
-function Stripes(x,edges)
-% EraseMode: [ {normal} | background | xor | none ]
-xc = [ x(1) x(2) x(2) x(1) x(1) ];
+function h = Stripes(x,edges)
+%% 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.
+
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.
+axis(axlims)
+
+xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ];
+
for i = 1:2:(length(edges)-1)
yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
- fill(xc,yc,[0.8 0.8 0.8],'EraseMode','background','EdgeColor','none');
+ hf = fill(xc,yc,[0.8 0.8 0.8],'EdgeColor','none');
+ set(hf,'FaceAlpha',0.3,'EdgeAlpha',0.3)
+ set(hf,'AlphaDataMapping','none','FaceVertexAlphaData',0.3)
end
+set(gca,'XGrid','on')
hold off;
+set(h,'Visible','off')
+
+
+
+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);
+
Modified: DART/trunk/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_profile.m 2010-04-02 22:43:12 UTC (rev 4340)
+++ DART/trunk/diagnostics/matlab/plot_profile.m 2010-04-05 21:51:11 UTC (rev 4341)
@@ -19,7 +19,7 @@
% plotdat = plot_profile(fname,copystring);
-%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+%% DART software - Copyright � 2004 - 2010 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
%
@@ -36,50 +36,54 @@
% Harvest plotting info/metadata from netcdf file.
-plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.fname = fname;
+plotdat.copystring = copystring;
-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.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
-plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
-plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
-plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
-plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
-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.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.bincenters = nc_varget(fname, 'time');
-plotdat.binedges = nc_varget(fname, 'time_bounds');
-plotdat.mlevel = nc_varget(fname, 'mlevel');
-plotdat.plevel = nc_varget(fname, 'plevel');
-plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
-plotdat.hlevel = nc_varget(fname, 'hlevel');
-plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
+plotdat.bincenters = nc_varget(fname, 'time');
+plotdat.binedges = nc_varget(fname, 'time_bounds');
+plotdat.mlevel = nc_varget(fname, 'mlevel');
+plotdat.plevel = nc_varget(fname, 'plevel');
+plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
+plotdat.hlevel = nc_varget(fname, 'hlevel');
+plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
-diminfo = nc_getdiminfo(fname,'region');
-plotdat.nregions = diminfo.Length;
-region_names = nc_varget(fname,'region_names');
-plotdat.region_names = deblank(region_names);
+diminfo = nc_getdiminfo(fname,'region');
+plotdat.nregions = diminfo.Length;
+region_names = nc_varget(fname,'region_names');
+plotdat.region_names = deblank(region_names);
% Coordinate between time types and dates
-calendar = nc_attget(fname,'time','calendar');
-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;
+calendar = nc_attget(fname,'time','calendar');
+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));
+timefloats = zeros(size(time_to_skip)); % stupid int32 type conversion
+timefloats(:) = time_to_skip(:);
+skip_seconds = timefloats(4)*3600 + timefloats(5)*60 + timefloats(6);
+iskip = timefloats(3) + skip_seconds/86400.0;
-plotdat.bincenters = plotdat.bincenters + timeorigin;
-plotdat.binedges = plotdat.binedges + timeorigin;
-plotdat.Nbins = length(plotdat.bincenters);
-plotdat.toff = plotdat.bincenters(1) + iskip;
+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(:))));
+
% set up a structure with all static plotting components
-plotdat.xlabel = sprintf('%s',copystring);
+plotdat.xlabel = {sprintf('%s',copystring), plotdat.timespan};
plotdat.linewidth = 2.0;
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
@@ -203,7 +207,7 @@
if (plotdat.nregions > 2)
CenterAnnotation(plotdat.myvarname)
end
- % BottomAnnotation(ges)
+ BottomAnnotation(fname)
% create a postscript file
print(gcf,'-dpsc','-append',psfname);
@@ -262,61 +266,54 @@
end
Stripes(plotdat.Xrange, plotdat.level_edges);
+ set(ax1,'YDir', plotdat.YDir,'YTick',plotdat.level)
+ ylabel(plotdat.level_units)
+
+ %% draw the result of the experiment
+
hold on;
h1 = plot(CG,plotdat.level,'k+-',CA,plotdat.level,'k+:');
- hold off;
+
set(h1,'LineWidth',plotdat.linewidth);
- h = legend(h1,str_other_pr, str_other_po, 'Location', 'East');
+ h = legend(h1, str_other_pr, str_other_po);
legend(h,'boxoff')
set(h,'Interpreter','none')
- axlims = [plotdat.Xrange plotdat.Yrange];
- axis(axlims)
-
- set(gca,'YDir', plotdat.YDir)
switch plotdat.copystring
case 'bias'
- hold on; plot([0 0], plotdat.Yrange,'k-')
- plotdat.xlabel = sprintf('bias (%s)',plotdat.biasconv);
+ biasline = line([0 0],plotdat.Yrange,'Color','k','Parent',ax1);
+ set(biasline,'LineWidth',2.0,'LineStyle','-.')
+ plotdat.xlabel = {sprintf('bias (%s)',plotdat.biasconv), plotdat.timespan};
otherwise
end
- set(gca,'YTick',plotdat.level,'Ylim',plotdat.Yrange)
- ylabel(plotdat.level_units)
-
- % use same X,Y limits for all plots in this region
- nXticks = length(get(ax1,'XTick'));
- xlimits = plotdat.Xrange;
- xinc = (xlimits(2)-xlimits(1))/(nXticks-1);
- xticks = xlimits(1):xinc:xlimits(2);
- set(ax1,'XTick',xticks,'Xlim',xlimits)
-
- % create a separate scale for the number of observations
+ hold off;
+
+ %% Create another axes to use for plotting the observation counts
+
ax2 = axes('position',get(ax1,'Position'), ...
'XAxisLocation','top', ...
'YAxisLocation','right',...
'Color','none',...
'XColor','b','YColor','b',...
- 'YLim',plotdat.Yrange, ...
- 'YDir',plotdat.YDir);
+ 'YLim',get(ax1,'YLim'), ...
+ 'YDir',get(ax1,'YDir'));
+
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','+');
- % use same number of X ticks and the same Y ticks
-
- xlimits = get(ax2,'XLim');
- xinc = (xlimits(2)-xlimits(1))/(nXticks-1);
- xticks = xlimits(1):xinc:xlimits(2);
- nicexticks = round(10*xticks')/10;
- set(ax2,'YTick',get(ax1,'YTick'),'YTicklabel',[], ...
- 'XTick', xticks,'XTicklabel',num2str(nicexticks))
-
+ % use same Y ticks
+ set(ax2,'YTick', get(ax1,'YTick'), ...
+ 'YTicklabel',get(ax1,'YTicklabel'));
+
+ % use the same X ticks, but find the right label values
+ [xticks, newticklabels] = matchingXticks(ax1,ax2);
+ set(ax2,'XTick', xticks, 'XTicklabel', newticklabels)
+
set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
- set(ax1,'Position',get(ax2,'Position'))
- grid
if (plotdat.nregions <=2 )
title({plotdat.myvarname, plotdat.myregion}, ...
@@ -341,16 +338,18 @@
function BottomAnnotation(main)
-% annotates the directory containing the data being plotted
-subplot('position',[0.48 0.01 0.04 0.04])
+% annotates the filename containing the data being plotted
+subplot('position',[0.10 0.01 0.8 0.04])
axis off
-bob = which(main);
-[pathstr,name,ext,versn] = fileparts(bob);
-h = text(0.0,0.5,pathstr);
-set(h,'HorizontalAlignment','center', ...
- 'VerticalAlignment','middle',...
- 'Interpreter','none',...
- 'FontSize',8)
+if ( main(1) == '/' ) % must be a 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');
@@ -361,6 +360,8 @@
end
j = 0;
+basenames = struct([]);
+basedims = struct([]);
for i = 1:length(x.allvarnames)
indx = findstr('time',x.allvardims{i});
@@ -373,9 +374,11 @@
end
[b,i,j] = unique(basenames);
+y = struct([]);
+ydims = struct([]);
for j = 1:length(i)
- disp(sprintf('%2d is %s',j,basenames{j}))
+ fprintf('%2d is %s\n',j,basenames{j})
y{j} = basenames{j};
ydims{j} = basedims{j};
end
@@ -463,17 +466,69 @@
Yrange = [ymin ymax];
+ % Make sure a zero bias is visible on plot
+ switch y.copystring
+ case {'bias'}
+ if ymax < 0
+ Yrange = [ ymin 0.0 ];
+ elseif ymin > 0
+ Yrange = [ 0.0 ymax ];
+ end
+ end
+
x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
end
-function Stripes(x,edges)
-% EraseMode: [ {normal} | background | xor | none ]
-xc = [ x(1) x(2) x(2) x(1) x(1) ];
+function h = Stripes(x,edges)
+%% 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.
+
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.
+axis(axlims)
+
+xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ];
+
for i = 1:2:(length(edges)-1)
yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
- fill(xc,yc,[0.8 0.8 0.8],'EraseMode','background','EdgeColor','none');
+ hf = fill(xc,yc,[0.8 0.8 0.8],'EdgeColor','none');
+ set(hf,'FaceAlpha',0.3,'EdgeAlpha',0.3)
+ set(hf,'AlphaDataMapping','none','FaceVertexAlphaData',0.3)
end
+set(gca,'XGrid','on')
hold off;
+set(h,'Visible','off')
+
+
+
+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);
+
Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m 2010-04-02 22:43:12 UTC (rev 4340)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m 2010-04-05 21:51:11 UTC (rev 4341)
@@ -18,7 +18,8 @@
% copystring = 'totalspread'; % 'copy' string == quantity of interest
% plotdat = plot_rmse_xxx_profile(fname,copystring);
-%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+
+%% DART software - Copyright � 2004 - 2010 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
%
@@ -35,50 +36,54 @@
% Harvest plotting info/metadata from netcdf file.
-plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.fname = fname;
+plotdat.copystring = copystring;
-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.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
-plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
-plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
-plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
-plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
-plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
-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.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.bincenters = nc_varget(fname, 'time');
-plotdat.binedges = nc_varget(fname, 'time_bounds');
-plotdat.mlevel = nc_varget(fname, 'mlevel');
-plotdat.plevel = nc_varget(fname, 'plevel');
-plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
-plotdat.hlevel = nc_varget(fname, 'hlevel');
-plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
+plotdat.bincenters = nc_varget(fname, 'time');
+plotdat.binedges = nc_varget(fname, 'time_bounds');
+plotdat.mlevel = nc_varget(fname, 'mlevel');
+plotdat.plevel = nc_varget(fname, 'plevel');
+plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
+plotdat.hlevel = nc_varget(fname, 'hlevel');
+plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
-diminfo = nc_getdiminfo(fname,'region');
-plotdat.nregions = diminfo.Length;
-region_names = nc_varget(fname,'region_names');
-plotdat.region_names = deblank(region_names);
+diminfo = nc_getdiminfo(fname,'region');
+plotdat.nregions = diminfo.Length;
+region_names = nc_varget(fname,'region_names');
+plotdat.region_names = deblank(region_names);
% Coordinate between time types and dates
-calendar = nc_attget(fname,'time','calendar');
-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;
+calendar = nc_attget(fname,'time','calendar');
+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));
+timefloats = zeros(size(time_to_skip)); % stupid int32 type conversion
+timefloats(:) = time_to_skip(:);
+skip_seconds = timefloats(4)*3600 + timefloats(5)*60 + timefloats(6);
+iskip = timefloats(3) + skip_seconds/86400.0;
-plotdat.bincenters = plotdat.bincenters + timeorigin;
-plotdat.binedges = plotdat.binedges + timeorigin;
-plotdat.Nbins = length(plotdat.bincenters);
-plotdat.toff = plotdat.bincenters(1) + iskip;
+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(:))));
+
% set up a structure with all static plotting components
-plotdat.xlabel = sprintf('rmse and %s',copystring);
+plotdat.xlabel = {sprintf('rmse and %s',copystring), plotdat.timespan};
plotdat.linewidth = 2.0;
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
@@ -206,7 +211,7 @@
if (plotdat.nregions > 2)
CenterAnnotation(plotdat.myvarname)
end
- % BottomAnnotation(ges)
+ BottomAnnotation(fname)
% create a postscript file
print(gcf,'-dpsc','-append',psfname);
@@ -274,58 +279,45 @@
end
Stripes(plotdat.Xrange, plotdat.level_edges);
+ set(ax1,'YDir', plotdat.YDir,'YTick',plotdat.level)
+ ylabel(plotdat.level_units)
+
+ %% 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:');
+ plot([0 0],plotdat.Yrange,'k-')
hold off;
set(h1,'LineWidth',plotdat.linewidth);
- h = legend(h1, str_rmse_pr, str_rmse_po, str_other_pr, str_other_po, ...
- 'Location','East');
- legend(h,'boxoff')
- set(h,'Interpreter','none')
+ h = legend(h1, str_rmse_pr, str_rmse_po, str_other_pr, str_other_po);
+ set(h,'Interpreter','none','Box','off')
- axlims = [plotdat.Xrange plotdat.Yrange];
- axis(axlims)
+ %% Create another axes to use for plotting the observation counts
- set(gca,'YDir', plotdat.YDir)
- hold on; plot([0 0],plotdat.Yrange,'k-')
-
- set(gca,'YTick',plotdat.level,'Ylim',plotdat.Yrange)
- ylabel(plotdat.level_units)
-
- % use same X,Y limits for all plots in this region
- nXticks = length(get(ax1,'XTick'));
- xlimits = plotdat.Xrange;
- xinc = (xlimits(2)-xlimits(1))/(nXticks-1);
- xticks = xlimits(1):xinc:xlimits(2);
- set(ax1,'XTick',xticks,'Xlim',xlimits)
-
- % 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',...
- 'YLim',plotdat.Yrange, ...
- 'YDir',plotdat.YDir);
+ 'YLim',get(ax1,'YLim'), ...
+ 'YDir',get(ax1,'YDir'));
+
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','+');
- % use same number of X ticks and the same Y ticks
-
- xlimits = get(ax2,'XLim');
- xinc = (xlimits(2)-xlimits(1))/(nXticks-1);
- xticks = xlimits(1):xinc:xlimits(2);
- nicexticks = round(10*xticks')/10;
- set(ax2,'YTick',get(ax1,'YTick'),'YTicklabel',[], ...
- 'XTick', xticks,'XTicklabel',num2str(nicexticks))
-
+ % use same Y ticks
+ set(ax2,'YTick', get(ax1,'YTick'), ...
+ 'YTicklabel',get(ax1,'YTicklabel'));
+
+ % use the same X ticks, but find the right label values
+ [xticks, newticklabels] = matchingXticks(ax1,ax2);
+ set(ax2,'XTick', xticks, 'XTicklabel', newticklabels)
+
set(get(ax2,'Xlabel'),'String','# of obs (o=poss, +=used)')
set(get(ax1,'Xlabel'),'String',plotdat.xlabel,'Interpreter','none')
- set(ax1,'Position',get(ax2,'Position'))
- grid
if (plotdat.nregions <=2 )
title({plotdat.myvarname, plotdat.myregion}, ...
@@ -350,16 +342,18 @@
function BottomAnnotation(main)
-% annotates the directory containing the data being plotted
-subplot('position',[0.48 0.01 0.04 0.04])
+% annotates the filename containing the data being plotted
+subplot('position',[0.10 0.01 0.8 0.04])
axis off
-bob = which(main);
-[pathstr,name,ext,versn] = fileparts(bob);
-h = text(0.0,0.5,pathstr);
-set(h,'HorizontalAlignment','center', ...
- 'VerticalAlignment','middle',...
- 'Interpreter','none',...
- 'FontSize',8)
+if ( main(1) == '/' ) % must be a 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');
@@ -370,6 +364,8 @@
end
j = 0;
+basenames = struct([]);
+basedims = struct([]);
for i = 1:length(x.allvarnames)
indx = findstr('time',x.allvardims{i});
@@ -382,9 +378,11 @@
end
[b,i,j] = unique(basenames);
+y = struct([]);
+ydims = struct([]);
for j = 1:length(i)
- disp(sprintf('%2d is %s',j,basenames{j}))
+ fprintf('%2d is %s\n',j,basenames{j})
y{j} = basenames{j};
ydims{j} = basedims{j};
end
@@ -475,14 +473,56 @@
x = sort([min([Yrange(1) 0.0]) Yrange(2)] ,'ascend');
end
-function Stripes(x,edges)
-% EraseMode: [ {normal} | background | xor | none ]
-xc = [ x(1) x(2) x(2) x(1) x(1) ];
+function h = Stripes(x,edges)
+%% 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.
+
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.
+axis(axlims)
+
+xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ];
+
for i = 1:2:(length(edges)-1)
yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ];
- fill(xc,yc,[0.8 0.8 0.8],'EraseMode','background','EdgeColor','none');
+ hf = fill(xc,yc,[0.8 0.8 0.8],'EdgeColor','none');
+ set(hf,'FaceAlpha',0.3,'EdgeAlpha',0.3)
+ set(hf,'AlphaDataMapping','none','FaceVertexAlphaData',0.3)
end
+set(gca,'XGrid','on')
hold off;
+set(h,'Visible','off')
+
+
+
+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);
+
Modified: DART/trunk/diagnostics/matlab/two_experiments_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/two_experiments_profile.m 2010-04-02 22:43:12 UTC (rev 4340)
+++ DART/trunk/diagnostics/matlab/two_experiments_profile.m 2010-04-05 21:51:11 UTC (rev 4341)
@@ -7,9 +7,9 @@
% files = {'/ptmp/nancy/CSL/Base5/032-061s0_def_reg/obs_diag_output.nc',
% '/ptmp/thoar/GPS+AIRS/Sep_032-061/obs_diag_output.nc'};
% titles = {'Base5', 'GPS+AIRS'};
-% files = {'/fs/image/home/hliu/DART/models/wrf/work/ernesto/gpsonly/obs_diag_output.nc',
-% '/fs/image/home/hliu/DART/models/wrf/work/ernesto/ctl/obs_diag_output.nc'};
-% titles = {'GPS only', 'Control'};
+% files = {'/fs/image/home/hliu/DART/models/wrf/work/ernesto/ctl/obs_diag_output.nc',
+% '/fs/image/home/hliu/DART/models/wrf/work/ernesto/airsQ/obs_diag_output.nc'};
+% titles = {'Control', 'airsQ'};
% varnames = {'RADIOSONDE_U_WIND_COMPONENT', 'RADIOSONDE_TEMPERATURE'};
% qtty = 'rmse'; % rmse, spread, totalspread, bias, etc.
% prpo = 'analysis'; % [analy, analysis, posterior ] == posterior
@@ -44,9 +44,10 @@
% set up all the stuff that is common.
commondata = check_compatibility(files, varnames, qtty);
+% DEBUG - next line is for testing only
+% commondata.nregions = 2;
figuredata = setfigure(commondata);
-
%%--------------------------------------------------------------------
% Set some static data
%---------------------------------------------------------------------
@@ -94,6 +95,8 @@
myplot( plotobj, Nrange, Drange, Yrange, figuredata);
end % of loop around regions
+
+ BottomAnnotation(files)
if ( ivar ~= nvars )
disp('Pausing, hit any key to continue ...')
@@ -216,31 +219,32 @@
calendar = nc_attget(fname,'time','calendar');
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;
+timefloats = zeros(size(time_to_skip)); % stupid int32 type conversion
+timefloats(:) = time_to_skip(:);
+skip_seconds = timefloats(4)*3600 + timefloats(5)*60 + timefloats(6);
+iskip = timefloats(3) + skip_seconds/86400;
plotdat.bincenters = plotdat.bincenters + timeorigin;
plotdat.binedges = plotdat.binedges + timeorigin;
plotdat.Nbins = length(plotdat.bincenters);
plotdat.toff = plotdat.bincenters(1) + iskip;
-plotdat.timespan = sprintf('%s through %s', ...
- datestr(min(plotdat.binedges(:))), ...
+plotdat.timespan = sprintf('%s through %s', datestr(plotdat.toff), ...
datestr(max(plotdat.binedges(:))));
% Get the right indices for the intended variable, regardless of the storage order
-plotdat.copyindex = get_copy_index(fname, copystring);
-plotdat.priorvar = sprintf('%s_VPguess',plotdat.varname);
-plotdat.postevar = sprintf('%s_VPanaly',plotdat.varname);
+plotdat.copyindex = get_copy_index(fname, copystring);
+plotdat.priorvar = sprintf('%s_VPguess',plotdat.varname);
+plotdat.postevar = sprintf('%s_VPanaly',plotdat.varname);
-myinfo.diagn_file = fname;
-myinfo.copyindex = plotdat.copyindex;
-myinfo.regionindex = plotdat.region;
-[start, count] = GetNCindices(myinfo,'diagn',plotdat.priorvar);
-plotdat.prior = nc_varget(fname, plotdat.priorvar, start, count);
-[start, count] = GetNCindices(myinfo,'diagn',plotdat.postevar);
-plotdat.poste = nc_varget(fname, plotdat.postevar, start, count);
+myinfo.diagn_file = fname;
+myinfo.copyindex = plotdat.copyindex;
+myinfo.regionindex = plotdat.region;
+[start, count] = GetNCindices(myinfo,'diagn',plotdat.priorvar);
+plotdat.prior = nc_varget(fname, plotdat.priorvar, start, count);
+[start, count] = GetNCindices(myinfo,'diagn',plotdat.postevar);
+plotdat.poste = nc_varget(fname, plotdat.postevar, start, count);
% Now that we know the variable ... get the appropriate vertical information
@@ -310,7 +314,13 @@
function h = Stripes(x,edges)
-% EraseMode: [ {normal} | background | xor | none ]
+%% 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.
hold on;
@@ -319,9 +329,9 @@
% tick labels ... KEEP THEM. Later, make the dots invisible.
h = plot([min(x) max(x)],[min(edges) max(edges)]);
-axlims = axis;
+axlims = axis;
axlims(4) = max(edges);
-axlims(3) = -100;
+axlims(3) = -100; % This gives extra space for the legend.
axis(axlims)
xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ];
@@ -401,7 +411,7 @@
for i = 1:Nexp
h2 = line(plotobj{i}.nposs, plotobj{i}.levels,'Color',figdata.expcolors{i},'Parent',ax2);
h3 = line(plotobj{i}.nused, plotobj{i}.levels,'Color',figdata.expcolors{i},'Parent',ax2);
- set(h2,'LineStyle','none','Marker','o');
+ set(h2,'LineStyle','none','Marker',figdata.expsymbols{i});
set(h3,'LineStyle','none','Marker','+');
end
@@ -511,6 +521,30 @@
+function BottomAnnotation(filenames)
+% annotates the filename containing the data being plotted
+nfiles = length(filenames);
+dy = 1.0/(nfiles+1); % vertical spacing for text
+
+subplot('position',[0.10 0.01 0.8 0.05])
+axis off
+
+for ifile = 1:nfiles
+ main = filenames{ifile};
+ if ( main(1) == '/' ) % must be a absolute pathname
+ string1 = sprintf('data file: %s',main);
+ else
+ mydir = pwd;
+ string1 = sprintf('data file: %s/%s',mydir,main);
+ end
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list