[Dart-dev] [6327] DART/trunk/models/POP/matlab: Incorporating changes to be consistent with the DART/diagnostics/matlab
nancy at ucar.edu
nancy at ucar.edu
Mon Jul 29 15:10:13 MDT 2013
Revision: 6327
Author: thoar
Date: 2013-07-29 15:10:13 -0600 (Mon, 29 Jul 2013)
Log Message:
-----------
Incorporating changes to be consistent with the DART/diagnostics/matlab
versions of these files. These now use the new tick-axis labelling method.
These differ from the model-agnostic versions in the treatement of
the number of possible observations. If the namelist is set to ignore
observations, or the observations are too close to land for the model
to apply the forward operator, then those observations are not really
'possible'. FYI - the only time get_expected_obs() can fail (i.e. DART QC == 4)
with POP is when it is trying to extrapolate to a 'dry' location.
If this changes, then my definition of 'possible' needs to change.
Modified Paths:
--------------
DART/trunk/models/POP/matlab/plot_evolution.m
DART/trunk/models/POP/matlab/plot_rmse_xxx_evolution.m
DART/trunk/models/POP/matlab/plot_spread_xxx_evolution.m
DART/trunk/models/POP/matlab/two_experiments_evolution.m
-------------- next part --------------
Modified: DART/trunk/models/POP/matlab/plot_evolution.m
===================================================================
--- DART/trunk/models/POP/matlab/plot_evolution.m 2013-07-29 21:06:18 UTC (rev 6326)
+++ DART/trunk/models/POP/matlab/plot_evolution.m 2013-07-29 21:10:13 UTC (rev 6327)
@@ -10,23 +10,24 @@
%
% 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 evolution of the bias for all observation types, all levels
%
% fname = 'obs_diag_output.nc'; % netcdf file produced by 'obs_diag'
@@ -57,7 +58,6 @@
error('wrong number of arguments ... ')
end
-
if (exist(fname,'file') ~= 2)
error('file/fname <%s> does not exist',fname)
end
@@ -68,27 +68,30 @@
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');
-if (plotdat.nregions == 1)
+% 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.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.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
@@ -96,8 +99,14 @@
timeunits = nc_attget(fname,'time','units');
timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
timeorigin = datenum(timebase(1),timebase(2),timebase(3));
-skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
-iskip = time_to_skip(3) + skip_seconds/86400;
+if ( isempty(time_to_skip) == 1)
+ iskip = 0;
+elseif ( numel(time_to_skip) == 6)
+ skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
+ iskip = time_to_skip(3) + skip_seconds/86400;
+else
+ error('time_to_skip variable has unusual length. Should be either 0 or 6.')
+end
plotdat.bincenters = plotdat.bincenters + timeorigin;
plotdat.binedges = plotdat.binedges + timeorigin;
@@ -117,8 +126,7 @@
plotdat.nvars = nvars;
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');
@@ -131,24 +139,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);
@@ -158,10 +166,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
@@ -170,16 +181,16 @@
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
nextrap = sum(guess(:,plotdat.NQC4index ,:,:)) + ...
@@ -187,32 +198,32 @@
nposs = sum(guess(:,plotdat.Npossindex, :,:)) - nextrap;
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,:);
@@ -222,12 +233,12 @@
plotdat.ges_Nposs = plotdat.ges_Nposs - plotdat.ges_Nqc4 - plotdat.ges_Nqc6;
plotdat.anl_Nposs = plotdat.anl_Nposs - plotdat.anl_Nqc4 - plotdat.anl_Nqc6;
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,19 +248,23 @@
% plot by region
if (plotdat.nregions > 2)
- clf; orient tall
- else
- clf; orient landscape
+ clf; orient tall; wysiwyg
+ else
+ clf; orient landscape; wysiwyg
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
@@ -264,129 +279,132 @@
end
end
-%----------------------------------------------------------------------
+%=====================================================================
% 'Helper' functions
-%----------------------------------------------------------------------
+%=====================================================================
+
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);
+%% The prior and posterior are plotted as separate items.
+% By this point, the middle two dimensions are singletons.
+cg = plotdat.ges_copy(:,:,:,plotdat.region);
+ca = plotdat.anl_copy(:,:,:,plotdat.region);
- 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')
+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
+% get the range of the existing axis and replace with
+% replace y axis values
+% reset the axes limits
- % only put x axis on last/bottom plot
- if (plotdat.region == plotdat.nregions)
- xlabel(xlabelstring)
- end
+axlims = axis;
+axlims = [axlims(1:2) plotdat.Yrange];
+axis(axlims)
- % more annotation ...
+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.0,'LineSTyle','-.')
+ plotdat.ylabel{2} = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
+ otherwise
+ plotdat.ylabel{2} = sprintf('%s',plotdat.copystring);
+end
- 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','+');
+% hokey effort to decide to plot months/days vs. daynum vs.
+ttot = plotdat.bincenters(plotdat.Nbins) - plotdat.bincenters(1) + 1;
- % 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)
+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
- set(get(ax2,'Ylabel'),'String','# of obs : o=poss, +=used')
- set(get(ax1,'Ylabel'),'String',plotdat.ylabel)
- set(ax1,'Position',get(ax2,'Position'))
- grid
+% only put x axis on last/bottom plot
+if (plotdat.region == plotdat.nregions)
+ xlabel(xlabelstring)
+end
+% more annotation ...
+if (plotdat.region == 1)
+ title({plotdat.title, plotdat.subtitle}, ...
+ 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
+ BottomAnnotation(plotdat.fname)
+else
+ title(plotdat.subtitle, ...
+ 'Interpreter', 'none', 'Fontsize', 12, 'FontWeight', 'bold')
+end
+% 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
+% Make sure that values at ticks are whole numbers.
+matchingYticks(ax1,ax2);
+
+set(get(ax2,'Ylabel'),'String','# of obs : o=poss, +=used')
+set(get(ax1,'Ylabel'),'String',plotdat.ylabel,'Interpreter','none')
+
+
+%=====================================================================
+
+
function BottomAnnotation(main)
-% annotates the directory containing the data being plotted
+%% annotates the full path of the file being plotted
subplot('position',[0.48 0.01 0.04 0.04])
axis off
fullname = which(main); % Could be in MatlabPath
@@ -403,14 +421,16 @@
h = text(0.0, 0.5, string1);
set(h,'HorizontalAlignment','center', ...
- 'VerticalAlignment','middle',...
- 'Interpreter','none',...
- 'FontSize',10)
+ 'VerticalAlignment','middle',...
+ 'Interpreter','none',...
+ 'FontSize',8)
+%=====================================================================
+
function [y,ydims] = FindTemporalVars(x)
-% Returns UNIQUE (i.e. base) temporal variable names
+%% Returns UNIQUE (i.e. base) temporal variable names
if ( ~(isfield(x,'allvarnames') && isfield(x,'allvardims')))
error('Doh! no ''allvarnames'' and ''allvardims'' components')
end
@@ -418,56 +438,62 @@
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});
- basedims{j} = x.allvardims{i};
+ basedims{ j} = x.allvardims{i};
end
end
-[b,i,j] = unique(basenames);
+[~,i,j] = unique(basenames);
y = cell(length(i),1);
ydims = cell(length(i),1);
-for j = 1:length(i)
- disp(sprintf('%2d is %s',j,basenames{j}))
- y{j} = basenames{j};
-ydims{j} = basedims{j};
+for k = 1:length(i)
+ 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);
+%% Pick off the variable name.
+inds = strfind(string1,'_guess');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-
-inds = findstr('_analy',string1);
+inds = strfind(string1,'_analy');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-
-inds = findstr('_VPguess',string1);
+inds = strfind(string1,'_VPguess');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
-
-inds = findstr('_VPanaly',string1);
+inds = strfind(string1,'_VPanaly');
if (inds > 0 )
s = string1(1:inds-1);
+ return
end
+%=====================================================================
+
function x = FindRange(y)
-% Trying to pick 'nice' limits for plotting.
+%% Trying to pick 'nice' limits for plotting.
% Completely ad hoc ... and not well posed.
%
% In this scope, y is bounded from below by 0.0
%
-% If the numbers are very small ...
+% If the numbers are very small ...
bob = [y.ges_copy(:) ; ...
y.anl_copy(:)];
@@ -480,10 +506,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
@@ -493,23 +519,46 @@
end
+%=====================================================================
-function [yticks newticklabels] = matchingYticks(ax1, ax2)
-%% This takes the existing Y ticks from ax1 (presumed nice)
-% and determines the matching labels for ax2 so we can keep
-% at least one of the axes looking nice.
-Dlimits = get(ax1,'YLim');
-DYticks = get(ax1,'YTick');
-nYticks = length(DYticks);
-ylimits = get(ax2,'YLim');
+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.
-slope = (ylimits(2) - ylimits(1))/(Dlimits(2) - Dlimits(1));
-xtrcpt = ylimits(2) -slope*Dlimits(2);
+[variable_present, varid] = nc_var_exists(fname,varname);
+if (variable_present)
+ ncid = netcdf.open(fname);
+ value = netcdf.getVar(ncid, varid);
+ netcdf.close(ncid)
+else
+ value = [];
+end
-yticks = slope*DYticks + xtrcpt;
-newticklabels = num2str(round(10*yticks')/10);
+%=====================================================================
+
+
+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
+
+
+
% <next few lines under version control, do not edit>
% $URL$
% $Revision$
Modified: DART/trunk/models/POP/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/trunk/models/POP/matlab/plot_rmse_xxx_evolution.m 2013-07-29 21:06:18 UTC (rev 6326)
+++ DART/trunk/models/POP/matlab/plot_rmse_xxx_evolution.m 2013-07-29 21:10:13 UTC (rev 6327)
@@ -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.
@@ -67,27 +67,30 @@
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');
-if (plotdat.nregions == 1)
+% 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.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.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
@@ -95,8 +98,14 @@
timeunits = nc_attget(fname,'time','units');
timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
timeorigin = datenum(timebase(1),timebase(2),timebase(3));
-skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
-iskip = time_to_skip(3) + skip_seconds/86400;
+if ( isempty(time_to_skip) == 1)
+ iskip = 0;
+elseif ( numel(time_to_skip) == 6)
+ skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
+ iskip = time_to_skip(3) + skip_seconds/86400;
+else
+ error('time_to_skip variable has unusual length. Should be either 0 or 6.')
+end
plotdat.bincenters = plotdat.bincenters + timeorigin;
plotdat.binedges = plotdat.binedges + timeorigin;
@@ -116,8 +125,7 @@
plotdat.nvars = nvars;
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');
@@ -131,24 +139,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))
+ 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);
@@ -158,10 +166,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
@@ -170,16 +181,16 @@
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
nextrap = sum(guess(:,plotdat.NQC4index ,:,:)) + ...
@@ -187,32 +198,32 @@
nposs = sum(guess(:,plotdat.Npossindex, :,:)) - nextrap;
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,:);
@@ -222,12 +233,12 @@
plotdat.ges_Nposs = plotdat.ges_Nposs - plotdat.ges_Nqc4 - plotdat.ges_Nqc6;
plotdat.anl_Nposs = plotdat.anl_Nposs - plotdat.anl_Nqc4 - plotdat.anl_Nqc6;
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,19 +250,23 @@
% plot by region
if (plotdat.nregions > 2)
- clf; orient tall
- else
- clf; orient landscape
+ clf; orient tall; wysiwyg
+ else
+ clf; orient landscape; wysiwyg
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
@@ -266,139 +281,142 @@
end
end
-%----------------------------------------------------------------------
+%=====================================================================
% 'Helper' functions
-%----------------------------------------------------------------------
+%=====================================================================
+
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')
+ax1 = subplot(plotdat.nregions,1,plotdat.region);
- axlims = axis;
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list