[Dart-dev] [6124] DART/branches/development: The obs_diag for the 1D models/ observations now outputs a netCDF file
nancy at ucar.edu
nancy at ucar.edu
Thu May 9 15:28:47 MDT 2013
Revision: 6124
Author: thoar
Date: 2013-05-09 15:28:46 -0600 (Thu, 09 May 2013)
Log Message:
-----------
The obs_diag for the 1D models/observations now outputs a netCDF file
and can use the same plot_evolution.m and plot_rmse_xxx_evolution.m
scripts as the 3D observations.
This renders fit_ens_mean_time.m
fit_ens_spread_time.m
obs_num_time.m
fit_mean_spread_time.m obsolete. Indeed, they do not work anymore.
This is not so bad, because they have not worked in years!
Tutorial section 18 must be updated to reflect this.
Modified Paths:
--------------
DART/branches/development/diagnostics/matlab/plot_evolution.m
DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m
DART/branches/development/diagnostics/oned/obs_diag.f90
Added Paths:
-----------
DART/branches/development/matlab/nc_var_exists.m
-------------- next part --------------
Modified: DART/branches/development/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_evolution.m 2013-05-09 16:38:04 UTC (rev 6123)
+++ DART/branches/development/diagnostics/matlab/plot_evolution.m 2013-05-09 21:28:46 UTC (rev 6124)
@@ -10,18 +10,18 @@
%
% 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)
% obstypestring : 'observation type' string. Optional.
-% Must match something in the netcdf
+% Must match something in the netcdf
% file 'ObservationTypes' variable.
% (ncdump -v ObservationTypes obs_diag_output.nc)
% If specified, only this observation type will be plotted.
% If not specified, all observation types incluced in the netCDF file
% will be plotted.
%
-% OUTPUT: two files will result for each observation type plotted. One is a
+% OUTPUT: two files will result for each observation type plotted. One is a
% postscript file containing a page for each level - all regions.
% The other file is a simple text file containing summary information
% about how many observations were assimilated, how many were available, etc.
@@ -61,7 +61,6 @@
error('wrong number of arguments ... ')
end
-
if (exist(fname,'file') ~= 2)
error('file/fname <%s> does not exist',fname)
end
@@ -72,29 +71,31 @@
plotdat.copystring = copystring;
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.mlevel = local_nc_varget(fname,'mlevel');
+plotdat.plevel = local_nc_varget(fname,'plevel');
+plotdat.plevel_edges = local_nc_varget(fname,'plevel_edges');
+plotdat.hlevel = local_nc_varget(fname,'hlevel');
+plotdat.hlevel_edges = local_nc_varget(fname,'hlevel_edges');
plotdat.ncopies = length(nc_varget(fname,'copy'));
plotdat.nregions = length(nc_varget(fname,'region'));
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
-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');
+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.rat_cri = local_nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = local_nc_attget(fname, nc_global, 'input_qc_threshold');
+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
@@ -124,7 +125,7 @@
end
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
plotdat.NQC4index = get_copy_index(fname,'N_DARTqc_4');
@@ -137,24 +138,24 @@
%----------------------------------------------------------------------
for ivar = 1:plotdat.nvars
-
+
% create the variable names of interest.
-
- plotdat.myvarname = plotdat.varnames{ivar};
+
+ plotdat.myvarname = plotdat.varnames{ivar};
plotdat.guessvar = sprintf('%s_guess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_analy',plotdat.varnames{ivar});
% 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);
- disp(sprintf('Removing %s from the current directory.',psfname))
+ fprintf('Removing %s from the current directory.\n',psfname)
system(sprintf('rm %s',psfname));
- % remove any existing log file -
-
+ % remove any existing log file -
+
lgfname = sprintf('%s_%s_obscount.txt',plotdat.varnames{ivar},plotdat.copystring);
- disp(sprintf('Removing %s from the current directory.',lgfname))
+ fprintf('Removing %s from the current directory.\n',lgfname)
system(sprintf('rm %s',lgfname));
logfid = fopen(lgfname,'wt');
fprintf(logfid,'%s\n',lgfname);
@@ -164,10 +165,13 @@
guessdims = nc_var_dims(fname, plotdat.guessvar);
analydims = nc_var_dims(fname, plotdat.analyvar);
- if ( findstr('surface',guessdims{3}) > 0 )
+ if ( dimensionality == 1 ) % observations on a unit circle, no level
+ plotdat.level = 1;
+ plotdat.level_units = [];
+ elseif ( strfind(guessdims{3},'surface') > 0 )
plotdat.level = 1;
plotdat.level_units = 'surface';
- elseif ( findstr('undef',guessdims{3}) > 0 )
+ elseif ( strfind(guessdims{3},'undef') > 0 )
plotdat.level = 1;
plotdat.level_units = 'undefined';
else
@@ -176,57 +180,57 @@
end
plotdat.nlevels = length(plotdat.level);
- % Here is the tricky part. Singleton dimensions are auto-squeezed ...
+ % Here is the tricky part. Singleton dimensions are auto-squeezed ...
% single levels, single regions ...
- guess_raw = nc_varget(fname, plotdat.guessvar);
+ guess_raw = nc_varget(fname, plotdat.guessvar);
guess = reshape(guess_raw, plotdat.Nbins, plotdat.ncopies, ...
- plotdat.nlevels, plotdat.nregions);
+ plotdat.nlevels, plotdat.nregions);
- analy_raw = nc_varget(fname, plotdat.analyvar);
+ analy_raw = nc_varget(fname, plotdat.analyvar);
analy = reshape(analy_raw, plotdat.Nbins, plotdat.ncopies, ...
- plotdat.nlevels, plotdat.nregions);
+ plotdat.nlevels, plotdat.nregions);
% check to see if there is anything to plot
nposs = sum(guess(:,plotdat.Npossindex,:,:));
if ( sum(nposs(:)) < 1 )
- disp(sprintf('%s no obs for %s... skipping', plotdat.varnames{ivar}))
+ fprintf('%s no obs for %s... skipping\n', plotdat.varnames{ivar})
continue
end
-
+
for ilevel = 1:plotdat.nlevels
fprintf(logfid,'\nlevel %d %f %s\n',ilevel,plotdat.level(ilevel),plotdat.level_units);
plotdat.ges_Nqc4 = guess(:,plotdat.NQC4index ,ilevel,:);
plotdat.anl_Nqc4 = analy(:,plotdat.NQC4index ,ilevel,:);
fprintf(logfid,'DART QC == 4, prior/post %d %d\n',sum(plotdat.ges_Nqc4(:)), ...
- sum(plotdat.anl_Nqc4(:)));
+ sum(plotdat.anl_Nqc4(:)));
plotdat.ges_Nqc5 = guess(:,plotdat.NQC5index ,ilevel,:);
plotdat.anl_Nqc5 = analy(:,plotdat.NQC5index ,ilevel,:);
fprintf(logfid,'DART QC == 5, prior/post %d %d\n',sum(plotdat.ges_Nqc5(:)), ...
- sum(plotdat.anl_Nqc5(:)));
+ sum(plotdat.anl_Nqc5(:)));
plotdat.ges_Nqc6 = guess(:,plotdat.NQC6index ,ilevel,:);
plotdat.anl_Nqc6 = analy(:,plotdat.NQC6index ,ilevel,:);
fprintf(logfid,'DART QC == 6, prior/post %d %d\n',sum(plotdat.ges_Nqc6(:)), ...
- sum(plotdat.anl_Nqc6(:)));
+ sum(plotdat.anl_Nqc6(:)));
plotdat.ges_Nqc7 = guess(:,plotdat.NQC7index ,ilevel,:);
plotdat.anl_Nqc7 = analy(:,plotdat.NQC7index ,ilevel,:);
fprintf(logfid,'DART QC == 7, prior/post %d %d\n',sum(plotdat.ges_Nqc7(:)), ...
- sum(plotdat.anl_Nqc7(:)));
+ sum(plotdat.anl_Nqc7(:)));
plotdat.ges_Nposs = guess(:,plotdat.Npossindex, ilevel,:);
plotdat.anl_Nposs = analy(:,plotdat.Npossindex, ilevel,:);
fprintf(logfid,'# obs poss, prior/post %d %d\n',sum(plotdat.ges_Nposs(:)), ...
- sum(plotdat.anl_Nposs(:)));
+ sum(plotdat.anl_Nposs(:)));
plotdat.ges_Nused = guess(:,plotdat.Nusedindex, ilevel,:);
plotdat.anl_Nused = analy(:,plotdat.Nusedindex, ilevel,:);
fprintf(logfid,'# obs used, prior/post %d %d\n',sum(plotdat.ges_Nused(:)), ...
- sum(plotdat.anl_Nused(:)));
+ sum(plotdat.anl_Nused(:)));
plotdat.ges_copy = guess(:,plotdat.copyindex, ilevel,:);
plotdat.anl_copy = analy(:,plotdat.copyindex, ilevel,:);
@@ -237,18 +241,22 @@
if (plotdat.nregions > 2)
clf; orient tall
- else
+ else
clf; orient landscape
end
for iregion = 1:plotdat.nregions
- plotdat.region = iregion;
+ plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
- plotdat.title = sprintf('%s @ %d %s', ...
- plotdat.myvarname, ...
- plotdat.level(ilevel), ...
- plotdat.level_units);
+ if ( isempty(plotdat.level_units) )
+ plotdat.title = plotdat.myvarname;
+ else
+ plotdat.title = sprintf('%s @ %d %s', ...
+ plotdat.myvarname, ...
+ plotdat.level(ilevel), ...
+ plotdat.level_units);
+ end
myplot(plotdat);
end
@@ -269,121 +277,121 @@
function myplot(plotdat)
- % Interlace the [ges,anl] to make a sawtooth plot.
- % By this point, the middle two dimensions are singletons.
- cg = plotdat.ges_copy(:,:,:,plotdat.region);
- ca = plotdat.anl_copy(:,:,:,plotdat.region);
+% Interlace the [ges,anl] to make a sawtooth plot.
+% By this point, the middle two dimensions are singletons.
+cg = plotdat.ges_copy(:,:,:,plotdat.region);
+ca = plotdat.anl_copy(:,:,:,plotdat.region);
- g = plotdat.ges_Nposs(:,:,:,plotdat.region);
- a = plotdat.anl_Nposs(:,:,:,plotdat.region);
- nobs_poss = reshape([g a]',2*plotdat.Nbins,1);
+g = plotdat.ges_Nposs(:,:,:,plotdat.region);
+a = plotdat.anl_Nposs(:,:,:,plotdat.region);
+nobs_poss = reshape([g a]',2*plotdat.Nbins,1);
- g = plotdat.ges_Nused(:,:,:,plotdat.region);
- a = plotdat.anl_Nused(:,:,:,plotdat.region);
- nobs_used = reshape([g a]',2*plotdat.Nbins,1);
+g = plotdat.ges_Nused(:,:,:,plotdat.region);
+a = plotdat.anl_Nused(:,:,:,plotdat.region);
+nobs_used = reshape([g a]',2*plotdat.Nbins,1);
- tg = plotdat.bincenters;
- ta = plotdat.bincenters;
- t = reshape([tg ta]',2*plotdat.Nbins,1);
+tg = plotdat.bincenters;
+ta = plotdat.bincenters;
+t = reshape([tg ta]',2*plotdat.Nbins,1);
- % Determine some quantities for the legend
- nobs = sum(nobs_used);
- if ( nobs > 1 )
- mean_prior = mean(cg(isfinite(cg)));
- mean_post = mean(ca(isfinite(ca)));
- else
- mean_prior = NaN;
- mean_post = NaN;
- end
+% Determine some quantities for the legend
+nobs = sum(nobs_used);
+if ( nobs > 1 )
+ mean_prior = mean(cg(isfinite(cg)));
+ mean_post = mean(ca(isfinite(ca)));
+else
+ mean_prior = NaN;
+ mean_post = NaN;
+end
- 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);
+string_guess = sprintf('forecast: mean=%.5g', mean_prior);
+string_analy = sprintf('analysis: mean=%.5g', mean_post);
+plotdat.subtitle = sprintf('%s %s %s',plotdat.myregion, string_guess, string_analy);
- % Plot the requested quantity on the left axis.
- % The observation count will use the axis on the right.
- % We want to suppress the 'auto' feature of the axis labelling,
- % so we manually set some values that normally
- % don't need to be set.
-
- ax1 = subplot(plotdat.nregions,1,plotdat.region);
+% Plot the requested quantity on the left axis.
+% The observation count will use the axis on the right.
+% 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.
- h1 = plot(tg,cg,'k+-',ta,ca,'ro-','LineWidth',plotdat.linewidth);
- h = legend('forecast', 'analysis');
- legend(h,'boxoff','Interpreter','none')
+ax1 = subplot(plotdat.nregions,1,plotdat.region);
- axlims = axis;
- axlims = [axlims(1:2) plotdat.Yrange];
- axis(axlims)
+h1 = plot(tg,cg,'k+-',ta,ca,'ro-','LineWidth',plotdat.linewidth);
+h = legend('forecast', 'analysis');
+legend(h,'boxoff','Interpreter','none')
- plotdat.ylabel{1} = plotdat.myregion;
- switch lower(plotdat.copystring)
- case 'bias'
- % plot a zero-bias line
- h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
- set(h4,'LineWidth',1.5,'LineSTyle',':')
- plotdat.ylabel{2} = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
- otherwise
- plotdat.ylabel{2} = sprintf('%s',plotdat.copystring);
- end
-
- % hokey effort to decide to plot months/days vs. daynum vs.
- ttot = plotdat.bincenters(plotdat.Nbins) - plotdat.bincenters(1) + 1;
-
- if ((plotdat.bincenters(1) > 1000) && (ttot > 5))
- datetick('x',6,'keeplimits','keepticks');
- monstr = datestr(plotdat.bincenters(1),21);
- xlabelstring = sprintf('month/day - %s start',monstr);
- elseif (plotdat.bincenters(1) > 1000)
- datetick('x',15,'keeplimits','keepticks')
- monstr = datestr(plotdat.bincenters(1),21);
- xlabelstring = sprintf('%s start',monstr);
- else
- xlabelstring = 'days';
- end
+axlims = axis;
+axlims = [axlims(1:2) plotdat.Yrange];
+axis(axlims)
- % only put x axis on last/bottom plot
- if (plotdat.region == plotdat.nregions)
- xlabel(xlabelstring)
- end
+plotdat.ylabel{1} = plotdat.myregion;
+switch lower(plotdat.copystring)
+ case 'bias'
+ % plot a zero-bias line
+ h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
+ set(h4,'LineWidth',1.5,'LineSTyle',':')
+ plotdat.ylabel{2} = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
+ otherwise
+ plotdat.ylabel{2} = sprintf('%s',plotdat.copystring);
+end
- % more annotation ...
+% hokey effort to decide to plot months/days vs. daynum vs.
+ttot = plotdat.bincenters(plotdat.Nbins) - plotdat.bincenters(1) + 1;
- 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');
- 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','+');
+if ((plotdat.bincenters(1) > 1000) && (ttot > 5))
+ datetick('x',6,'keeplimits','keepticks');
+ monstr = datestr(plotdat.bincenters(1),21);
+ xlabelstring = sprintf('month/day - %s start',monstr);
+elseif (plotdat.bincenters(1) > 1000)
+ datetick('x',15,'keeplimits','keepticks')
+ monstr = datestr(plotdat.bincenters(1),21);
+ xlabelstring = sprintf('%s start',monstr);
+else
+ xlabelstring = 'days';
+end
- % use same X ticks - with no labels
- 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)
+% only put x axis on last/bottom plot
+if (plotdat.region == plotdat.nregions)
+ xlabel(xlabelstring)
+end
- 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
+% 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');
+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','+');
+% use same X ticks - with no labels
+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)
+
+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
subplot('position',[0.48 0.01 0.04 0.04])
@@ -402,9 +410,9 @@
h = text(0.0, 0.5, string1);
set(h,'HorizontalAlignment','center', ...
- 'VerticalAlignment','middle',...
- 'Interpreter','none',...
- 'FontSize',10)
+ 'VerticalAlignment','middle',...
+ 'Interpreter','none',...
+ 'FontSize',10)
@@ -417,8 +425,8 @@
j = 0;
for i = 1:length(x.allvarnames)
- indx = findstr('time',x.allvardims{i});
- if (indx > 0)
+ indx = strfind(x.allvardims{i},'time');
+ if (indx > 0)
j = j + 1;
basenames{j} = ReturnBase(x.allvarnames{i});
@@ -426,34 +434,34 @@
end
end
-[b,i,j] = unique(basenames);
+[~,i,j] = unique(basenames);
y = cell(length(i),1);
ydims = cell(length(i),1);
for k = 1:length(i)
- disp(sprintf('%2d is %s',k,basenames{i(k)}))
- y{k} = basenames{i(k)};
-ydims{k} = basedims{i(k)};
+ fprintf('%2d is %s\n',k,basenames{i(k)})
+ y{k} = basenames{i(k)};
+ ydims{k} = basedims{i(k)};
end
function s = ReturnBase(string1)
-inds = findstr('_guess',string1);
+inds = strfind(string1,'_guess');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_analy',string1);
+inds = strfind(string1,'_analy');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_VPguess',string1);
+inds = strfind(string1,'_VPguess');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_VPanaly',string1);
+inds = strfind(string1,'_VPanaly');
if (inds > 0 )
s = string1(1:inds-1);
end
@@ -466,10 +474,10 @@
%
% 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.anl_copy(:)];
+ y.anl_copy(:)];
inds = find(isfinite(bob));
if ( isempty(inds) )
@@ -479,10 +487,10 @@
ymin = min(glommed);
ymax = max(glommed);
- if ( ymax > 1.0 )
+ if ( ymax > 1.0 )
ymin = floor(min(glommed));
ymax = ceil(max(glommed));
- elseif ( ymax < 0.0 & y.copystring == 'bias' )
+ elseif ( ymax < 0.0 && strcmp(y.copystring,'bias') )
ymax = 0.0;
end
@@ -504,8 +512,43 @@
ylimits = get(ax2,'YLim');
slope = (ylimits(2) - ylimits(1))/(Dlimits(2) - Dlimits(1));
-xtrcpt = ylimits(2) -slope*Dlimits(2);
+xtrcpt = ylimits(2) - slope*Dlimits(2);
yticks = slope*DYticks + xtrcpt;
newticklabels = num2str(round(10*yticks')/10);
+
+
+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
+% out.
+
+[variable_present, varid] = nc_var_exists(fname,varname);
+if (variable_present)
+ ncid = netcdf.open(fname);
+ value = netcdf.getVar(ncid, varid);
+else
+ value = [];
+end
+
+
+
+function value = local_nc_attget(fname,varid,varname)
+%% If the (global) attribute exists, return the value.
+% If it does not, do not throw a hissy-fit.
+
+value = [];
+if (varid == nc_global)
+ finfo = ncinfo(fname);
+ for iatt = 1:length(finfo.Attributes)
+ if (strcmp(finfo.Attributes(iatt).Name, deblank(varname)))
+ value = finfo.Attributes(iatt).Value;
+ return
+ end
+ end
+else
+ fprintf('function not supported for local variables, only global atts.\n')
+end
+
+
Modified: DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m 2013-05-09 16:38:04 UTC (rev 6123)
+++ DART/branches/development/diagnostics/matlab/plot_rmse_xxx_evolution.m 2013-05-09 21:28:46 UTC (rev 6124)
@@ -10,24 +10,23 @@
%
% 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)
% obstypestring : 'observation type' string. Optional.
-% Must match something in the netcdf
+% Must match something in the netcdf
% file 'ObservationTypes' variable.
% (ncdump -v ObservationTypes obs_diag_output.nc)
% If specified, only this observation type will be plotted.
% If not specified, all observation types incluced in the netCDF file
% will be plotted.
%
-% OUTPUT: two files will result for each observation type plotted. One is a
+% OUTPUT: two files will result for each observation type plotted. One is a
% postscript file containing a page for each level - all regions.
% The other file is a simple text file containing summary information
% about how many observations were assimilated, how many were available, etc.
% Both of these filenames contain the observation type as part of the name.
%
-%
% EXAMPLE 1 - plot the RMSE and totalspread on the same axis.
%
% fname = 'obs_diag_output.nc'; % netcdf file produced by 'obs_diag'
@@ -71,29 +70,31 @@
plotdat.copystring = copystring;
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.mlevel = local_nc_varget(fname,'mlevel');
+plotdat.plevel = local_nc_varget(fname,'plevel');
+plotdat.plevel_edges = local_nc_varget(fname,'plevel_edges');
+plotdat.hlevel = local_nc_varget(fname,'hlevel');
+plotdat.hlevel_edges = local_nc_varget(fname,'hlevel_edges');
plotdat.ncopies = length(nc_varget(fname,'copy'));
plotdat.nregions = length(nc_varget(fname,'region'));
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
-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');
+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.rat_cri = local_nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = local_nc_attget(fname, nc_global, 'input_qc_threshold');
+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
@@ -123,7 +124,7 @@
end
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.rmseindex = get_copy_index(fname,'rmse');
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
@@ -137,24 +138,24 @@
%----------------------------------------------------------------------
for ivar = 1:plotdat.nvars
-
+
% create the variable names of interest.
-
- plotdat.myvarname = plotdat.varnames{ivar};
+
+ plotdat.myvarname = plotdat.varnames{ivar};
plotdat.guessvar = sprintf('%s_guess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_analy',plotdat.varnames{ivar});
% remove any existing postscript file - will simply append each
% level as another 'page' in the .ps file.
-
- psfname = sprintf('%s_rmse_%s_evolution.ps',plotdat.varnames{ivar},plotdat.copystring);
- disp(sprintf('Removing %s from the current directory.',psfname))
+
+ 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));
- % remove any existing log file -
-
+ % remove any existing log file -
+
lgfname = sprintf('%s_rmse_%s_obscount.txt',plotdat.varnames{ivar},plotdat.copystring);
- disp(sprintf('Removing %s from the current directory.',lgfname))
+ fprintf('Removing %s from the current directory.\n',lgfname)
system(sprintf('rm %s',lgfname));
logfid = fopen(lgfname,'wt');
fprintf(logfid,'%s\n',lgfname);
@@ -164,10 +165,13 @@
guessdims = nc_var_dims(fname, plotdat.guessvar);
analydims = nc_var_dims(fname, plotdat.analyvar);
- if ( findstr('surface',guessdims{3}) > 0 )
+ if ( dimensionality == 1 ) % observations on a unit circle, no level
+ plotdat.level = 1;
+ plotdat.level_units = [];
+ elseif ( strfind(guessdims{3},'surface') > 0 )
plotdat.level = 1;
plotdat.level_units = 'surface';
- elseif ( findstr('undef',guessdims{3}) > 0 )
+ elseif ( strfind(guessdims{3},'undef') > 0 )
plotdat.level = 1;
plotdat.level_units = 'undefined';
else
@@ -176,57 +180,57 @@
end
plotdat.nlevels = length(plotdat.level);
- % Here is the tricky part. Singleton dimensions are auto-squeezed ...
+ % Here is the tricky part. Singleton dimensions are auto-squeezed ...
% single levels, single regions ...
- guess_raw = nc_varget(fname, plotdat.guessvar);
+ guess_raw = nc_varget(fname, plotdat.guessvar);
guess = reshape(guess_raw, plotdat.Nbins, plotdat.ncopies, ...
- plotdat.nlevels, plotdat.nregions);
+ plotdat.nlevels, plotdat.nregions);
- analy_raw = nc_varget(fname, plotdat.analyvar);
+ analy_raw = nc_varget(fname, plotdat.analyvar);
analy = reshape(analy_raw, plotdat.Nbins, plotdat.ncopies, ...
- plotdat.nlevels, plotdat.nregions);
+ plotdat.nlevels, plotdat.nregions);
% check to see if there is anything to plot
nposs = sum(guess(:,plotdat.Npossindex,:,:));
if ( sum(nposs(:)) < 1 )
- disp(sprintf('%s no obs for %s... skipping', plotdat.varnames{ivar}))
+ fprintf('%s no obs for %s... skipping\n', plotdat.varnames{ivar})
continue
end
-
+
for ilevel = 1:plotdat.nlevels
fprintf(logfid,'\nlevel %d %f %s\n',ilevel,plotdat.level(ilevel),plotdat.level_units);
plotdat.ges_Nqc4 = guess(:,plotdat.NQC4index ,ilevel,:);
plotdat.anl_Nqc4 = analy(:,plotdat.NQC4index ,ilevel,:);
fprintf(logfid,'DART QC == 4, prior/post %d %d\n',sum(plotdat.ges_Nqc4(:)), ...
- sum(plotdat.anl_Nqc4(:)));
+ sum(plotdat.anl_Nqc4(:)));
plotdat.ges_Nqc5 = guess(:,plotdat.NQC5index ,ilevel,:);
plotdat.anl_Nqc5 = analy(:,plotdat.NQC5index ,ilevel,:);
fprintf(logfid,'DART QC == 5, prior/post %d %d\n',sum(plotdat.ges_Nqc5(:)), ...
- sum(plotdat.anl_Nqc5(:)));
+ sum(plotdat.anl_Nqc5(:)));
plotdat.ges_Nqc6 = guess(:,plotdat.NQC6index ,ilevel,:);
plotdat.anl_Nqc6 = analy(:,plotdat.NQC6index ,ilevel,:);
fprintf(logfid,'DART QC == 6, prior/post %d %d\n',sum(plotdat.ges_Nqc6(:)), ...
- sum(plotdat.anl_Nqc6(:)));
+ sum(plotdat.anl_Nqc6(:)));
plotdat.ges_Nqc7 = guess(:,plotdat.NQC7index ,ilevel,:);
plotdat.anl_Nqc7 = analy(:,plotdat.NQC7index ,ilevel,:);
fprintf(logfid,'DART QC == 7, prior/post %d %d\n',sum(plotdat.ges_Nqc7(:)), ...
- sum(plotdat.anl_Nqc7(:)));
+ sum(plotdat.anl_Nqc7(:)));
plotdat.ges_Nposs = guess(:,plotdat.Npossindex, ilevel,:);
plotdat.anl_Nposs = analy(:,plotdat.Npossindex, ilevel,:);
fprintf(logfid,'# obs poss, prior/post %d %d\n',sum(plotdat.ges_Nposs(:)), ...
- sum(plotdat.anl_Nposs(:)));
+ sum(plotdat.anl_Nposs(:)));
plotdat.ges_Nused = guess(:,plotdat.Nusedindex, ilevel,:);
plotdat.anl_Nused = analy(:,plotdat.Nusedindex, ilevel,:);
fprintf(logfid,'# obs used, prior/post %d %d\n',sum(plotdat.ges_Nused(:)), ...
- sum(plotdat.anl_Nused(:)));
+ sum(plotdat.anl_Nused(:)));
plotdat.ges_copy = guess(:,plotdat.copyindex, ilevel,:);
plotdat.anl_copy = analy(:,plotdat.copyindex, ilevel,:);
@@ -239,18 +243,22 @@
if (plotdat.nregions > 2)
clf; orient tall
- else
+ else
clf; orient landscape
end
for iregion = 1:plotdat.nregions
- plotdat.region = iregion;
+ plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
- plotdat.title = sprintf('%s @ %d %s', ...
- plotdat.myvarname, ...
- plotdat.level(ilevel), ...
- plotdat.level_units);
+ if ( isempty(plotdat.level_units) )
+ plotdat.title = plotdat.myvarname;
+ else
+ plotdat.title = sprintf('%s @ %d %s', ...
+ plotdat.myvarname, ...
+ plotdat.level(ilevel), ...
+ plotdat.level_units);
+ end
myplot(plotdat);
end
@@ -271,131 +279,131 @@
function myplot(plotdat)
- % Interlace the [ges,anl] to make a sawtooth plot.
- % By this point, the middle two dimensions are singletons.
- cg = plotdat.ges_copy(:,:,:,plotdat.region);
- ca = plotdat.anl_copy(:,:,:,plotdat.region);
- other = reshape([cg ca]',2*plotdat.Nbins,1);
+% Interlace the [ges,anl] to make a sawtooth plot.
+% By this point, the middle two dimensions are singletons.
+cg = plotdat.ges_copy(:,:,:,plotdat.region);
+ca = plotdat.anl_copy(:,:,:,plotdat.region);
+other = reshape([cg ca]',2*plotdat.Nbins,1);
- mg = plotdat.ges_rmse(:,:,:,plotdat.region);
- ma = plotdat.anl_rmse(:,:,:,plotdat.region);
- rmse = reshape([mg ma]',2*plotdat.Nbins,1);
+mg = plotdat.ges_rmse(:,:,:,plotdat.region);
+ma = plotdat.anl_rmse(:,:,:,plotdat.region);
+rmse = reshape([mg ma]',2*plotdat.Nbins,1);
- g = plotdat.ges_Nposs(:,:,:,plotdat.region);
- a = plotdat.anl_Nposs(:,:,:,plotdat.region);
- nobs_poss = reshape([g a]',2*plotdat.Nbins,1);
+g = plotdat.ges_Nposs(:,:,:,plotdat.region);
+a = plotdat.anl_Nposs(:,:,:,plotdat.region);
+nobs_poss = reshape([g a]',2*plotdat.Nbins,1);
- g = plotdat.ges_Nused(:,:,:,plotdat.region);
- a = plotdat.anl_Nused(:,:,:,plotdat.region);
- nobs_used = reshape([g a]',2*plotdat.Nbins,1);
+g = plotdat.ges_Nused(:,:,:,plotdat.region);
+a = plotdat.anl_Nused(:,:,:,plotdat.region);
+nobs_used = reshape([g a]',2*plotdat.Nbins,1);
- tg = plotdat.bincenters;
- ta = plotdat.bincenters;
- t = reshape([tg ta]',2*plotdat.Nbins,1);
+tg = plotdat.bincenters;
+ta = plotdat.bincenters;
+t = reshape([tg ta]',2*plotdat.Nbins,1);
- % Determine some quantities for the legend
- nobs = sum(nobs_used);
- if ( nobs > 1 )
- mean_pr_rmse = mean(mg(isfinite(mg)));
- mean_po_rmse = mean(ma(isfinite(ma)));
- mean_pr_other = mean(cg(isfinite(cg)));
- mean_po_other = mean(ca(isfinite(ca)));
- else
- mean_pr_rmse = NaN;
- mean_po_rmse = NaN;
- mean_pr_other = NaN;
- mean_po_other = NaN;
- end
+% Determine some quantities for the legend
+nobs = sum(nobs_used);
+if ( nobs > 1 )
+ mean_pr_rmse = mean(mg(isfinite(mg)));
+ mean_po_rmse = mean(ma(isfinite(ma)));
+ mean_pr_other = mean(cg(isfinite(cg)));
+ mean_po_other = mean(ca(isfinite(ca)));
+else
+ mean_pr_rmse = NaN;
+ mean_po_rmse = NaN;
+ mean_pr_other = NaN;
+ mean_po_other = NaN;
+end
- string_rmse = sprintf('%s pr=%.5g, po=%.5g','rmse', mean_pr_rmse, mean_po_rmse);
- string_other = sprintf('%s pr=%.5g, po=%.5g', plotdat.copystring, ...
- mean_pr_other, mean_po_other);
- plotdat.subtitle = sprintf('%s %s %s',plotdat.myregion,string_rmse, string_other);
+string_rmse = sprintf('%s pr=%.5g, po=%.5g','rmse', mean_pr_rmse, mean_po_rmse);
+string_other = sprintf('%s pr=%.5g, po=%.5g', plotdat.copystring, ...
+ mean_pr_other, mean_po_other);
+plotdat.subtitle = sprintf('%s %s %s',plotdat.myregion,string_rmse, string_other);
- % Plot the rmse and 'xxx' on the same (left) axis.
- % The observation count will use the axis on the right.
- % We want to suppress the 'auto' feature of the axis labelling,
- % so we manually set some values that normally
- % don't need to be set.
-
- ax1 = subplot(plotdat.nregions,1,plotdat.region);
+% Plot the rmse and 'xxx' on the same (left) axis.
+% The observation count will use the axis on the right.
+% 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.
- h1 = plot(t,rmse,'k+-',t,other,'ro-','LineWidth',plotdat.linewidth);
- h = legend(h1,'rmse', plotdat.copystring);
- legend(h,'boxoff','Interpreter','none')
+ax1 = subplot(plotdat.nregions,1,plotdat.region);
- axlims = axis;
- axlims = [axlims(1:2) plotdat.Yrange];
- axis(axlims)
+h1 = plot(t,rmse,'k+-',t,other,'ro-','LineWidth',plotdat.linewidth);
+h = legend(h1,'rmse', plotdat.copystring);
+legend(h,'boxoff','Interpreter','none')
- plotdat.ylabel{1} = plotdat.myregion;
- switch lower(plotdat.copystring)
- case 'bias'
- % plot a zero-bias line
- h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
- set(h4,'LineWidth',1.5,'LineSTyle',':')
- plotdat.ylabel{2} = sprintf('rmse and %s (%s)',plotdat.copystring,plotdat.biasconv);
- otherwise
- plotdat.ylabel{2} = sprintf('rmse and %s',plotdat.copystring);
- end
-
- % hokey effort to decide to plot months/days vs. daynum vs.
- ttot = plotdat.bincenters(plotdat.Nbins) - plotdat.bincenters(1) + 1;
-
- if ((plotdat.bincenters(1) > 1000) && (ttot > 5))
- datetick('x',6,'keeplimits','keepticks');
- monstr = datestr(plotdat.bincenters(1),21);
- xlabelstring = sprintf('month/day - %s start',monstr);
- elseif (plotdat.bincenters(1) > 1000)
- datetick('x',15,'keeplimits','keepticks')
- monstr = datestr(plotdat.bincenters(1),21);
- xlabelstring = sprintf('%s start',monstr);
- else
- xlabelstring = 'days';
- end
+axlims = axis;
+axlims = [axlims(1:2) plotdat.Yrange];
+axis(axlims)
- % only put x axis on last/bottom plot
- if (plotdat.region == plotdat.nregions)
- xlabel(xlabelstring)
- end
+plotdat.ylabel{1} = plotdat.myregion;
+switch lower(plotdat.copystring)
+ case 'bias'
+ % plot a zero-bias line
+ h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
+ set(h4,'LineWidth',1.5,'LineSTyle',':')
+ plotdat.ylabel{2} = sprintf('rmse and %s (%s)',plotdat.copystring,plotdat.biasconv);
+ otherwise
+ plotdat.ylabel{2} = sprintf('rmse and %s',plotdat.copystring);
+end
- % more annotation ...
+% hokey effort to decide to plot months/days vs. daynum vs.
+ttot = plotdat.bincenters(plotdat.Nbins) - plotdat.bincenters(1) + 1;
- 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');
- 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','+');
+if ((plotdat.bincenters(1) > 1000) && (ttot > 5))
+ datetick('x',6,'keeplimits','keepticks');
+ monstr = datestr(plotdat.bincenters(1),21);
+ xlabelstring = sprintf('month/day - %s start',monstr);
+elseif (plotdat.bincenters(1) > 1000)
+ datetick('x',15,'keeplimits','keepticks')
+ monstr = datestr(plotdat.bincenters(1),21);
+ xlabelstring = sprintf('%s start',monstr);
+else
+ xlabelstring = 'days';
+end
- % use same X ticks - with no labels
- 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)
+% only put x axis on last/bottom plot
+if (plotdat.region == plotdat.nregions)
+ xlabel(xlabelstring)
+end
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list