根据需要，修改了一些参数：1. 字体全部统一为新罗马；2. 颜色改为彩色。
function varargout = wind_rose(D,F,varargin)%WIND_ROSE Wind rose of direction and intensity%% Syntax:%[HANDLES,DATA] = WIND_ROSE(D,I,VARARGIN)%% Inputs:%D Directions%I Intensities%VARARGIN:% -dtype, type of input directions D, standard or meteo,% if meteo, the conversion dnew=mod(-90-D,360) is done;% if not meteo, standard is used (default)% -n, number of D subdivisons% -di, intensities subdivisons, default is automatic% -ci, percentage circles to draw, default is automatic% -labtitle, main title% -lablegend, legend title% -cmap, colormap [jet]% -colors, to use instead of colormap, for each di% -quad, Quadrant to show percentages [1]% -ri, empty internal radius, relative to size of higher% percentage [1/30]% -legtype, legend type: 1, continuous, 2, separated boxes [2]% -bcolor, full rectangle border color ['none']% -lcolor, line colors for axes and circles ['k']% -percbg, percentage labels bg ['w']% -ax, to place wind rose on pervious axes, the input for ax% must be [theax x y width], where theax is the previous% axes, x and y are the location and width is the wind% rose width relative to theax width (default=1/5)% -parent, by default a new axes is created unless parent is%given, ex, parent may be a subplot% -iflip, flip the intensities as they go outward radially, ie,%highest values are placed nearest the origin [{0} 1]% -inorm, normalize intensities, means all angles will have 100%% -incout, if 0, data outside di limits will not be used [0 {1}]%% Output:%HANDLES Handles of all lines, fills, texts%DATA Wind rose occurences per direction and intensity%% Examle:%d=0:10:350;%D=[];%V=[];%for i=1:length(d)% n=d(i)/10;% D=[D ones(1,n)*d(i)];% V=[V 1:n];%end%%figure%wind_rose(D,V)%%figure%wind_rose(D,V,'iflip',1)%%figure%wind_rose(D,V,'ci',[1 2 7],'dtype','meteo')%%% place it on a previous axes:%ax=axes;%plot(lon,lat)%wind_rose(D,V,'ax',[ax x y 1/3])%% MMA 26-11-, mma@%% IEO, Instituto Espaol de Oceanografa% La Corua, Espaa% 10-12- - Added varargin ci and n (nAngles removed as input)% 17-12- - Added varargin ax, colors% 22-02- - Added varargin dtype% 08-05- - Bug fix (bar at dir=0 could be incorrect in some cases)% 14-05- - Added varargin iflip% 16-06- - Added varargin parent% 10-06- - Added varargin incout% 27-04- - Added output DATA% 17-06- - Bug fix (E(i,end)=length(find(b>=Ag(end-1))),%previously was ...b>Ag...). So the percentages where%wrong only when using intensities equal to the lower%value of the highest intensity subdivision, basically%an academic case.handles=[];% varargin options:dtype='standard';nAngles=36;ri=1/30;quad=1;legType=2;percBg='w';titStr='';legStr='';% cmap=jet;cmap=jet;colors=[];Ag=[]; % intensity subdivs.ci=[]; % percentage circleslineColors='k';borderColor='none';onAxes=false;iflip=0;inorm=0;parent=0;IncHiLow=1; % include values higher and lower that the limits of Ag.vin=varargin;for i=1:length(vin)if isequal(vin{i},'dtype')dtype=vin{i+1};elseif isequal(vin{i},'n')nAngles=vin{i+1};elseif isequal(vin{i},'ri')ri=vin{i+1};elseif isequal(vin{i},'quad')quad=vin{i+1};elseif isequal(vin{i},'legtype')legType=vin{i+1};elseif isequal(vin{i},'percbg')percBg=vin{i+1};elseif isequal(vin{i},'labtitle')titStr=vin{i+1};elseif isequal(vin{i},'lablegend')legStr=vin{i+1};elseif isequal(vin{i},'cmap')cmap=vin{i+1};elseif isequal(vin{i},'colors')colors=vin{i+1};elseif isequal(vin{i},'di')Ag=vin{i+1};elseif isequal(vin{i},'ci')ci=vin{i+1};elseif isequal(vin{i},'lcolor')lineColors=vin{i+1};elseif isequal(vin{i},'bcolor')borderColor=vin{i+1};elseif isequal(vin{i},'ax')ax=vin{i+1};tryonAxes=ax(1);onAxesX=ax(2);onAxesY=ax(3);onAxesR=ax(4);catchdisp(':: cannot place wind rose on axes, bad argument for ax')returnendelseif isequal(vin{i},'iflip')iflip=vin{i+1};elseif isequal(vin{i},'inorm')inorm=vin{i+1};elseif isequal(vin{i},'parent')parent=vin{i+1};elseif isequal(vin{i},'incout')IncHiLow=vin{i+1};endend% other options:% size of the full rectangle:rs=1.2;rl=1.7;% directions conversion:if isequal(dtype,'meteo')% D=mod(-90-D,360);D=mod(90-D,360);end% angles subdivisons:D=mod(D,360);Ay=linspace(0,360,nAngles+1)-0.5*360/nAngles;% calc instensity subdivisions:if isempty(Ag)% gen Ag:f=figure('visible','off');plot(F); axis tight;yl=get(gca,'ytick');close(f)dyl=diff(yl); dyl=dyl(1);if min(F)<yl(1), yl=[yl(1)-dyl yl]; endif max(F)>yl(end), yl=[yl yl(end)+dyl]; endAg=yl;endfor i=1:length(Ay)-1if i==1I=find( (D>=Ay(i) & D<Ay(i+1)) | D>=Ay(end));elseI=find(D>=Ay(i) & D<Ay(i+1));endb=F(I);for j=1:length(Ag)-1if j==length(Ag)-1J=find(b>=Ag(j) & b<=Ag(j+1)); % include data with last AggelseJ=find(b>=Ag(j) & b<Ag(j+1));endE(i,j)=length(J);endif IncHiLowE(i,1)=length(find(b<Ag(2)));E(i,end)=length(find(b>=Ag(end-1)));endendb=sum(E,2)/length(D)*100;% normalize data:if inormn=sum(E,2);for i=1:length(n)E(i,:)=E(i,:)/n(i);endb=100*ones(size(b));end% check if has values higher or lower than the Ag limitshasH=length(find(F>Ag(end)));hasL=length(find(F<Ag(1)));% calc number of percentage circles to draw:if isempty(ci)if inormci=[25 50 75];g=120;ncircles=3;elsedcircles=[1 2 5 10 15 20 25 30 50];ncircles=3;d=abs(1./(dcircles/max(b))-ncircles);i=find(d==min(d));d=dcircles(i(1));if d*ncircles<max(b)ncircles=ncircles+1;endci=[1:ncircles]*d;g=ncircles*d;endelsencircles=length(ci);g=max(max(ci),max(b));end% plot axes, percentage circles and percent. data:if parentwrAx=parent;set(wrAx,'units','normalized');elsewrAx=axes('units','normalized');endri=g*ri;handles(end+1)=fill([-rs*g rl*g rl*g -rs*g],[-rs*g -rs*g rs*g rs*g],'w',...'EdgeColor',borderColor);if onAxesset(handles(end),'facecolor','none')endhold onhandles(end+1)=plot([-g-ri -ri nan ri g+ri nan 0 0 nan 0 0],...[0 0 nan 0 0 nan -g-ri -ri nan ri g+ri],':','color',lineColors);t0=[0:360]*pi/180;labs=[];Ang=[1/4 3/4 5/4 7/4]*pi;Valign={'top' 'top' 'bottom' 'bottom'};Halign={'right' 'left' 'left' 'right'};for i=1:ncirclesx=(ci(i)+ri)*cos(t0);y=(ci(i)+ri)*sin(t0);circles(i)=plot(x,y,':','color',lineColors);handles(end+1)=circles(i);labs(i)=text((ci(i)+ri)*cos(Ang(quad)),(ci(i)+ri)*sin(Ang(quad)),[num2str(ci(i)),'%'],...'VerticalAlignment',Valign{quad},'HorizontalAlignment',Halign{quad},...'BackgroundColor',percBg,'FontSize',8,'Fontname', 'Times New Roman');endhandles=[handles labs];% calc colors:if isempty(colors)cor={};for j=1:length(Ag)-1cor{j}=caxcolor(Ag(j),[Ag(1) Ag(end-1)],cmap);endelsecor=colors;end% fill data:n=sum(E,2);if iflip, E=fliplr(E); endfor i=1:length(Ay)-1if n(i)t=linspace(Ay(i),Ay(i+1),20)*pi/180;r1=ri;for j=1:length(Ag)-1r2=E(i,j)/n(i) *b(i) +r1;x=[r1*cos(t(1)) r2*cos(t) r1*cos(fliplr(t))];y=[r1*sin(t(1)) r2*sin(t) r1*sin(fliplr(t))];if iflip, jcor=length(Ag)-1-j+1;else, jcor=j;endif E(i,j)>0, handles(end+1)=fill(x,y,cor{jcor}); endr1=r2;endendendaxis equalaxis off% uistack has problems in some matlab versions, so:%uistack(labs,'top')%uistack(circles,'top')ch=get(wrAx,'children');if inorm% only bring circles up in inorm case.for i=1:length(circles)ch(ch==circles(i))=[]; ch=[circles(i); ch];endendfor i=1:length(labs)ch(ch==labs(i))=[]; ch=[labs(i); ch];endset(wrAx,'children',ch);% N S E W labels:%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bg='none';args={'BackgroundColor',bg,'FontSize',8,'Fontname', 'Times New Roman'};%%%%%%%%%%原来是8，现在是16%%%%%%%%%%%%%h(1)=text(-g-ri, 0,'WEST', 'VerticalAlignment','top','HorizontalAlignment','left', args{:});h(2)=text( g+ri, 0,'EAST', 'VerticalAlignment','top','HorizontalAlignment','right',args{:});h(3)=text( 0,-g-ri,'SOUTH','VerticalAlignment','bottom','HorizontalAlignment','left', args{:});h(4)=text( 0, g+ri,'NORTH','VerticalAlignment','top','HorizontalAlignment','left', args{:});handles=[handles h];% scale legend:L=(g*rl-g-ri)/7;h=(g+ri)/10;dy=h/3;x0=g+ri+(g*rl-g-ri)/7;x1=x0+L;y0=-g-ri;if legType==1 % contimuous.图例for j=1:length(Ag)-1lab=num2str(Ag(j));if j==1 & hasL & IncHiLowlab='';endy1=y0+h;handles(end+1)=fill([x0 x1 x1 x0],[y0 y0 y1 y1],cor{j});handles(end+1)=text(x1+L/4,y0,lab,'VerticalAlignment','middle','fontsize',8,'Fontname', 'Times New Roman');y0=y1;endif ~ (hasH & IncHiLow)handles(end+1)=text(x1+L/4,y0,num2str(Ag(end)),'VerticalAlignment','middle','fontsize',8,'Fontname', 'Times New Roman');endelseif legType==2 % separated boxes.for j=1:length(Ag)-1lab=[num2str(Ag(j)) ' - ' num2str(Ag(j+1))];if j==1 & hasL & IncHiLowlab=['<',num2str(Ag(2))];endif j==length(Ag)-1 & hasH & IncHiLowlab=['>=',num2str(Ag(j))];endy1=y0+h;handles(end+1)=fill([x0 x1 x1 x0],[y0+dy y0+dy y1 y1],cor{j});handles(end+1)=text(x1+L/4,(y0+dy+y1)/2,lab,'VerticalAlignment','middle','fontsize',8,'Fontname', 'Times New Roman');y0=y1;endend% title and legend label: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%x=mean([-g*rs,g*rl]);y=mean([g+ri,g*rs]);handles(end+1)=text(x,y,titStr,'HorizontalAlignment','center');x=x0;y=y1+dy;handles(end+1)=text(x,y,legStr,'HorizontalAlignment','left','VerticalAlignment','bottom','FontName','Times New Roman');if onAxesplace_wr(onAxes,wrAx,onAxesX,onAxesY,onAxesR);endif nargout>=1varargout{1}=handles;endif nargout>=2varargout{2}=E;endfunction place_wr(ax,ax2,x,y,width)if nargin < 5width=1/5;enduax=get(ax,'units');pax=get(ax,'position');set(ax,'units',uax)axXlim=get(ax,'xlim');axYlim=get(ax,'ylim');x_ax2=pax(1)+pax(3)*(x-axXlim(1))/diff(axXlim);y_ax2=pax(2)+pax(4)*(y-axYlim(1))/diff(axYlim);pax2=get(ax2,'position');width=pax(3)*width;height=pax2(4)*width/pax2(3);pax2=[x_ax2 y_ax2 width height];if 1% place at centre of the wr, not the bottom left corner:ax2Xlim=get(ax2,'xlim');ax2Ylim=get(ax2,'ylim');dx=(0-ax2Xlim(1))/diff(ax2Xlim)*pax2(3);dy=(0-ax2Ylim(1))/diff(ax2Ylim)*pax2(4);x_ax2=x_ax2-dx;y_ax2=y_ax2-dy;pax2=[x_ax2 y_ax2 width height];endset(ax2,'position',pax2)function cor = caxcolor(val,cax,cmap)%CAXCOLOR Caxis color for value% Find the color for a given value in a colormap.%% Syntax:%COLOR = CAXCOLOR(VALUE,CAXIS,COLORMAP)%% Inputs:%VALUE%CAXIS Default is current caxis%COLORMAP Default is current colormap%% Output:%COLOR RGB color vector%% Example:%figure%pcolor(peaks)%color=caxcolor(0);%set(gcf,'color',color)%% MMA 28-5-, martinho@fis.ua.pt% Department of Physics% University of Aveiro, Portugalif nargin < 3cmap = get(gcf,'colormap');endif nargin < 2cax = caxis;endn=size(cmap,1);i= (val-cax(1))/diff(cax) * (n-1) +1;a=i-floor(i);i=floor(i);i=min(i,n);i=max(i,1);if i==ncor=cmap(n,:);elseif i==1cor=cmap(1,:);elsecor=cmap(i,:)*(1-a) + cmap(i+1,:)*a;end
