Code covered by the BSD License

# ncquiverref

18 Nov 2007 (Updated )

Quiver (vector) with a reference vector or color code on a map or Cartesian axes.

ncquiverref(x,y,u,v,units,reftype,refvec,veccol,cont)
function ncquiverref(x,y,u,v,units,reftype,refvec,veccol,cont)

% NCQUIVERREF: Vector plotting function for either map or Cartesian axes.
%
% function ncquiverref(x,y,u,v,units,reftype,refvec,veccol,cont)
%
% This function is a substitute for the standard version of quiver
% and quiverm available using a vanilla release of matlab. This version
% assumes a 2D vector field being plotted using a gridded flow field
% from numerical models with a regular geometry. The function enables
% the scaling of vectors according to a reference vector plotted in the
% lower right hand corner of the plot axes. The function works for both
% map and Cartesian axes and allows the color of vectors to be changed.
%
% If a reference value is not provided, the reference value is calculated
% by rounding the median or maximum lengths of the quiver vectors to the
% first significant digit. Scaling of vectors still occurs even
% if the reference vector plotting is switched off. This enables
% different subplots to share identical scaling so that the relative
% magnitude of vectors can be compared between subplots (provided they
% share the same grid).
%
% The function includes the ability to plot color vectors, all of equal
% length but color coded according to their magnitude. In this case, a
% colorbar is provided, complete with units, and no scaling vector is
% plotted.
%
% Input:
% x       - x-coordinate or latitude
%
% y       - y-coordinate or longitude
%
% u       - u-component (Cartesian +x-direction, map +longitude-direction)
%
% v       - v-component (Cartesian +y-direction, map +latitude-direction)
%
% units   - a string providing the units of the vector field.  This assumes
%           the default Tex interpreter is being used as set
%           elsewhere using set(0,'DefaultTextInterpreter','tex').
%
% reftype - character variable specifying the type of reference vector.
%           Allowable values are 'median' for giving a reference
%           vector based on the median, or 'max' for giving the
%           reference vector based on the maximum.  This argument may
%           be omitted with 'max' as the default. If reftype is entered
%           as a number, then this is the value of the reference vector
%           in the units of the data in u and v. If veccol is set to
%           'col', this argument has no effect and should be entered
%           only as a dummy argument.
%
% refvec  - logical that turns off or on the plotting of the reference
%           vector scale for plotting vectors.  The default is 'true'.
%
% veccol  - color of the vectors to be plotted.  This may either be in the
%           form of RGB or as a single letter, such as 'b' for blue
%           using standard Matlab color specifications. This may be set
%           to 'col' if the vectors are to be color coded by magnitude
%           instead of sized by magnitude. In this case, all vectors
%           have the same size based on the optimum value for the grid
%           provided, and a colorbar is provided to reference the values.
%
% cont    - contour levels to be used if color shading the vectors.  For
%           this to work, veccol='col', and as a result all vectors
%           are made of equal length but color coded according to their
%           magnitude. The contour levels must be non-zero, and must include
%           at least one value.   The values at the divider between
%           vector color categories.
%
% Output:
% Output is graphical to the current active figure axes. All vectors
% are centered on the grid points, rather than the starting point of the
% vector being positioned on the grid point. This can easily be changed
% in the code if required, but is the preference of the author.
%
% Written by Andrew Roberts
% Tested using MATLAB Version 7.11.0.584 (R2010b)
%
% $Id: ncquiverref.m 254 2011-02-09 05:26:08Z aroberts$

% make sure the mapping toolbox is present
h=ver('map') ; if length(h)==0 ; error('Mapping toolbox not installed') ; end

% error checking of inputs
if nargin<4; error('Missing vector field input data'); end

% set default units
if nargin<5 ; units=''; end

% set default reftype value
if nargin<6 ; reftype='max' ; end

% default it to turn on the reference vector
if nargin<7 ; refvec=true ; end

% default vector color is blue
col=false;
if nargin<8 ;
veccol='b' ;
elseif ischar(veccol) & strcmp(veccol,'col')
col=true ;
refvec=false;
if nargin<9; error('No contour values provided'); end
end

% get current axis
h=get(gcf,'CurrentAxes');

% use meshgrid if needed
sx=size(x); sy=size(y);
if (sx(1)==1 & sy(1)==1) | (sx(2)==1 & sy(2)==1);
[x,y]=meshgrid(x,y);
elseif sx(1)==1 | sx(2)==1;
error('Dimensions of x and y are inconsistent')
elseif sy(1)==1 | sy(2)==1;
error('Dimensions of x and y are inconsistent')
elseif sx~=sy
error('Dimensions of x and y are inconsistent')
end

% check sizes all agree in input data
if size(x)~=size(y); error('x and y sizes disagree'); end
if size(x)~=size(u); error('x and u sizes disagree'); end
if size(y)~=size(v); error('y and v sizes disagree'); end

% If plotting on a matlab map, determine if the axes are map or cartesian
% coordinates, and if the former calculate mapping to plot axis, and
% then do vector field otherwise just plot the vector field.
if ismap(h)

disp('Quivering on current map axes')

% set lat and lon
lat=x;
lon=y;

% get x and y location on the map
sz=size(x);
mstruct=gcm;
[x,y] = mfwdtran(mstruct,lat,lon,h,'none');
xz=size(x);
if sz~=xz
error('Change in size of x using mfwdtran. Try changing surface to none in the code')
end

% get angle on the map, but do not distort the length according to the projection
% so that all vectors can use the same reference vector.  DO NOT project
% the length of the vector to be different in x and y directions.
[th,z] = cart2pol(u,v);
[thproj,len] = vfwdtran(mstruct,lat,lon,90*ones(size(lat)));

elseif isempty(h)

disp('Creating new Cartesian axes')

% get the magnitude of the vector field
[th,z] = cart2pol(u,v);

% set up axes
axis([min(x(:)) max(x(:)) min(y(:)) max(y(:))])
axis xy;
axis tight;
set(gca,'Layer','bottom');
box on;

else

disp('Quivering on current Cartesian axes')

% get the magnitude of the vector field
[th,z] = cart2pol(u,v);

end

% remove masked grid points from the input by filling coordinates with NaN;
x(isnan(u))=NaN;
y(isnan(u))=NaN;

% Scale the vectors according to the reference arrow vector length based on
% the mean distance between grid points. This is a good measure, as it remains
% constant for multiple plots using the same grid with different values.
x1=abs(diff(x')); x2=abs(diff(x));
y1=abs(diff(y')); y2=abs(diff(y));
[th,z1] = cart2pol(x1,y1); [th,z2] = cart2pol(x2,y2);
scalelength=min(mean(z1(~isnan(z1))),mean(z2(~isnan(z2))));

% Calculate reference vector length based on rounded median
% or maximum value of plot.  The default is median based.
if isnumeric(reftype) & ~col
disp('Calculating reference vector based on input number');
refval=reftype;
elseif strcmp(reftype,'median') & ~col
disp('Calculating reference vector based on median');
z(z==0)=NaN;
refval=median(z(~isnan(z)))
elseif strcmp(reftype,'max') & ~col
disp('Calculating reference vector based on maximum');
refval=max(z(~isnan(z)));
elseif ~col
error('reftype must be either "max" or "median"');
else
disp('Color vectors being used with constant length');
end

% Remove NaN values that will not be plotted
% and turn points into a row of coordinates
u=u(~isnan(x))';
v=v(~isnan(x))';
y=y(~isnan(x))';
x=x(~isnan(x))';

% Set arrow size (1= full length of vector)
arrow=0.40;

% set the length of the vectors to be constant and generate
% a color for each particular vector
if col

% get the magnitude and direction of each vector
[th,z] = cart2pol(u,v);

% check cont has at least one value
if isempty(cont) | length(cont)==0 | any(cont==0)
error('cont must be non-zero and contain at least one value')
end

% arrange contour values from min to max, and
% add an extra value for the purpose of processing
cont=sort(abs(cont));
cont(end+1)=cont(end)+1;

% set the colormap
hc=colormap(jet(length(cont)));

for i=1:length(cont)

if i==1
elseif i==length(cont)
else
end

% Center vectors over grid points

% Get x coordinates of each vector plotted
lx = [xstart; ...
xstart+(1-arrow/3)*(xend-xstart); ...
xend-arrow*(u+arrow*v); ...
xend; ...
xend-arrow*(u-arrow*v); ...
xstart+(1-arrow/3)*(xend-xstart); ...
repmat(NaN,size(xstart))];

% Get y coordinates of each vector plotted
ly = [ystart; ...
ystart+(1-arrow/3)*(yend-ystart); ...
yend-arrow*(v-arrow*u); ...
yend; ...
yend-arrow*(v+arrow*u); ...
ystart+(1-arrow/3)*(yend-ystart); ...
repmat(NaN,size(ystart))];

% Plot the vectors
line(lx,ly,'Color',hc(i,:));

end

% generate the colorbar and get its dimensions
hcb=colorbar('peer',gca,'location','EastOutside');
xlimcb=get(hcb,'xlim');
ylimcb=get(hcb,'ylim');

% initialize the cell array for the colorbar labels
tickn=cell(1,length(cont)+1);

for i=1:length(cont)

% set lower and upper bounds on the y-axis of the colorbar
ycb1=ylimcb(1)+(i-1)*(ylimcb(2)-ylimcb(1))/length(cont);
ycb2=ylimcb(1)+(i)*(ylimcb(2)-ylimcb(1))/length(cont);

% patch the color bar (this is required because of an error in the
% matlab postscript generator which becomes evident when trying to print)
patch([xlimcb(1);xlimcb(1);xlimcb(2);xlimcb(2)],[ycb1;ycb2;ycb2;ycb1],...
hc(i,:),'Parent',hcb);

% set tick labels of colorbar
tickn{i+1}=num2str(cont(i));

end

if max(abs(cont))>=10000 | max(abs(cont))<=1
factor=fix(log10(max(abs(cont))))-1;
for i=1:length(cont)-1
tickn{i+1}=num2str(cont(i)/(10^factor));
end
text(xlimcb(1),ylimcb(1),{' ',['{\times} 10^{',num2str(factor),'} ',units]},...
'Parent',hcb,...
'HorizontalAlign','left','VerticalAlign','top')
else
text(xlimcb(1),ylimcb(1),{' ',units},...
'Parent',hcb,...
'HorizontalAlign','left','VerticalAlign','top')
end

% Set the upper tick to units and turn on the new labels
tickn{end}='';
tickn{1}='0';
tickpos=ylimcb(1)+diff(ylimcb)*[0:length(cont)]/length(cont);
set(hcb,'YTickLabel',tickn,'YTickMode','manual','YTick',tickpos)

% Add up arrow for top of colorbar
%text(xlimcb(2),ylimcb(2),[' $\uparrow$'],'Interpreter','LaTex','Parent',hcb,...
%     'HorizontalAlign','left','VerticalAlign','top')

yl=diff(ylimcb)*0.03;
yt=ylimcb(2);
ys=diff(ylimcb)*0.012;
ya=yt-ys;
yb=ylimcb(2)-yl;
xc=xlimcb(2)+0.25*diff(xlimcb);
xl=xc-ys; xr=xc+ys;
line([xc xc xl xc xr],[yb yt ya yt ya],'color','k',...
'Clipping','off','Parent',hcb,'linewidth',0.3)

else

% set scale value based on refval and scale length
roundp=floor(log10(refval));
refval=floor(refval/(10^roundp))*(10^roundp);
scale=scalelength/refval;

% Center vectors over grid points
xstart=x-0.5*scale*u;
xend=x+0.5*scale*u;
ystart=y-0.5*scale*v;
yend=y+0.5*scale*v;

% Get x coordinates of each vector plotted
lx = [xstart; x; ...
xstart+(1-arrow/3)*(xend-xstart); ...
xend-arrow*(scale*u+arrow*(scale*v)); ...
xend; ...
xend-arrow*(scale*u-arrow*(scale*v)); ...
xstart+(1-arrow/3)*(xend-xstart); ...
repmat(NaN,size(x))];

% Get y coordinates of each vector plotted
ly = [ystart; y; ...
ystart+(1-arrow/3)*(yend-ystart); ...
yend-arrow*(scale*v-arrow*(scale*u)); ...
yend; ...
yend-arrow*(scale*v+arrow*(scale*u)); ...
ystart+(1-arrow/3)*(yend-ystart); ...
repmat(NaN,size(y))];

% Plot the vectors
line(lx,ly,'Color',veccol);

end

% Draw the reference vector key at altitude 2 above the map and grid
if refvec

% Get the reference text string, formatted to powers of ten if required
if refval < 0.1 | refval > 100
factor=floor(log10(refval));
reftext=[num2str(refval/(10^factor)),' \times 10^{',num2str(factor),'} ',units,' '];
else
reftext=[num2str(refval),' ',units,' '];
end

% Get the current axis limits
xlim=get(gca,'xlim'); xp1=xlim(1); xp2=xlim(2);
ylim=get(gca,'ylim'); yp1=ylim(1); yp2=ylim(2);

% set padding around the reference vector

% Set x position of reference vector
xstart=xend-scalelength;

% Plot reference text in lower right hand corner
'VerticalAlignment','Bottom','HorizontalAlignment','Right');
textextent=get(ht,'Extent');

% Draw patch over area of vector key
xr=xp2;
yb=yp1;
hp=patch([xl; xl; xr; xr],[yb; yt; yt; yb],[2; 2; 2; 2],'w',...
'LineWidth',0.5,'Parent',gca);
uistack(hp,'top');

% Redraw reference text on top of patch
ht=text(xstart,(yb+yt)/2,2.1,reftext,'Parent',gca,'FontSize',8.5,...
'VerticalAlignment','Middle','HorizontalAlignment','Right');

% Set y position of reference vector
yend=textextent(2)+textextent(4)/2;
ystart=yend;

% Get x coordinates of reference vector plotted
lx = [xstart; ...
xstart+(1-arrow/3)*(xend-xstart); ...
xend-arrow*scalelength; ...
xend; ...
xend-arrow*scalelength; ...
xstart+(1-arrow/3)*(xend-xstart); ...
NaN];

% Get y coordinates of reference vector plotted
ly = [ystart; ...
ystart+(1-arrow/3)*(yend-ystart); ...
yend+arrow*(arrow*scalelength); ...
yend; ...
yend-arrow*(arrow*scalelength); ...
ystart+(1-arrow/3)*(yend-ystart); ...
NaN];

% Get z coordinates of reference vector
lz = 2*ones(size(ly));

% Plot the reference vector
line(lx,ly,lz,'Color',veccol);

else

disp('No reference vector is being plotted')

end

end % function