[Dart-dev] [7605] DART/trunk/diagnostics: The matlab plotting scripts annotate the fact that trusted observations
nancy at ucar.edu
nancy at ucar.edu
Fri Feb 20 14:57:32 MST 2015
Revision: 7605
Author: thoar
Date: 2015-02-20 14:57:32 -0700 (Fri, 20 Feb 2015)
Log Message:
-----------
The matlab plotting scripts annotate the fact that trusted observations
were used. The number of observations assimilated remains unchanged, but
the values plotted reflect the fact that outlier observations were used
by obs_diag.
obs_diag now puts an attribute on each of the trusted variables to
identify which ones are trusted or not.
obs_diag also counts up the number of pathological observations that
have a DART QC of 7 but a missing posterior mean (would be a DART QC 4).
The count for the N_DARTqc_4 and N_DARTqc_7 are more accurate now than
previously. The metrics (rmse, bias, spread, etc) are unchanged.
This feature only affects the calculation of the metrics for TRUSTED observations.
If the observation is TRUSTED, then even the outlier obs are used to
calculate the metrics. For those few with missing posterior means, every now
and again there are CRAZY posterior metrics because a MISSING posterior mean
is leaking into the calculation. To compensate for this effect, there are
two lines of code that can be enabled to force the QC to be a 4. Those obs
cannot be used in the calculation of TRUSTED observations.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
DART/trunk/diagnostics/matlab/plot_evolution.m
DART/trunk/diagnostics/matlab/plot_profile.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
DART/trunk/diagnostics/matlab/two_experiments_evolution.m
DART/trunk/diagnostics/matlab/two_experiments_profile.m
DART/trunk/diagnostics/threed_sphere/obs_diag.f90
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2015-02-20 21:57:03 UTC (rev 7604)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2015-02-20 21:57:32 UTC (rev 7605)
@@ -1,22 +1,52 @@
-function plotdat = plot_bias_xxx_profile(fname,copystring)
+function plotdat = plot_bias_xxx_profile(fname, copy, varargin)
%% plot_bias_xxx_profile plots the vertical profile of the observation-space quantities for all possible levels, all possible variables.
% Part of the observation-space diagnostics routines.
%
% 'obs_diag' produces a netcdf file containing the diagnostics.
+% obs_diag condenses the obs_seq.final information into summaries for a few specified
+% regions - on a level-by-level basis.
%
-% USAGE: plotdat = plot_bias_xxx_profile(fname,copystring);
+% The number of observations possible reflects only those observations
+% that have incoming QC values of interest. Any observation with a DART
+% QC of 5 or 6 is not considered 'possible' for the purpose of this graphic.
%
-% fname : netcdf file produced by 'obs_diag'
-% copystring : 'copy' string == quantity of interest. These
-% can be any of the ones available in the netcdf
-% file 'CopyMetaData' variable.
+% NOTE: if the observation was designated as a TRUSTED observation in the
+% obs_diag program, the observations that were rejected by the outlier
+% threshhold STILL PARTICIPATE in the calculation of the rmse, spread, etc.
+% The _values_ plotted by plot_profile reflect that. The number of observations
+% "used" becomes unclear. The number of observations used (designated by the
+% asterisk) is ALWAYS the number of observations successfully assimilated.
+% For TRUSTED observations, this is different than the number used to calculate
+% bias, rmse, spread, etc.
+%
+% USAGE: plotdat = plot_bias_xxx_profile(fname, copy);
+%
+% fname : netcdf file produced by 'obs_diag'
+%
+% copy : string defining the metric of interest. 'rmse', 'spread', etc.
+% Possible values are available in the netcdf 'CopyMetaData' variable.
% (ncdump -v CopyMetaData obs_diag_output.nc)
%
-% EXAMPLE:
+% obsname : Optional. If present, The strings of each observation type to plot.
+% Each observation type will be plotted in a separate graphic.
+% Default is to plot all available observation types.
%
-% fname = 'obs_diag_output.nc'; % netcdf file produced by 'obs_diag'
-% copystring = 'totalspread'; % 'copy' string == quantity of interest
-% plotdat = plot_bias_xxx_profile(fname,copystring);
+% OUTPUT: 'plotdat' is a structure containing what was plotted.
+% A .pdf of each graphic is created. Each .pdf has a name that
+% reflects the variable, quantity, and region being plotted.
+%
+% EXAMPLE 1: All the observation types possible are plotted in separate figures.
+%
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% plotdat = plot_bias_xxx_profile(fname, copy);
+%
+% EXAMPLE 2: Just a single observation type.
+%
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% obsname = 'RADIOSONDE_U_WIND_COMPONENT';
+% plotdat = plot_bias_xxx_profile(fname, copy, 'obsname', obsname);
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -28,6 +58,24 @@
% Decode,Parse,Check the input
%---------------------------------------------------------------------
+default_obsname = 'none';
+p = inputParser;
+
+addRequired(p,'fname', at ischar);
+addRequired(p,'copy', at ischar);
+addParamValue(p,'obsname',default_obsname, at ischar);
+parse(p, fname, copy, varargin{:});
+
+% if you want to echo the input
+% disp(['fname : ', p.Results.fname])
+% disp(['copy : ', p.Results.copy])
+% disp(['obsname : ', p.Results.obsname])
+
+if ~isempty(fieldnames(p.Unmatched))
+ disp('Extra inputs:')
+ disp(p.Unmatched)
+end
+
if (exist(fname,'file') ~= 2)
error('file/fname <%s> does not exist',fname)
end
@@ -37,7 +85,7 @@
%---------------------------------------------------------------------
plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.copystring = copy;
plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
@@ -85,13 +133,13 @@
plotdat.toff = plotdat.binedges(1) + iskip;
plotdat.timespan = sprintf('%s through %s', datestr(plotdat.toff), ...
datestr(max(plotdat.binedges(:))));
-plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copystring);
+plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copy);
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindVerticalVars(plotdat);
plotdat.nvars = length(plotdat.varnames);
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copy);
plotdat.biasindex = get_copy_index(fname,'bias');
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
@@ -104,14 +152,28 @@
% Loop around (copy-level-region) observation types
%----------------------------------------------------------------------
-for ivar = 1:plotdat.nvars
+% Either use all the variables or just the one optionally specified.
+if strcmp(p.Results.obsname,'none')
+ varlist = 1:plotdat.nvars;
+else
+ varlist = find (strcmpi(p.Results.obsname,plotdat.varnames));
+ if isempty(varlist)
+ error('%s is not in the list of observations',p.Results.obsname)
+ end
+end
+
+for ivar = varlist
+
% create the variable names of interest.
plotdat.myvarname = plotdat.varnames{ivar};
plotdat.guessvar = sprintf('%s_VPguess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_VPanaly',plotdat.varnames{ivar});
+ plotdat.trusted = nc_read_att(fname, plotdat.guessvar, 'TRUSTED');
+ if (isempty(plotdat.trusted)), plotdat.trusted = 'NO'; end
+
% get appropriate vertical coordinate variable
guessdims = nc_var_dims( fname, plotdat.guessvar);
@@ -132,7 +194,7 @@
continue
end
- [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
+ [level_org, level_units, nlevels, level_edges, Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
plotdat.level_org = level_org;
plotdat.level_units = level_units;
plotdat.nlevels = nlevels;
@@ -225,7 +287,7 @@
% plot by region - each in its own figure.
for iregion = 1:plotdat.nregions
- figure(iregion); clf; orient(figuredata.orientation); wysiwyg
+ figure(iregion); clf(iregion); orient(figuredata.orientation); wysiwyg
plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
myplot(plotdat, figuredata);
@@ -313,6 +375,21 @@
h = legend(h1, str_bias_pr, str_bias_po, str_other_pr, str_other_po, 'Location', 'NorthWest');
set(h,'Interpreter','none','Box','off')
+% If the observation is trusted, reference that somehow
+switch lower(plotdat.trusted)
+ case 'true'
+ axlims = axis;
+ tx = axlims(2) + (axlims(2) - axlims(1))/20;
+ if strcmpi('normal',plotdat.YDir)
+ ty = plotdat.Yrange(1);
+ else
+ ty = plotdat.Yrange(2);
+ end
+ h = text(tx,ty,'TRUSTED. Values include outlying observations.');
+ set(h,'FontSize',20,'Rotation',90,'VerticalAlignment','middle')
+ otherwise
+end
+
% Create another axes to use for plotting the observation counts
ax2 = axes('position',get(ax1,'Position'), ...
@@ -341,7 +418,7 @@
set(get(ax1,'Xlabel'),'String',{plotdat.xlabel, plotdat.timespan}, ...
'Interpreter','none','FontSize',figdata.fontsize)
set(get(ax2,'Xlabel'),'String', ...
- ['# of obs (o=poss, \ast=used) x' int2str(uint32(xscale))],'FontSize',figdata.fontsize)
+ ['# of obs (o=possible, \ast=assimilated) x' int2str(uint32(xscale))],'FontSize',figdata.fontsize)
title({plotdat.myregion, plotdat.myvarname}, ...
'Interpreter', 'none', 'FontSize', figdata.fontsize, 'FontWeight', 'bold')
@@ -409,7 +486,7 @@
%=====================================================================
-function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname)
+function [level_org, level_units, nlevels, level_edges, Yrange] = FindVerticalInfo(fname,varname)
%% Find the vertical dimension and harvest some info
varinfo = nc_getvarinfo(fname,varname);
@@ -577,7 +654,7 @@
position = [0.15 0.12 0.7 0.75];
linewidth = 2.0;
-figdata = struct('expcolors', {{'k','r','g','m','b','c','y'}}, ...
+figdata = struct('expcolors', {{'k','r','b','m','g','c','y'}}, ...
'expsymbols', {{'o','s','d','p','h','s','*'}}, ...
'prpolines', {{'-','--'}}, 'position', position, ...
'fontsize',fontsize, 'orientation',orientation, ...
Modified: DART/trunk/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_evolution.m 2015-02-20 21:57:03 UTC (rev 7604)
+++ DART/trunk/diagnostics/matlab/plot_evolution.m 2015-02-20 21:57:32 UTC (rev 7605)
@@ -1,4 +1,4 @@
-function plotdat = plot_evolution(fname, copystring, varargin)
+function plotdat = plot_evolution(fname, copy, varargin)
%% plot_evolution plots the temporal evolution of the observation-space quantities for all possible levels, all possible variables.
% Part of the observation-space diagnostics routines.
%
@@ -6,42 +6,64 @@
% obs_diag condenses the obs_seq.final information into summaries for a few specified
% regions - on a level-by-level basis.
%
-% USAGE: plotdat = plot_evolution(fname, copystring [,varargin]);
+% The number of observations possible reflects only those observations
+% that have incoming QC values of interest. Any observation with a DART
+% QC of 5 or 6 is not considered 'possible' for the purpose of this graphic.
%
-% fname : netcdf file produced by 'obs_diag'
-% copystring : 'copy' string == quantity of interest. These
-% 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
-% 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.
+% NOTE: if the observation was designated as a TRUSTED observation in the
+% obs_diag program, the observations that were rejected by the outlier
+% threshhold STILL PARTICIPATE in the calculation of the rmse, spread, etc.
+% The _values_ plotted by plot_profile reflect that. The number of observations
+% "used" becomes unclear. The number of observations used (designated by the
+% asterisk) is ALWAYS the number of observations successfully assimilated.
+% For TRUSTED observations, this is different than the number used to calculate
+% bias, rmse, spread, etc.
%
-% OUTPUT: two files will result for each observation type plotted. One is a
-% postscript file containing a page for each level - all regions.
+% USAGE: plotdat = plot_evolution(fname, copy);
+%
+% fname : netcdf file produced by 'obs_diag'
+%
+% copy : string defining the metric of interest. 'rmse', 'spread', etc.
+% Possible values are available in the netcdf 'CopyMetaData' variable.
+% (ncdump -v CopyMetaData obs_diag_output.nc)%
+%
+% obsname : Optional. If present, The strings of each observation type to plot.
+% Each observation type will be plotted in a separate graphic.
+% Default is to plot all available observation types.
+%
+% level : Optional. 'level' index. Default is to plot all levels.
+%
+% range : Optional. 'range' of the value being plotted. Default is to
+% automatically determine range based on the data values.
+%
+% OUTPUT: 'plotdat' is a structure containing what was last plotted.
+% A postscript file containing a page for each level - each region.
% 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.
+% Both of these filenames contain the observation type,
+% copy and region 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'
-% copystring = 'bias'; % 'copy' string == quantity of interest
-% plotdat = plot_evolution(fname, copystring);
+% fname = 'obs_diag_output.nc';
+% copy = 'bias';
+% plotdat = plot_evolution(fname, copy);
%
%
-% EXAMPLE 2 - plot the evolution of the rmse for just the radiosonde temperature observations
-% This requires that the 'RADIOSONDE_TEMPERATURE' is one of the known observation
-% types in the netCDF file.
+% EXAMPLE 2 - plot the evolution of the rmse for just the radiosonde temperature obs
+% This requires that the 'RADIOSONDE_TEMPERATURE' is one of the known
+% observation types in the netCDF file.
%
-% fname = 'obs_diag_output.nc'; % netcdf file produced by 'obs_diag'
-% copystring = 'rmse'; % 'copy' string == quantity of interest
-% plotdat = plot_evolution(fname, copystring, 'RADIOSONDE_TEMPERATURE');
+% fname = 'obs_diag_output.nc';
+% copy = 'rmse';
+% plotdat = plot_evolution(fname, copy, 'obsname', 'RADIOSONDE_TEMPERATURE');
+%
+%
+% EXAMPLE 3 - plot the evolution of the rmse for just the radiosonde temperature obs
+% for the 4th level and force the vertical axis of the 'rmse' to be 0,10
+%
+% plotdat = plot_evolution(fname, copy, 'obsname', 'RADIOSONDE_TEMPERATURE', ...
+% 'level', 4, 'range', [0 10]);
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -49,13 +71,39 @@
%
% DART $Id$
-if nargin == 2
+default_level = -1;
+default_obsname = 'none';
+default_range = [NaN NaN];
+p = inputParser;
+
+addRequired(p,'fname', at ischar);
+addRequired(p,'copy', at ischar);
+addParamValue(p,'obsname',default_obsname, at ischar);
+addParamValue(p,'range',default_range, at isnumeric);
+addParamValue(p,'level',default_level, at isnumeric);
+parse(p, fname, copy, varargin{:});
+
+% if you want to echo the input
+% fprintf('fname : %s\n', p.Results.fname)
+% fprintf('copy : %s\n', p.Results.copy)
+% fprintf('obsname : %s\n', p.Results.obsname)
+% fprintf('level : %d\n', p.Results.level)
+% fprintf('range : %f %f \n', p.Results.range)
+
+if ~isempty(fieldnames(p.Unmatched))
+ disp('Extra inputs:')
+ disp(p.Unmatched)
+end
+
+if (numel(p.Results.range) ~= 2)
+ error('range must be an array of length two ... [bottom top]')
+end
+
+if strcmp(p.Results.obsname,'none')
nvars = 0;
-elseif nargin == 3
- varname = varargin{1};
+else
+ obsname = p.Results.obsname;
nvars = 1;
-else
- error('wrong number of arguments ... ')
end
if (exist(fname,'file') ~= 2)
@@ -67,7 +115,7 @@
%---------------------------------------------------------------------
plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.copystring = copy;
plotdat.bincenters = nc_varget(fname,'time');
plotdat.binedges = nc_varget(fname,'time_bounds');
plotdat.mlevel = local_nc_varget(fname,'mlevel');
@@ -85,15 +133,15 @@
plotdat.region_names = deblank(plotdat.region_names');
end
-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');
+dimensionality = nc_read_att(fname, nc_global, 'LocationRank');
+plotdat.binseparation = nc_read_att(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_read_att(fname, nc_global, 'bin_width');
+time_to_skip = nc_read_att(fname, nc_global, 'time_to_skip');
+plotdat.lonlim1 = nc_read_att(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_read_att(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_read_att(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_read_att(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_read_att(fname, nc_global, 'bias_convention');
% Coordinate between time types and dates
@@ -122,11 +170,11 @@
[plotdat.varnames, plotdat.vardims] = FindTemporalVars(plotdat);
plotdat.nvars = length(plotdat.varnames);
else
- plotdat.varnames{1} = varname;
+ plotdat.varnames{1} = obsname;
plotdat.nvars = nvars;
end
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copy);
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
plotdat.NQC4index = get_copy_index(fname,'N_DARTqc_4');
@@ -148,6 +196,9 @@
plotdat.guessvar = sprintf('%s_guess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_analy',plotdat.varnames{ivar});
+ plotdat.trusted = nc_read_att(fname, plotdat.guessvar, 'TRUSTED');
+ if (isempty(plotdat.trusted)), plotdat.trusted = 'NO'; end
+
% remove any existing postscript file - will simply append each
% level as another 'page' in the .ps file.
@@ -218,8 +269,14 @@
continue
end
- for ilevel = 1:plotdat.nlevels
+ if (p.Results.level < 0)
+ wantedlevels = [1:plotdat.nlevels];
+ else
+ wantedlevels = p.Results.level;
+ end
+ for ilevel = wantedlevels
+
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,:);
@@ -256,12 +313,16 @@
plotdat.ges_copy = guess(:,plotdat.copyindex, ilevel,:);
plotdat.anl_copy = analy(:,plotdat.copyindex, ilevel,:);
- plotdat.Yrange = FindRange(plotdat);
+ if isnan(p.Results.range(1))
+ plotdat.Yrange = FindRange(plotdat);
+ else
+ plotdat.Yrange = p.Results.range;
+ end
% plot each region, each level to a separate figure
for iregion = 1:plotdat.nregions
- figure(iregion); clf; orient(figuredata.orientation); wysiwyg
+ figure(iregion); clf(iregion); orient(figuredata.orientation); wysiwyg
plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
@@ -346,6 +407,7 @@
axlims = [axlims(1:2) plotdat.Yrange];
axis(axlims)
+
switch lower(plotdat.copystring)
case 'bias'
% plot a zero-bias line
@@ -375,7 +437,7 @@
title({plotdat.myregion, plotdat.title, plotdat.subtitle}, ...
'Interpreter', 'none', 'Fontsize', figdata.fontsize, 'FontWeight', 'bold')
-BottomAnnotation(plotdat.fname)
+BottomAnnotation(plotdat)
% create a separate scale for the number of observations
ax2 = axes('position',get(ax1,'Position'), ...
@@ -398,7 +460,7 @@
set(get(ax1,'Ylabel'), 'String', plotdat.ylabel, ...
'Interpreter','none','FontSize',figdata.fontsize)
-set(get(ax2,'Ylabel'),'String','# of obs : o=poss, \ast=used', ...
+set(get(ax2,'Ylabel'),'String','# of obs : o=possible, \ast=assimilated', ...
'FontSize',figdata.fontsize)
%=====================================================================
@@ -408,13 +470,13 @@
%% 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
+fullname = which(main.fname); % Could be in MatlabPath
if( isempty(fullname) )
- if ( main(1) == '/' ) % must be a absolute pathname
- string1 = sprintf('data file: %s',main);
+ if ( main.fname(1) == '/' ) % must be a absolute pathname
+ string1 = sprintf('data file: %s',main.fname);
else % must be a relative pathname
mydir = pwd;
- string1 = sprintf('data file: %s/%s',mydir,main);
+ string1 = sprintf('data file: %s/%s',mydir,main.fname);
end
else
string1 = sprintf('data file: %s',fullname);
@@ -426,6 +488,15 @@
'Interpreter','none',...
'FontSize',8)
+switch lower(main.trusted)
+ case 'true'
+ h = text(0.0, 1.0,'TRUSTED OBSERVATION. Values include outlying obs. ');
+ set(h,'HorizontalAlignment','center', ...
+ 'VerticalAlignment','middle',...
+ 'Interpreter','none',...
+ 'FontSize',20)
+ otherwise
+end
%=====================================================================
@@ -534,7 +605,7 @@
position = [0.10 0.15 0.8 0.7];
linewidth = 2.0;
-figdata = struct('expcolors', {{'k','r','g','m','b','c','y'}}, ...
+figdata = struct('expcolors', {{'k','r','b','m','g','c','y'}}, ...
'expsymbols', {{'o','s','d','p','h','s','*'}}, ...
'prpolines', {{'-','--'}}, 'position', position, ...
'fontsize',fontsize, 'orientation',orientation, ...
@@ -559,28 +630,6 @@
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
-
-
-
% <next few lines under version control, do not edit>
% $URL$
% $Revision$
Modified: DART/trunk/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_profile.m 2015-02-20 21:57:03 UTC (rev 7604)
+++ DART/trunk/diagnostics/matlab/plot_profile.m 2015-02-20 21:57:32 UTC (rev 7605)
@@ -1,22 +1,51 @@
-function plotdat = plot_profile(fname,copystring)
+function plotdat = plot_profile(fname, copy, varargin)
%% plot_profile plots the vertical profile of the observation-space quantities for all possible levels, all possible variables.
% Part of the observation-space diagnostics routines.
%
% 'obs_diag' produces a netcdf file containing the diagnostics.
+% obs_diag condenses the obs_seq.final information into summaries for a few specified
+% regions - on a level-by-level basis.
%
-% USAGE: plotdat = plot_profile(fname,copystring);
+% The number of observations possible reflects only those observations
+% that have incoming QC values of interest. Any observation with a DART
+% QC of 5 or 6 is not considered 'possible' for the purpose of this graphic.
%
-% fname : netcdf file produced by 'obs_diag'
-% copystring : 'copy' string == quantity of interest. These
-% can be any of the ones available in the netcdf
-% file 'CopyMetaData' variable.
+% NOTE: if the observation was designated as a TRUSTED observation in the
+% obs_diag program, the observations that were rejected by the outlier
+% threshhold STILL PARTICIPATE in the calculation of the rmse, spread, etc.
+% The _values_ plotted by plot_profile reflect that. The number of observations
+% "used" becomes unclear. The number of observations used (designated by the
+% asterisk) is ALWAYS the number of observations successfully assimilated.
+% For TRUSTED observations, this is different than the number used to calculate
+% bias, rmse, spread, etc.
+%
+% USAGE: plotdat = plot_profile(fname, copy);
+%
+% fname : netcdf file produced by 'obs_diag'
+%
+% copy : string defining the metric of interest. 'rmse', 'spread', etc.
+% Possible values are available in the netcdf 'CopyMetaData' variable.
% (ncdump -v CopyMetaData obs_diag_output.nc)
%
-% EXAMPLE:
+% obsname : Optional. If present, The strings of each observation type to plot.
+% Each observation type will be plotted in a separate graphic.
+% Default is to plot all available observation types.
%
-% fname = 'obs_diag_output.nc'; % netcdf file produced by 'obs_diag'
-% copystring = 'totalspread'; % 'copy' string == quantity of interest
-% plotdat = plot_profile(fname,copystring);
+% OUTPUT: 'plotdat' is a structure containing what was plotted.
+% A .pdf of each graphic is created. Each .pdf has a name that
+% reflects the variable, quantity, and region being plotted.
+%
+% EXAMPLE 1: All the observation types possible are plotted in separate figures.
+%
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% plotdat = plot_profile(fname, copy);
+%
+% EXAMPLE 2: Just a single observation type.
+%
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% plotdat = plot_profile(fname, copy, 'obsname', 'RADIOSONDE_TEMPERATURE');
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -28,6 +57,24 @@
% Decode,Parse,Check the input
%---------------------------------------------------------------------
+default_obsname = 'none';
+p = inputParser;
+
+addRequired(p,'fname', at ischar);
+addRequired(p,'copy', at ischar);
+addParamValue(p,'obsname',default_obsname, at ischar);
+parse(p, fname, copy, varargin{:});
+
+% if you want to echo the input
+% disp(['fname : ', p.Results.fname])
+% disp(['copy : ', p.Results.copy])
+% disp(['obsname : ', p.Results.obsname])
+
+if ~isempty(fieldnames(p.Unmatched))
+ disp('Extra inputs:')
+ disp(p.Unmatched)
+end
+
if (exist(fname,'file') ~= 2)
error('file/fname <%s> does not exist',fname)
end
@@ -37,7 +84,7 @@
%---------------------------------------------------------------------
plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.copystring = copy;
plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
@@ -85,13 +132,13 @@
plotdat.toff = plotdat.binedges(1) + iskip;
plotdat.timespan = sprintf('%s through %s', datestr(plotdat.toff), ...
datestr(max(plotdat.binedges(:))));
-plotdat.xlabel = sprintf('%s',copystring);
+plotdat.xlabel = sprintf('%s',copy);
[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindVerticalVars(plotdat);
plotdat.nvars = length(plotdat.varnames);
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copy);
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
plotdat.NQC5index = get_copy_index(fname,'N_DARTqc_5');
@@ -103,14 +150,28 @@
% Loop around (copy-level-region) observation types
%----------------------------------------------------------------------
-for ivar = 1:plotdat.nvars
+% Either use all the variables or just the one optionally specified.
+if strcmp(p.Results.obsname,'none')
+ varlist = 1:plotdat.nvars;
+else
+ varlist = find (strcmpi(p.Results.obsname,plotdat.varnames));
+ if isempty(varlist)
+ error('%s is not in the list of observations',p.Results.obsname)
+ end
+end
+
+for ivar = varlist
+
% create the variable names of interest.
plotdat.myvarname = plotdat.varnames{ivar};
plotdat.guessvar = sprintf('%s_VPguess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_VPanaly',plotdat.varnames{ivar});
+ plotdat.trusted = nc_read_att(fname, plotdat.guessvar, 'TRUSTED');
+ if (isempty(plotdat.trusted)), plotdat.trusted = 'NO'; end
+
% get appropriate vertical coordinate variable
guessdims = nc_var_dims( fname, plotdat.guessvar);
@@ -131,7 +192,7 @@
continue
end
- [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
+ [level_org, level_units, nlevels, level_edges, Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
plotdat.level_org = level_org;
plotdat.level_units = level_units;
plotdat.nlevels = nlevels;
@@ -222,7 +283,7 @@
% plot by region - each in its own figure.
for iregion = 1:plotdat.nregions
- figure(iregion); clf; orient(figuredata.orientation); wysiwyg
+ figure(iregion); clf(iregion); orient(figuredata.orientation); wysiwyg
plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
myplot(plotdat, figuredata);
@@ -251,6 +312,7 @@
g = plotdat.ges_Nposs(:,:,plotdat.region); G = g(plotdat.indices);
a = plotdat.anl_Nposs(:,:,plotdat.region); A = a(plotdat.indices);
+
nobs_poss = G;
nposs_delta = G - A;
@@ -304,6 +366,21 @@
h = legend(h1, str_other_pr, str_other_po, 'Location', 'NorthWest');
set(h,'Interpreter','none','Box','off')
+% If the observation is trusted, reference that somehow
+switch lower(plotdat.trusted)
+ case 'true'
+ axlims = axis;
+ tx = axlims(2) + (axlims(2) - axlims(1))/20;
+ if strcmpi('normal',plotdat.YDir)
+ ty = plotdat.Yrange(1);
+ else
+ ty = plotdat.Yrange(2);
+ end
+ h = text(tx,ty,'TRUSTED. Values include outlying observations.');
+ set(h,'FontSize',20,'Rotation',90,'VerticalAlignment','middle')
+ otherwise
+end
+
% Create another axes to use for plotting the observation counts
ax2 = axes('position',get(ax1,'Position'), ...
@@ -332,7 +409,7 @@
set(get(ax1,'Xlabel'),'String',{plotdat.xlabel, plotdat.timespan}, ...
'Interpreter','none','FontSize',figdata.fontsize)
set(get(ax2,'Xlabel'),'String', ...
- ['# of obs (o=poss, \ast=used) x' int2str(uint32(xscale))],'FontSize',figdata.fontsize)
+ ['# of obs (o=possible, \ast=assimilated) x' int2str(uint32(xscale))],'FontSize',figdata.fontsize)
title({plotdat.myregion, plotdat.myvarname}, ...
'Interpreter', 'none', 'FontSize', figdata.fontsize, 'FontWeight', 'bold')
@@ -400,7 +477,7 @@
%=====================================================================
-function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname)
+function [level_org, level_units, nlevels, level_edges, Yrange] = FindVerticalInfo(fname,varname)
%% Find the vertical dimension and harvest some info
varinfo = nc_getvarinfo(fname,varname);
@@ -571,7 +648,7 @@
position = [0.15 0.12 0.7 0.75];
linewidth = 2.0;
-figdata = struct('expcolors', {{'k','r','g','m','b','c','y'}}, ...
+figdata = struct('expcolors', {{'k','r','b','m','g','c','y'}}, ...
'expsymbols', {{'o','s','d','p','h','s','*'}}, ...
'prpolines', {{'-','--'}}, 'position', position, ...
'fontsize',fontsize, 'orientation',orientation, ...
Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m 2015-02-20 21:57:03 UTC (rev 7604)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m 2015-02-20 21:57:32 UTC (rev 7605)
@@ -1,4 +1,4 @@
-function plotdat = plot_rmse_xxx_evolution(fname, copystring, varargin)
+function plotdat = plot_rmse_xxx_evolution(fname, copy, varargin)
%% plot_rmse_xxx_evolution plots the temporal evolution of the observation-space quantity RMSE and any other for all possible levels, all possible variables.
% Part of the observation-space diagnostics routines.
%
@@ -6,41 +6,64 @@
% obs_diag condenses the obs_seq.final information into summaries for a few specified
% regions - on a level-by-level basis.
%
-% USAGE: plotdat = plot_rmse_xxx_evolution(fname, copystring [,varargin]);
+% The number of observations possible reflects only those observations
+% that have incoming QC values of interest. Any observation with a DART
+% QC of 5 or 6 is not considered 'possible' for the purpose of this graphic.
%
-% fname : netcdf file produced by 'obs_diag'
-% copystring : 'copy' string == quantity of interest. These
-% 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
-% 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.
+% NOTE: if the observation was designated as a TRUSTED observation in the
+% obs_diag program, the observations that were rejected by the outlier
+% threshhold STILL PARTICIPATE in the calculation of the rmse, spread, etc.
+% The _values_ plotted by plot_profile reflect that. The number of observations
+% "used" becomes unclear. The number of observations used (designated by the
+% asterisk) is ALWAYS the number of observations successfully assimilated.
+% For TRUSTED observations, this is different than the number used to calculate
+% bias, rmse, spread, etc.
%
-% OUTPUT: two files will result for each observation type plotted. One is a
-% postscript file containing a page for each level - all regions.
+% USAGE: plotdat = plot_evolution(fname, copy);
+%
+% fname : netcdf file produced by 'obs_diag'
+%
+% copy : string defining the metric of interest. 'rmse', 'spread', etc.
+% Possible values are available in the netcdf 'CopyMetaData' variable.
+% (ncdump -v CopyMetaData obs_diag_output.nc)%
+%
+% obsname : Optional. If present, The strings of each observation type to plot.
+% Each observation type will be plotted in a separate graphic.
+% Default is to plot all available observation types.
+%
+% level : Optional. 'level' index. Default is to plot all levels.
+%
+% range : Optional. 'range' of the value being plotted. Default is to
+% automatically determine range based on the data values.
+%
+% OUTPUT: 'plotdat' is a structure containing what was last plotted.
+% A postscript file containing a page for each level - each region.
% 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.
+% Both of these filenames contain the observation type,
+% copy and region 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'
-% copystring = 'totalspread'; % 'copy' string == quantity of interest
-% plotdat = plot_rmse_xxx_evolution(fname,copystring);
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% plotdat = plot_rmse_xxx_evolution(fname, copy);
%
%
% EXAMPLE 2 - plot the RMSE and spread on the same axis - for just one observation type.
%
-% fname = 'obs_diag_output.nc';
-% copystring = 'totalspread';
-% obstype = 'RADIOSONDE_TEMPERATURE';
-% plotdat = plot_rmse_xxx_evolution(fname, copystring, obstype);
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% obsname = 'RADIOSONDE_TEMPERATURE';
+% plotdat = plot_rmse_xxx_evolution(fname, copy, 'obsname', obsname);
+%
+%
+% EXAMPLE 3 - plot the RMSE and spread on the same axis - for just one observation type, 1 level.
+%
+% fname = 'obs_diag_output.nc';
+% copy = 'totalspread';
+% bob = 'RADIOSONDE_TEMPERATURE';
+% plotdat = plot_rmse_xxx_evolution(fname,copy,'obsname',bob,'level',3,'range',[-1 5]);
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -48,13 +71,39 @@
%
% DART $Id$
-if nargin == 2
+default_level = -1;
+default_obsname = 'none';
+default_range = [NaN NaN];
+p = inputParser;
+
+addRequired(p,'fname', at ischar);
+addRequired(p,'copy', at ischar);
+addParamValue(p,'obsname',default_obsname, at ischar);
+addParamValue(p,'range',default_range, at isnumeric);
+addParamValue(p,'level',default_level, at isnumeric);
+parse(p, fname, copy, varargin{:});
+
+% if you want to echo the input
+% fprintf('fname : %s\n', p.Results.fname)
+% fprintf('copy : %s\n', p.Results.copy)
+% fprintf('obsname : %s\n', p.Results.obsname)
+% fprintf('level : %d\n', p.Results.level)
+% fprintf('range : %f %f \n', p.Results.range)
+
+if ~isempty(fieldnames(p.Unmatched))
+ disp('Extra inputs:')
+ disp(p.Unmatched)
+end
+
+if (numel(p.Results.range) ~= 2)
+ error('range must be an array of length two ... [bottom top]')
+end
+
+if strcmp(p.Results.obsname,'none')
nvars = 0;
-elseif nargin == 3
- varname = varargin{1};
+else
+ obsname = p.Results.obsname;
nvars = 1;
-else
- error('wrong number of arguments ... ')
end
if (exist(fname,'file') ~= 2)
@@ -66,7 +115,7 @@
%---------------------------------------------------------------------
plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.copystring = copy;
plotdat.bincenters = nc_varget(fname,'time');
plotdat.binedges = nc_varget(fname,'time_bounds');
plotdat.mlevel = local_nc_varget(fname,'mlevel');
@@ -84,15 +133,15 @@
plotdat.region_names = deblank(plotdat.region_names');
end
-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');
+dimensionality = nc_read_att(fname, nc_global, 'LocationRank');
+plotdat.binseparation = nc_read_att(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_read_att(fname, nc_global, 'bin_width');
+time_to_skip = nc_read_att(fname, nc_global, 'time_to_skip');
+plotdat.lonlim1 = nc_read_att(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_read_att(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_read_att(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_read_att(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_read_att(fname, nc_global, 'bias_convention');
% Coordinate between time types and dates
@@ -121,11 +170,11 @@
[plotdat.varnames, plotdat.vardims] = FindTemporalVars(plotdat);
plotdat.nvars = length(plotdat.varnames);
else
- plotdat.varnames{1} = varname;
+ plotdat.varnames{1} = obsname;
plotdat.nvars = nvars;
end
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copy);
plotdat.rmseindex = get_copy_index(fname,'rmse');
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
@@ -148,6 +197,9 @@
plotdat.guessvar = sprintf('%s_guess',plotdat.varnames{ivar});
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list