| 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
% Naval Postgraduate School
% 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)));
[u,v] = pol2cart(th+deg2rad(thproj),z);
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
mask=find(z<cont(i));
elseif i==length(cont)
mask=find(z>=cont(i-1));
else
mask=find(z<cont(i) & z>=cont(i-1));
end
% Center vectors over grid points
[u,v] = pol2cart(th(mask),scalelength);
xstart=x(mask)-0.5*u;
xend=x(mask)+0.5*u;
ystart=y(mask)-0.5*v;
yend=y(mask)+0.5*v;
% 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
% Add exponent to colorbar if required and add units
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
padx=diff(xlim)/100;
pady=diff(ylim)/100;
% Set x position of reference vector
xend=xp2-padx;
xstart=xend-scalelength;
% Plot reference text in lower right hand corner
ht=text(xstart,yp1+pady,reftext,'Visible','off','Parent',gca,'FontSize',8.5,...
'VerticalAlignment','Bottom','HorizontalAlignment','Right');
textextent=get(ht,'Extent');
% Draw patch over area of vector key
xl=textextent(1)-padx;
xr=xp2;
yb=yp1;
yt=textextent(2)+textextent(4)+pady;
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
|
|
