function two_experiments_profile(files, titles, obsnames, copy, prpo, varargin) % Plot two or more experiments on the same axis. % 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. % % 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. % % 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: two_experiments_profile(files, titles, obsnames, copy, prpo, 'level', 1) % % files : Cell array containing the locations of the obs_diag_output.nc % files to compare. Each file is presumed to be the results from % a single experiment. % % titles : Cell array containing the titles used to annotate each of the experiments. % % obsnames : Cell array containing the strings of each observation type to plot. % Each observation type will be plotted in a separate graphic. % % 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) % % prpo : string defining whether to plot the prior or posterior metrics. % Due to the amount of information already plotted, we made a % conscious decision not to support plotting both prior and posterior % on the same plot. % % varargin: optional parameter-value pairs. Supported parameters are described below. % % level : The index of the level to plot. Defaults to level 1. % % verbose : true/false to control amount of run-time output % % MarkerSize : integer controlling the size of the symbols % % pause : true/false to conrol pausing after each figure is created. % true will require hitting any key to continue to next plot % % range : 'range' of the value being plotted. Default is to % automatically determine range based on the data values. % % OUTPUT: A .pdf of each graphic is created. Each .pdf has a name that % reflects the variable, quantity, and region being plotted. % % EXAMPLE % % files = {'/ptmp/nancy/CSL/Base5/032-061s0_def_reg/obs_diag_output.nc', % '/ptmp/thoar/GPS+AIRS/Sep_032-061/obs_diag_output.nc'}; % titles = {'Base5', 'GPS+AIRS'}; % files = {'/glade/scratch/nancy/fvdiags/obs_diag_2005_08.nc', ... % '/glade/scratch/raeder/SE_NCEP_assim1/Diag_hemi_poles_2005.8.1-30/obs_diag_output.nc'}; % titles = {'FV', 'SE'}; % obsnames = {'RADIOSONDE_TEMPERATURE', ... % 'RADIOSONDE_U_WIND_COMPONENT', 'RADIOSONDE_V_WIND_COMPONENT'}; % copy = 'rmse'; % rmse, spread, totalspread, bias, etc. % prpo = 'analysis'; % [analy, analysis, posterior ] == posterior % prpo = 'forecast'; % [guess, forecast, prior ] == prior % prpo = 'both'; % % two_experiments_profile(files, titles, obsnames, copy, prpo) % % Example 2: restrict the data limits to the data in a certain vertical area. % In this case, the observations using a pressure vertical coordinate % between +Inf (the surface) and 100hPa (inclusive) are used to determine % the scale. All values will be plotted, but the highest levels may be % clipped. The optional argument pairs at the end consist of a string % and a length 2 array specifying the [bottom top] levels to consider. % % two_experiments_profile(files, titles, obsnames, copy, prpo,'plevel',[Inf 100]) % % two_experiments_profile(files, titles, obsnames, copy, prpo, ... % 'plevel',[Inf 100],'mlevel',[1 10],'hlevel',[0 20000]) %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at % http://www.image.ucar.edu/DAReS/DART/DART_download % % DART $Id$ %%-------------------------------------------------------------------- % Decode,Parse,Check the input %--------------------------------------------------------------------- defaultPlevels = [ Inf 0 ]; defaultHlevels = [-Inf Inf]; defaultMlevels = [ 1 Inf]; default_verbosity = true; default_markersize = 12; default_pause = false; default_range = [NaN NaN]; default_level = 1; p = inputParser; addRequired(p,'files',@iscell); addRequired(p,'titles',@iscell); addRequired(p,'obsnames',@iscell); addRequired(p,'copy',@ischar); addRequired(p,'prpo',@ischar); if (exist('inputParser/addParameter','file') == 2) addParameter(p,'plevel',defaultPlevels,@isnumeric); addParameter(p,'hlevel',defaultHlevels,@isnumeric); addParameter(p,'mlevel',defaultMlevels,@isnumeric); addParameter(p,'verbose', default_verbosity, @islogical); addParameter(p,'MarkerSize', default_markersize, @isnumeric); addParameter(p,'pause', default_pause, @islogical); addParameter(p,'range', default_range, @isnumeric); addParameter(p,'level', default_level, @isnumeric); else addParamValue(p,'plevel',defaultPlevels,@isnumeric); %#ok addParamValue(p,'hlevel',defaultHlevels,@isnumeric); %#ok addParamValue(p,'mlevel',defaultMlevels,@isnumeric); %#ok addParamValue(p,'verbose', default_verbosity, @islogical); %#ok addParamValue(p,'MarkerSize',default_markersize, @isnumeric); %#ok addParamValue(p,'pause', default_pause, @islogical); %#ok addParamValue(p,'range', default_range, @isnumeric); %#ok addParamValue(p,'level', default_level, @isnumeric); %#ok end p.parse(files, titles, obsnames, copy, prpo, varargin{:}); % if you want to echo the input % disp(['files : ', p.Results.files]) % disp(['titles : ', p.Results.titles]) % disp(['obsnames: ', p.Results.obsnames]) % disp(['copy : ', p.Results.copy]) % disp(['prpo : ', p.Results.prpo]) % fprintf( 'plevel : %f %f \n', p.Results.plevel) % fprintf( 'hlevel : %f %f \n', p.Results.hlevel) % fprintf( 'mlevel : %f %f \n', p.Results.mlevel) if ~isempty(fieldnames(p.Unmatched)) disp('Extra inputs:') disp(p.Unmatched) end % if ~isempty(p.UsingDefaults) % disp('Using defaults: ') % disp(p.UsingDefaults) % end if (numel(p.Results.plevel) ~= 2) error('plevel must be an array of length two ... [bottom top]') end if (numel(p.Results.hlevel) ~= 2) error('hlevel must be an array of length two ... [bottom top]') end if (numel(p.Results.mlevel) ~= 2) error('mlevel must be an array of length two ... [bottom top]') end % Now that the input passes sanity checks ... if (length(files) ~= length(titles)) error('each file must have a title') end NumExp = length(files); for i = 1:NumExp if (exist(files{i},'file') ~= 2) error('File %s does not exist',files{i}) end end if (NumExp ~= length(titles)) error('each file must have an experiment title') end %% set up all the stuff that is common. global figuredata commondata = check_compatibility(files, prpo, obsnames, copy); figuredata = set_obsdiag_figure('tall', 'numexp', NumExp); figuredata.MarkerSize = p.Results.MarkerSize; %%-------------------------------------------------------------------- % Set some static data %--------------------------------------------------------------------- nvars = length(obsnames); for ivar = 1:nvars fprintf('Working on %s ...\n',obsnames{ivar}) for iregion = 1:commondata.nregions %--------------------------------------------------------------------- % Getting the data for each experiment %--------------------------------------------------------------------- Dlimits = zeros(NumExp,2); % range of the data plotobj = cell(1,NumExp); for iexp = 1:NumExp plotobj{iexp} = getvals(files{iexp}, commondata.targets{ivar}, copy, iregion, p); plotobj{iexp}.title = titles{iexp}; plotobj{iexp}.phase = commondata.phase; Dlimits(iexp,:) = plotobj{iexp}.Drange; end %--------------------------------------------------------------------- % Find nice limits that encompass all experiments % Note that Dlimits has been constructed by ignoring the top levels. %--------------------------------------------------------------------- Drange = [min(Dlimits(:,1)) max(Dlimits(:,2))]; span = abs(Drange(2) - Drange(1))* 0.05; Drange(1) = Drange(1) - span; Drange(2) = Drange(2) + span; %--------------------------------------------------------------------- % Plot all regions - one region to a page %--------------------------------------------------------------------- myplot(plotobj, Drange); BottomAnnotation(files) psfname = sprintf('%s_%s_region%d_profile_%dexp', ... obsnames{ivar}, plotobj{1}.copystring, iregion, NumExp); if verLessThan('matlab','R2016a') print(iregion, '-dpdf', psfname) else print(iregion, '-dpdf', '-bestfit', psfname) end end % of loop around regions if ( ivar ~= nvars && p.Results.pause ) disp('Pausing, hit any key to continue ...') pause end end % of loop around variable %===================================================================== % End of main function body. Helper functions below. %===================================================================== function common = check_compatibility(filenames, prpo, varnames, copystring) %% Trying to prevent the comparison of apples and oranges. % make sure the diagnostics were generated the same way. % need to check that the timeframe is the same for all files % need to check that the region definitions are the same for all files mystat = 0; nexp = length(filenames); commondata = cell(1,nexp); targets = struct([]); for i = 1:length(varnames) switch lower(prpo) case {'guess','forecast','prior'} targets{i} = sprintf('%s_VPguess',varnames{i}); commondata{i}.phase = 'prior'; case {'analy','analysis','posterior'} targets{i} = sprintf('%s_VPanaly',varnames{i}); commondata{i}.phase = 'posterior'; otherwise error('unknown prpo ... "%s"',prpo) end end for i = 1:nexp varexist(filenames{i}, {targets{:}, 'time', 'time_bounds'}) %#ok commondata{i}.targets = targets; commondata{i}.region_names = strtrim(ncread(filenames{i},'region_names')'); commondata{i}.times = ncread(filenames{i}, 'time'); commondata{i}.time_bnds = ncread(filenames{i}, 'time_bounds'); commondata{i}.copyindex = get_copy_index(filenames{i},copystring); commondata{i}.nobstypes = nc_dim_info(filenames{i}, 'obstypes'); commondata{i}.nregions = nc_dim_info(filenames{i}, 'region'); commondata{i}.time_to_skip = nc_read_att(filenames{i}, '/', 'time_to_skip'); commondata{i}.lonlim1 = nc_read_att(filenames{i}, '/', 'lonlim1'); commondata{i}.lonlim2 = nc_read_att(filenames{i}, '/', 'lonlim2'); commondata{i}.latlim1 = nc_read_att(filenames{i}, '/', 'latlim1'); commondata{i}.latlim2 = nc_read_att(filenames{i}, '/', 'latlim2'); end % error checking - compare everything to the first experiment for i = 2:nexp if (any(commondata{i}.lonlim1 ~= commondata{1}.lonlim1)) fprintf('The left longitudes of the regions (i.e. lonlim1) are not compatible.\n') mystat = 1; end if (any(commondata{i}.lonlim2 ~= commondata{1}.lonlim2)) fprintf('The right longitudes of the regions (i.e. lonlim2) are not compatible.\n') mystat = 1; end if (any(commondata{i}.latlim1 ~= commondata{1}.latlim1)) fprintf('The bottom latitudes of the regions (i.e. latlim1) are not compatible.\n') mystat = 1; end if (any(commondata{i}.latlim2 ~= commondata{1}.latlim2)) fprintf('The top latitudes of the regions (i.e. latlim2) are not compatible.\n') mystat = 1; end if (any(commondata{i}.time_bnds ~= commondata{1}.time_bnds)) fprintf('The time boundaries of the experiments (i.e. time_bnds) are not compatible.\n') mystat = 1; end if (any(commondata{i}.time_to_skip ~= commondata{1}.time_to_skip)) fprintf('The time skipped in the experiments (i.e. time_to_skip) are not compatible.\n') mystat = 1; end end if mystat > 0 error('The experiments are not compatible ... stopping.') end common = commondata{1}; % Coordinate between time types and dates timeunits = nc_read_att(filenames{1},'time','units'); timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD timeorigin = datenum(timebase(1),timebase(2),timebase(3)); timefloats = zeros(size(commondata{1}.time_to_skip)); % stupid int32 type conversion timefloats(:) = commondata{1}.time_to_skip(:); skip_seconds = timefloats(4)*3600 + timefloats(5)*60 + timefloats(6); iskip = timefloats(3) + skip_seconds/86400; common.bincenters = commondata{1}.times + timeorigin; common.binedges = commondata{1}.time_bnds + timeorigin; common.Nbins = length(common.bincenters); common.toff = common.binedges(1) + iskip; common.timespan = sprintf('%s through %s', datestr(common.toff), ... datestr(max(common.binedges(:)))); %===================================================================== function plotdat = getvals(fname, varname, copystring, regionindex, opt ) %% basic function to retrieve plotting data if (exist(fname,'file') ~= 2) error('%s does not exist',fname) end plotdat = read_obsdiag_staticdata(fname,copystring); plotdat.region = regionindex; plotdat.varname = varname; myinfo.diagn_file = fname; myinfo.copyindex = plotdat.copyindex; myinfo.regionindex = plotdat.region; [start, count] = GetNCindices(myinfo,'diagn',plotdat.varname); hyperslab = ncread(fname, plotdat.varname, start, count); plotdat.data = squeeze(hyperslab); plotdat.trusted = nc_read_att(fname, plotdat.varname, 'TRUSTED'); if (isempty(plotdat.trusted)), plotdat.trusted = 'NO'; end % Now that we know the variable ... get the appropriate vertical information varinfo = ncinfo(fname,plotdat.varname); plotdat.levels = ncread(fname,varinfo.Dimensions(2).Name); plotdat.level_units = nc_read_att(fname,varinfo.Dimensions(2).Name,'units'); plotdat.nlevels = length(plotdat.levels); plotdat.level_edges = ncread(fname,sprintf('%s_edges',varinfo.Dimensions(2).Name)); plotdat.YDir = 'normal'; inds = 1:plotdat.nlevels; %% Determine data limits % find the levels of interest for setting the data limits switch lower(varinfo.Dimensions(2).Name) case {'plevel'} plotdat.YDir = 'reverse'; inds = find((plotdat.levels <= opt.Results.plevel(1)) & ... (plotdat.levels >= opt.Results.plevel(2))); case {'hlevel'} inds = find((plotdat.levels >= opt.Results.hlevel(1)) & ... (plotdat.levels <= opt.Results.hlevel(2))); case {'mlevel'} inds = find((plotdat.levels >= opt.Results.mlevel(1)) & ... (plotdat.levels <= opt.Results.mlevel(2))); otherwise end bob = plotdat.data(inds); switch copystring case {'bias'} % always make sure we have a zero bias line ... dmin = min( [ min(bob) 0.0 ] ); dmax = max( [ max(bob) 0.0 ] ); plotdat.Drange = [ dmin dmax ]; plotdat.xlabel = sprintf('%s (%s)',copystring, plotdat.biasconv); case {'rmse'} plotdat.Drange = [0.0 max(bob)]; plotdat.xlabel = copystring; otherwise plotdat.Drange = [min(bob) max(bob)]; plotdat.xlabel = copystring; end qcvalues = get_qc_values(fname, plotdat.varname, ... 'regionindex', plotdat.region, ... 'fatal', false, ... 'verbose', false); plotdat.nposs = qcvalues.nposs; plotdat.nused = qcvalues.nused; plotdat.num_evaluated = qcvalues.num_evaluated; if sum(plotdat.num_evaluated(:) > 0) plotdat.assim_eval_string = 'evaluated'; else plotdat.assim_eval_string = 'assimilated'; end %===================================================================== function myplot( plotdat, Drange) %% Create graphic for one region - for all experiments. global figuredata Nexp = length(plotdat); iregion = plotdat{1}.region; figure(iregion); clf(iregion); orient(figuredata.orientation); ax1 = subplot('position',figuredata.position); Stripes(Drange, plotdat{1}.level_edges, plotdat{1}.level_units, Nexp); set(ax1,'YDir',plotdat{1}.YDir,'YTick',sort(plotdat{1}.levels),'Layer','top') set(ax1,'YAxisLocation','left','FontSize',figuredata.fontsize) % draw the results of the experiments - each with their own line type. hd = []; % handle to an unknown number of data lines legstr = {[]}; % strings for the legend for i = 1:Nexp hd(i) = line(plotdat{i}.data, plotdat{i}.levels, ... 'Color', figuredata.expcolors{i}, ... 'Marker', figuredata.expsymbols{i}, ... 'MarkerSize', figuredata.MarkerSize, ... 'MarkerFaceColor', figuredata.expcolors{i}, ... 'LineStyle', figuredata.prpolines{1}, ... 'LineWidth', figuredata.linewidth,'Parent',ax1); %#ok legstr{i} = sprintf('%s %s',plotdat{i}.title,plotdat{i}.phase); end switch plotdat{1}.copystring case {'bias','rmse'} zeroline = line([0 0],get(ax1,'YLim'),'Color',[200 200 200]/255,'Parent',ax1); set(zeroline,'LineWidth',2.5,'LineStyle','-') otherwise end set(ax1,'XLim',Drange) axlims = axis; dx = (axlims(2) - axlims(1))/20; if strcmpi('normal',get(ax1,'YDir')) ty = axlims(3); else ty = axlims(4); end for i = 1:Nexp % If the observation is trusted, reference that somehow switch lower(plotdat{i}.trusted) case 'true' tx = axlims(2) + (i*dx); h = text(tx,ty,sprintf('TRUSTED OBSERVATION in %s',plotdat{i}.title)); set(h, 'FontSize', 20, 'Rotation', 90, ... 'VerticalAlignment', 'middle', 'Interpreter', 'none') otherwise end end % Create another axes to use for plotting the observation counts ax2 = axes('position',get(ax1,'Position'), ... 'XAxisLocation','top', ... 'YAxisLocation','right',... 'Color','none', ... 'XColor','b', ... 'YColor',get(ax1,'YColor'), ... 'YLim',get(ax1,'YLim'), ... 'YDir',get(ax1,'YDir'), ... 'FontSize',get(ax1,'FontSize')); % Plot the data, which sets the range of the axis for i = 1:Nexp ax2h1 = line(plotdat{i}.nposs, plotdat{i}.levels, 'Parent', ax2); ax2h2 = line(plotdat{i}.nused, plotdat{i}.levels, 'Parent', ax2); set(ax2h1,'LineStyle','none', ... 'Color', figuredata.expcolors{i}, ... 'Marker', figuredata.obs_marker, ... 'MarkerSize',figuredata.MarkerSize); set(ax2h2,'LineStyle','none', ... 'Color', figuredata.expcolors{i}, ... 'Marker', figuredata.ges_marker, ... 'MarkerSize',figuredata.MarkerSize); end % use same Y ticks but no labels set(ax2,'YTick',get(ax1,'YTick'), 'YTicklabel',[]); % use the same X ticks, but find the right label values xscale = matchingXticks(ax1,ax2); % Annotate the whole thing - gets pretty complicated for multiple % regions on one page. Trying to maximize content, minimize clutter. % Any plot object will do for annotating region,levels,etc annotate(ax1, ax2, plotdat{1}, xscale) lh = legend(hd,legstr,'Location','NorthWest'); set(lh,'Interpreter','none','Box','off','FontSize',figuredata.fontsize); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. % Before 2017a, this was the default behavior, so do nothing. % We do not want to add the bias line to the legend, for example. else lh.AutoUpdate = 'off'; end %===================================================================== function h = Stripes(x, edges, units, nexp) %% plot the subtle background stripes that indicate the vertical % extent of what was averaged. % % axlims(3) should be conditional on the observation vertical coordinate: % values for pressure coordinates are inappropriate for height coord. % plot two little dots at the corners to make Matlab % choose some plot limits. Given those nice limits and % tick labels ... KEEP THEM. Later, make the dots invisible. h = plot([min(x) max(x)],[min(edges) max(edges)]); axlims = axis; legend_fraction = 0.05 * nexp + 0.005; % partial fix to legend space; add in option for vert coord = height. switch lower(units) case 'hpa' axlims(4) = max(edges); axlims(3) = min(edges) - legend_fraction*(axlims(4)-min(edges)); case 'm' axlims(3) = min(edges) ; axlims(4) = max(edges) + legend_fraction*(max(edges)-axlims(3)); otherwise end axis(axlims) % set up list of x,y values defining corner of every other stripe xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ]; hold on; for i = 1:2:(length(edges)-1) yc = [ edges(i) edges(i) edges(i+1) edges(i+1) edges(i) ]; fill(xc,yc,[0.8 0.8 0.8],'EdgeColor','none'); end hold off; set(gca,'XGrid','on') set(h,'Visible','off') % make the dots invisible %===================================================================== function annotate(ax1, ax2, plotobj, xscale) %% One figure ... everything gets annotated. global figuredata set(get(ax1,'Ylabel'),'String',plotobj.level_units, ... 'Interpreter','none','FontSize',figuredata.fontsize) set(get(ax1,'Xlabel'),'String',{plotobj.xlabel,plotobj.timespan}, ... 'Interpreter','none','FontSize',figuredata.fontsize) string1 = sprintf('# of obs (o=possible, %s=%s ) x ', '\ast', ... plotobj.assim_eval_string); set(get(ax2,'Xlabel'),'String', [string1 int2str(uint32(xscale))], ... 'FontSize',figuredata.fontsize) th = title({deblank(plotobj.region_names(plotobj.region,:)), plotobj.varname}); set(th,'Interpreter','none','FontSize',figuredata.fontsize,'FontWeight','bold'); %===================================================================== function BottomAnnotation(filenames) %% annotates the filenames containing the data being plotted nfiles = length(filenames); dy = 1.0/(nfiles+3); yheight = 0.0225*(nfiles+1); subplot('position',[0.10 0.01 0.8 yheight]) axis off % list all the files. for ifile = 1:nfiles main = filenames{ifile}; fullname = which(main); % Could be in MatlabPath if( isempty(fullname) ) if ( main(1) == '/' ) % must be a absolute pathname string1 = sprintf('data file: %s',main); else % must be a relative pathname mydir = pwd; string1 = sprintf('data file: %s/%s',mydir,main); end else string1 = sprintf('data file: %s',fullname); end ty = 1.0 - (ifile+2)*dy; h = text(0.5, ty, string1); set(h, 'Interpreter', 'none', 'FontSize', 8); set(h, 'HorizontalAlignment','center'); end %===================================================================== function varexist(filename, varnames) %% We already know the file exists by this point. % Lets check to make sure that file contains all needed variables. nvars = length(varnames); gotone = ones(1,nvars); for i = 1:nvars gotone(i) = nc_var_exists(filename,varnames{i}); if ( ~ gotone(i) ) fprintf('\n%s is not a variable in %s\n',varnames{i},filename) end end if ~ all(gotone) error('missing required variable ... exiting') end % % $URL$ % $Revision$ % $Date$