Home > atmlab > graphs > scatter_density_plot.m

scatter_density_plot

PURPOSE ^

% scatter_density_plot Create scatter density plot in style of Eliasson et. al (2013)

SYNOPSIS ^

function [h_ax, h_cb] = scatter_density_plot(x, y, varargin)

DESCRIPTION ^

% scatter_density_plot Create scatter density plot in style of Eliasson et. al (2013)

 Create a scatter-density plot in the style of Eliasson et. al (2013).
 This routine plots a 2-dimensional histogram between the provided
 variables, optionally plotting median y(x), median x(y) and the actual
 measurement points in the same graph.

 Adding labels, titles, axes-properties etc. is not done in this function,
 but should be done by the user manually, or by passing the 'axprops'
 option. See below for an example.

 IN

   x       vector  data on x-axis
   y       vector  data on y-axis
   opts    structure with options:
       bins        no. bins or vector with all bins (see hist2d)
       xbins       bins for x only (default: equal to bins)
       ybins       bins for y only (default: equal to bins)
       trans       transformation, e.g. @(x)log10(x) on both axis
       invtrans    inverse of trans, e.g. @(x)10.^(x) on both axis
       transx       transformation, e.g. @(x)log10(x)
       invtransx    inverse of trans, e.g. @(x)10.^(x)
       axprops     structure with axes properties
       ncolour     number of unique colours
       scatprops   if defined, structure with props for scatter-points
       diagonal    if defined, structure with props for diag-line
       medprops    if defined, structure or array of two structures
                   with props for median-lines.  If those are two
                   structures, the first is used for medians as a
                   function of x, the second for medians as a function
                   of y.
                   Additional fields that can be included in opts.medprops:
                   'atleastN'    Only calculate median where there are
                                 atleast N binned values. 
                   legend.names  = {names} to be displayed in a legend related to the median lines
              e.g. legend.Location = 'NorthWest'
                   etc. (except for legend.names, fieldnames and options must be valid for legend() )
       linreg      if defined, structure with props for linreg-line
       normalise   string representing how to normalise, if at all:
                   'none'  (default) just write collocs
                   'x'     normalise per x-bin
                   
       ...
   
 OUT

   h_ax     handle to main axes
   h_cb     handle to colourbar-axes

 EXAMPLE

   To compare IWP products with logarithmic axes, one may use e.g.:

   opts.trans = @(x)log10(x);   % provide transformation on both axis
   opts.invtrans = @(x)10.^(x); % inverse transformation on both axis
   opts.transx = @(x)log10(x);   % provide transformation
   opts.invtransx = @(x)10.^(x); % inverse transformation
   opts.axprops.xscale = 'log'; % logarithmic x-scale
   opts.axprops.yscale = 'log'; % idem dito
   opts.medprops = struct('LineWidth', 3, 'LineStyle', '-', 'Color', [0 0 0]);
   [h_ax, h_cb] = scatter_density_plot(x, y, opts);

 $Id: scatter_density_plot.m 8777 2014-02-12 12:45:13Z gerrit $
 Created by Gerrit Holl based on code from Salomon Eliasson

 See also: hist, hist2d, bin, median

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

scatter_density_plot.m

SOURCE CODE ^

0001 function [h_ax, h_cb] = scatter_density_plot(x, y, varargin)
0002 
0003 %% scatter_density_plot Create scatter density plot in style of Eliasson et. al (2013)
0004 %
0005 % Create a scatter-density plot in the style of Eliasson et. al (2013).
0006 % This routine plots a 2-dimensional histogram between the provided
0007 % variables, optionally plotting median y(x), median x(y) and the actual
0008 % measurement points in the same graph.
0009 %
0010 % Adding labels, titles, axes-properties etc. is not done in this function,
0011 % but should be done by the user manually, or by passing the 'axprops'
0012 % option. See below for an example.
0013 %
0014 % IN
0015 %
0016 %   x       vector  data on x-axis
0017 %   y       vector  data on y-axis
0018 %   opts    structure with options:
0019 %       bins        no. bins or vector with all bins (see hist2d)
0020 %       xbins       bins for x only (default: equal to bins)
0021 %       ybins       bins for y only (default: equal to bins)
0022 %       trans       transformation, e.g. @(x)log10(x) on both axis
0023 %       invtrans    inverse of trans, e.g. @(x)10.^(x) on both axis
0024 %       transx       transformation, e.g. @(x)log10(x)
0025 %       invtransx    inverse of trans, e.g. @(x)10.^(x)
0026 %       axprops     structure with axes properties
0027 %       ncolour     number of unique colours
0028 %       scatprops   if defined, structure with props for scatter-points
0029 %       diagonal    if defined, structure with props for diag-line
0030 %       medprops    if defined, structure or array of two structures
0031 %                   with props for median-lines.  If those are two
0032 %                   structures, the first is used for medians as a
0033 %                   function of x, the second for medians as a function
0034 %                   of y.
0035 %                   Additional fields that can be included in opts.medprops:
0036 %                   'atleastN'    Only calculate median where there are
0037 %                                 atleast N binned values.
0038 %                   legend.names  = {names} to be displayed in a legend related to the median lines
0039 %              e.g. legend.Location = 'NorthWest'
0040 %                   etc. (except for legend.names, fieldnames and options must be valid for legend() )
0041 %       linreg      if defined, structure with props for linreg-line
0042 %       normalise   string representing how to normalise, if at all:
0043 %                   'none'  (default) just write collocs
0044 %                   'x'     normalise per x-bin
0045 %
0046 %       ...
0047 %
0048 % OUT
0049 %
0050 %   h_ax     handle to main axes
0051 %   h_cb     handle to colourbar-axes
0052 %
0053 % EXAMPLE
0054 %
0055 %   To compare IWP products with logarithmic axes, one may use e.g.:
0056 %
0057 %   opts.trans = @(x)log10(x);   % provide transformation on both axis
0058 %   opts.invtrans = @(x)10.^(x); % inverse transformation on both axis
0059 %   opts.transx = @(x)log10(x);   % provide transformation
0060 %   opts.invtransx = @(x)10.^(x); % inverse transformation
0061 %   opts.axprops.xscale = 'log'; % logarithmic x-scale
0062 %   opts.axprops.yscale = 'log'; % idem dito
0063 %   opts.medprops = struct('LineWidth', 3, 'LineStyle', '-', 'Color', [0 0 0]);
0064 %   [h_ax, h_cb] = scatter_density_plot(x, y, opts);
0065 %
0066 % $Id: scatter_density_plot.m 8777 2014-02-12 12:45:13Z gerrit $
0067 % Created by Gerrit Holl based on code from Salomon Eliasson
0068 %
0069 % See also: hist, hist2d, bin, median
0070 
0071 % set default values and handle options
0072 defaults = struct(...
0073     'bins', 50, ...
0074     'xbins', [], ...
0075     'ybins', [], ...
0076     'trans', @(x)x, ...
0077     'transx',[],...
0078     'transy',[],...
0079     'invtrans', @(x)x, ...
0080     'invtransx', [], ...
0081     'invtransy', [], ...
0082     'axprops', struct(), ...
0083     'ncolour', 20, ...
0084     'medprops', [], ...
0085     'scatprops', [], ...
0086     'diagonal', [], ...
0087     'linreg', [], ...
0088     'normalise', '(none)');
0089 opts = optargs(varargin, {defaults});
0090 
0091 opts = optargs_struct(opts, defaults);
0092 
0093 if isempty(opts.xbins)
0094     opts.xbins = opts.bins;
0095 end
0096 
0097 if isempty(opts.ybins)
0098     opts.ybins = opts.bins;
0099 end
0100 
0101 if isempty(opts.transx)
0102    opts.transx = opts.trans;
0103    opts.invtransx = opts.invtrans;
0104 end
0105 
0106 if isempty(opts.transy)
0107    opts.transy = opts.trans;
0108    opts.invtransy = opts.invtrans;
0109 end
0110 
0111 
0112 % create the actual 2-d histogram. Apply forward transformation upon data,
0113 % then inverse transformation upon plotting. This allows for, e.g.
0114 % logarithmic representations.
0115 [N, BX, BY] = hist2d(opts.transx(x), opts.transy(y), opts.xbins, opts.ybins);
0116 xbins = BX(1, :);
0117 ybins = BY(:, 1);
0118 
0119 paddedx = [xbins(1)-diff(xbins(1:2)),xbins,xbins(end)+diff(xbins(1:2))];
0120 paddedy = [ybins(1)-diff(ybins(1:2));ybins;ybins(end)+diff(ybins(1:2))];
0121 paddedN = [zeros(1,size(N,2)+2);[zeros(size(N,1),1),N,zeros(size(N,1),1)];zeros(1,size(N,2)+2)];
0122 
0123 switch opts.normalise
0124     case '(none)' % do nothing
0125     case 'x'
0126         newpaddedN = bsxfun(@rdivide, paddedN, sum(paddedN, 1));
0127         newpaddedN(:, sum(paddedN, 1)==0) = 0;
0128         paddedN = newpaddedN;
0129     otherwise
0130         error(['atmlab:' mfilename ':invalid', ...
0131             'Invalid option for normalise: %s. See help.'], opts.normalise);
0132 end
0133 
0134 %% plot the actual 2d-histogram
0135 pcolor(opts.invtransx(paddedx), opts.invtransy(paddedy), paddedN);
0136 %pcolor(opts.invtransx(xbins), opts.invtransy(ybins), N);
0137 
0138 set(gca(), opts.axprops);
0139 shading('flat');
0140 
0141 % use a custom colour-table, most notably with white as the base colour
0142 largecolour = getColourTable(length(unique(N)),opts);
0143 set(gcf(), 'colormap', largecolour);
0144 
0145 %% plot medians
0146 
0147 %plot_quantiles();
0148 hold('on');
0149 if ~isempty(opts.medprops)
0150     plot_medians(x, y, paddedx, paddedy, opts);
0151 end
0152 
0153 %% plot actual measurements
0154 
0155 if ~isempty(opts.scatprops)
0156     plot(x, y, opts.scatprops);
0157 end
0158 
0159 %% plot diagonal
0160 
0161 if ~isempty(opts.diagonal)
0162     plot([1e-3 1e6], [1e-3 1e6], opts.diagonal);
0163 end
0164 %plot_medians(y, x, ybins, opts);
0165 
0166 %% plot linreg
0167 
0168 if ~isempty(opts.linreg)
0169     p = polyfit(x, y, 1);
0170     xx = linspace(min(x), max(x), 10);
0171     yy = polyval(p, xx);
0172     plot(xx, yy, opts.linreg);
0173 end
0174 
0175 
0176 
0177 %% tweak the looks
0178 
0179 axis('square');
0180 grid('on');
0181 set(gca(), 'layer', 'top'); % otherwise grid is below everything
0182 
0183 h_ax = gca();
0184 h_cb = colorbar();
0185 
0186 hold('off');
0187 end
0188 
0189 function plot_medians(x, y, xbins, ybins, opts)
0190 
0191 binned_x = bin(opts.transx(x), opts.transy(y), xbins);
0192 binned_y = bin(opts.transy(y), opts.transx(x), ybins);
0193 
0194 if isfield(opts.medprops,'atleastN')
0195     binned_x(cellfun(@numel,binned_x)<opts.medprops(1).atleastN)={NaN};
0196     binned_y(cellfun(@numel,binned_y)<opts.medprops(1).atleastN)={NaN};
0197     opts.medprops = rmfield(opts.medprops,'atleastN');
0198 end
0199 medians_x = cellfun(@median, binned_x);
0200 medians_y = cellfun(@median, binned_y);
0201 
0202 % do some extra things related to median lines before passing options to
0203 % the plot function
0204 if isfield(opts.medprops,'legend')
0205     leg = opts.medprops.legend;
0206     opts.medprops = rmfield(opts.medprops,'legend');
0207 end
0208 
0209 
0210 
0211 t(1) = plot(opts.invtrans(xbins), opts.invtrans(medians_x), opts.medprops);
0212 t(2) = plot(opts.invtrans(medians_y), opts.invtrans(ybins), opts.medprops(length(opts.medprops)));
0213 
0214 if exist('leg','var')
0215     lh = legend(t,leg.names);
0216     leg = rmfield(leg,'names');
0217     for F = fieldnames(leg)'
0218         set(lh,F{1},leg.(F{1}))
0219     end
0220 end
0221 
0222 end
0223 
0224 function largecolor = getColourTable(Nunique,opt)
0225 %% getColorTable Make custom colourtable with white for 0
0226 %
0227 % Custom colour-table
0228 % Nunique should be the number of unique values. This is so that every
0229 % thing except 0 gets a colour
0230 
0231 % greentoblue
0232 % blue to red
0233 %COLORS
0234 tot = ceil(Nunique);
0235 ncolour = min([opt.ncolour,tot]);
0236 white = [1 1 1];
0237 w2b=[(.93:-.93/(ceil(ncolour/3) -1):0)',...
0238     (.93:-.93/(ceil(ncolour/3) -1):0)',...
0239     ones(ceil(ncolour/3),1)]; % white2blue
0240 
0241 b2r=[(0:1/(2*ceil(ncolour/3)-1):1)',...
0242     zeros(2*ceil(ncolour/3),1),...
0243     (1:-1/(2*ceil(ncolour/3)-1):0)']; %blue2red
0244 
0245 colours=[w2b;b2r(2:end,:,:)];
0246 largecolor = zeros(tot,3);
0247 n=0;
0248 step = ceil(tot/ncolour);
0249 for t = 1:size(colours,1)%ncolor
0250     largecolor(n+1:n+step,:) = ...
0251         repmat(colours(t,:),length(n+1:n+step),1);
0252     n = step+n;
0253 end
0254 
0255 largecolor(1,:) = white;
0256 
0257 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005