1 篇文章 0 订阅

# 代码

## 1.极坐标转xy坐标（zjfunc_pol2cart.m）

function v_p = zjfunc_pol2cart(v)

[nx,ny] = size(v);

rn = round((min(nx,ny)-1)/2*sqrt(2));

tn = 360;

theta = 0:.1:359;

x0 = (nx-1)/2+1;

y0 = (ny-1)/2+1;

for i = 1:rn

x(:,i) = x0 + (i)*cosd(theta);

y(:,i) = y0 + (i)*sind(theta);

end

p = 1:nx;%flip(1:rn);

q = 1:ny;%1:rn;

v_p = interp2(p,q,v,x,y);

end

## 2.xy坐标转极坐标（zjfunc_cart2pol.m）

  function [v_xy,R] = zjfunc_cart2pol(v_p)

nx = round(size(v_p,2)*sqrt(2)+1);ny=nx;

rn = (min(nx,ny)-1)/2;

R=1:round(rn*sqrt(2));

x=1:nx;x=x-rn;

y=1:ny;y=y-rn;

theta = 0:.1:359;

for i=1:length(x)

for j=1:length(y)

if y(j)>=0

if x(i)>=0

theta_now(i,j)=round(atan(y(j)/x(i))*360/2/pi);

else

theta_now(i,j)=round(atan(y(j)/x(i))*360/2/pi)+180;

end

else

if x(i)>=0

theta_now(i,j)=round(atan(y(j)/x(i))*360/2/pi);

else

theta_now(i,j)=round(atan(y(j)/x(i))*360/2/pi)-180;

end

end

rho_now(i,j)=round(sqrt(x(i)^2+y(j)^2));

if x(i)==0

if y(j)==0

theta_now(i,j)=0;rho_now(i,j)=0;

elseif y(j)>0

theta_now(i,j)=90;

else

theta_now(i,j)=-90;

end

else

if y(j)==0

if x(i)>0

theta_now(i,j)=0;

else

theta_now(i,j)=-180;

end

end

end

end

end

theta_now=theta_now+180;theta_now(theta_now==360)=0;

v_xy=interp2(theta,R,v_p',theta_now,rho_now);

v_xy(isnan(v_xy))=nanmean(nanmean(v_xy));

v_xy=rot90(v_xy,2);

end

## 3.演示（demozj.m）

    load earth.mat

v=X(1:250,:);

v_p=zjfunc_cart2pol(v);

[v_xy,R]=zjfunc_pol2cart(v_p);

figure

theta = 0:.1:359;

polarPcolor(R,theta,v_p)

figure

contourf(v_xy)

## 4.MATLAB自带绘图程序（polarPcolor.m）

 function [varargout] = polarPcolor(R,theta,Z,varargin)

% [h,c] = polarPcolor1(R,theta,Z,varargin) is a pseudocolor plot of matrix

% Z for a vector radius R and a vector angle theta.

% The elements of Z specify the color in each cell of the

% plot. The goal is to apply pcolor function with a polar grid, which

% provides a better visualization than a cartesian grid.

%

%% Syntax

%

% [h,c] = polarPcolor(R,theta,Z)

% [h,c] = polarPcolor(R,theta,Z,'Ncircles',10)

% [h,c] = polarPcolor(R,theta,Z,'Nspokes',5)

% [h,c] = polarPcolor(R,theta,Z,'Nspokes',5,'colBar',0)

% [h,c] = polarPcolor(R,theta,Z,'Nspokes',5,'labelR','r (km)')

%

% INPUT

% * R :

%        - type: float

%        - size: [1 x Nrr ] where Nrr = numel(R).

% * theta :

%        - type: float

%        - size: [1 x Ntheta ] where Ntheta = numel(theta).

%        - dimension: azimuth or elevation angle (deg).

%        - N.B.: The zero is defined with respect to the North.

% * Z :

%        - type: float

%        - size: [Ntheta x Nrr]

%        - dimension: user's defined .

% * varargin:

%        - Ncircles: number  of circles for the grid definition.

%        - autoOrigin: 'on' (the first circle of the plar grid has a radius

%          equal to the lowest value of R) or 'off'.

%        - Nspokes: number of spokes for the grid definition.

%        - colBar: display the colorbar or not.

%        - labelR: title for radial axis.

%        - RtickLabel: Tick label for the radial axis.

%        - colormap: Colormap for the pcolor function

%        - ncolor: Number of colors in the colorbar and pcolor

%        - circlesPos: position of the circles with respect to the origin

%        (it overwrites Ncircles if necessary)

%

%

% OUTPUT

% h: returns a handle to a SURFACE object.

% c: returns a handle to a COLORBAR object.

%

%% Examples

% R = linspace(3,10,100);

% theta = linspace(0,180,360);

% Z = linspace(0,10,360)'*linspace(0,10,100);

% figure

% polarPcolor(R,theta,Z,'Ncircles',3)

%

%% Author

% Etienne Cheynet, University of Stavanger, Norway. 23/10/2019

%

%%  InputParseer

p = inputParser();

p.CaseSensitive = false;

p.parse(varargin{:});

Ncircles = p.Results.Ncircles;

Nspokes = p.Results.Nspokes ;

labelR = p.Results.labelR ;

RtickLabel = p.Results.RtickLabel ;

colBar = p.Results.colBar ;

Rscale = p.Results.Rscale ;

autoOrigin = p.Results.autoOrigin ;

myColorMap = p.Results.colormap ;

ncolor = p.Results.ncolor ;

circPos = p.Results.circlesPos ;

typeRose = p.Results.typeRose ;

if ~isempty(circPos)

Origin = max([min(circPos),min(R)]);

circPos(circPos<min(R))=[];

circPos(circPos>max(R))=[];

elseif strcmpi(autoOrigin,'on')

Origin = min(R);

elseif strcmpi(autoOrigin,'off')

Origin = 0;

else

error(' ''autoOrigin'' must be ''on'' or ''of'' ')

end

if Origin==0 && strcmpi(Rscale,'log')

warning(' The origin cannot be set to 0 if R is expressed on a logarithmic axis. The value ''Rmin'' is used instead')

Origin = min(R);

end

if isempty(circPos)

if ~isempty(RtickLabel)

if numel(RtickLabel)~=Ncircles

error(' The radial ticklabel must be equal to Ncircles');

end

if any(cellfun(@ischar,RtickLabel)==0)

error(' The radial ticklabel must be a cell array of characters');

end

end

end

if ~isempty(circPos)

circPos = unique([min(R),circPos,max(R)]);

end

%% Preliminary checks

% case where dimension is reversed

Nrr = numel(R);

Noo = numel(theta);

if isequal(size(Z),[Noo,Nrr]) && Noo~=Nrr,

Z=Z';

end

% case where dimension of Z is not compatible with theta and R

if ~isequal(size(Z),[Nrr,Noo])

fprintf('\n')

fprintf([ 'Size of Z is : [',num2str(size(Z)),'] \n']);

fprintf([ 'Size of R is : [',num2str(size(R)),'] \n']);

fprintf([ 'Size of theta is : [',num2str(size(theta)),'] \n\n']);

error(' dimension of Z does not agree with dimension of R and Theta')

end

%% data plot

rMin = min(R);

rMax = max(R);

thetaMin=min(theta);

thetaMax =max(theta);

if strcmpi(typeRose,'meteo')

theta = theta;

elseif strcmpi(typeRose,'default')

theta = 90-theta;

else

error('"type" must be "meteo" or "default" ');

end

% Definition of the mesh

cax = newplot;

Rrange = rMax - rMin; % get the range for the radius

[rNorm] = getRnorm(Rscale,Origin,R,Rrange); % getRnorm is a nested function

YY = (rNorm)'*cosd(theta);

XX = (rNorm)'*sind(theta);

h = pcolor(XX,YY,Z,'parent',cax);

if ~isempty(ncolor)

cmap = feval(myColorMap,ncolor);

colormap(gca,cmap);

else

colormap(gca,myColorMap);

end

% disp([max(R/Rrange),max(rNorm)])

set(cax,'dataaspectratio',[1 1 1]);axis off;

if ~ishold(cax);

hold(cax,'on')

% Draw circles and spokes

createSpokes(thetaMin,thetaMax,Ncircles,circPos,Nspokes);

createCircles(rMin,rMax,thetaMin,thetaMax,Ncircles,circPos,Nspokes)

end

%% PLot colorbar if specified

if colBar==1,

c =colorbar('location','WestOutside');

caxis([quantile(Z(:),0.01),quantile(Z(:),0.99)])

else

c = [];

end

%% Outputs

nargoutchk(0,2)

if nargout==1,

varargout{1}=h;

elseif nargout==2,

varargout{1}=h;

varargout{2}=c;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nested functions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function createSpokes(thetaMin,thetaMax,Ncircles,circlesPos,Nspokes)

spokeMesh = round(linspace(thetaMin,thetaMax,Nspokes));

if isempty(circlesPos)

circleMesh = linspace(rMin,rMax,Ncircles);

else

circleMesh  =  circlesPos;

end

contourD = abs((circleMesh - circleMesh(1))/Rrange+R(1)/Rrange);

if strcmpi(typeRose,'meteo')

cost = cosd(90-spokeMesh); % the zero angle is aligned with North

sint = sind(90-spokeMesh); % the zero angle is aligned with North

elseif strcmpi(typeRose,'default')

cost = cosd(spokeMesh); % the zero angle is aligned with east

sint = sind(spokeMesh); % the zero angle is aligned with east

else

error('"type" must be "meteo" or "default" ');

end

for kk = 1:Nspokes

X = cost(kk)*contourD;

Y = sint(kk)*contourD;

if  Origin==0

X(1)=Origin;

Y(1)=Origin;

end

plot(X,Y,'color',[0.5,0.5,0.5],'linewidth',0.75,...

'handlevisibility','off');

% avoid superimposition of 0 and 360

if and(thetaMin==0,thetaMax == 360),

if spokeMesh(kk)<360,

text(1.05.*contourD(end).*cost(kk),...

1.05.*contourD(end).*sint(kk),...

[num2str(spokeMesh(kk),3),char(176)],...

'horiz', 'center', 'vert', 'middle');

end

else

text(1.05.*contourD(end).*cost(kk),...

1.05.*contourD(end).*sint(kk),...

[num2str(spokeMesh(kk),3),char(176)],...

'horiz', 'center', 'vert', 'middle');

end

end

end

function createCircles(rMin,rMax,thetaMin,thetaMax,Ncircles,circlePos,Nspokes)

if isempty(circlePos)

if Origin ==0 % if the origin is set at rMin

contourD = linspace(0,1+R(1)/Rrange,Ncircles);

else % if the origin is automatically centered at 0

contourD = linspace(0,1,Ncircles)+R(1)/Rrange;

end

else

contourD = circlePos-circlePos(1);

contourD = contourD./max(contourD)*max(R/Rrange);

contourD =[contourD(1:end-1)./contourD(end),1]+R(1)/Rrange;

end

if isempty(circlePos)

if strcmpi(Rscale,'linear')||strcmpi(Rscale,'lin'),

tickMesh = linspace(rMin,rMax,Ncircles);

elseif strcmpi(Rscale,'log')||strcmpi(Rscale,'logarithmic'),

tickMesh  = logspace(log10(rMin),log10(rMax),Ncircles);

else

error('''Rscale'' must be ''log'' or ''linear'' ');

end

else

tickMesh  = circlePos;

Ncircles = numel(tickMesh);

end

% define the grid in polar coordinates

if strcmpi(typeRose,'meteo')

angleGrid = linspace(90-thetaMin,90-thetaMax,100);

elseif strcmpi(typeRose,'default')

angleGrid = linspace(thetaMin,thetaMax,100);

else

error('"type" must be "meteo" or "default" ');

end

xGrid = cosd(angleGrid);

yGrid = sind(angleGrid);

spokeMesh = linspace(thetaMin,thetaMax,Nspokes);

% plot circles

for kk=1:length(contourD)

X = xGrid*contourD(kk);

Y = yGrid*contourD(kk);

plot(X,Y,'color',[0.5,0.5,0.5],'linewidth',1);

end

position = 0.51.*(spokeMesh(min(Nspokes,round(Ncircles/2)))+...

spokeMesh(min(Nspokes,1+round(Ncircles/2))));

if strcmpi(typeRose,'meteo'),position = 90-position; end

if strcmpi(typeRose,'default') && min(90-theta)<5,position = 0; end

if min(round(theta))==90 && strcmpi(typeRose,'meteo'),  position = 0; end

if max(round(theta))==90 && strcmpi(typeRose,'meteo'),  position = 0; end

for kk=1:Ncircles

if isempty(RtickLabel),

rtick = num2str(tickMesh(kk),2);

else

rtick = RtickLabel(kk);

end

t = text(contourD(kk).*cosd(position),...

(contourD(kk)).*sind(position),...

rtick,'verticalalignment','BaseLine',...

'horizontalAlignment', 'right',...

'handlevisibility','off','parent',cax);

if min(round(abs(90-theta)))<5 && strcmpi(typeRose,'default'),

t.Position =  t.Position - [0,0.1,0];

t.Interpreter = 'latex';

clear t;

end

if min(round(theta))==90 && strcmpi(typeRose,'meteo')

t.Position =  t.Position + [0,0.02,0];

t.Interpreter = 'latex';

clear t;

elseif max(round(theta))==90 && strcmpi(typeRose,'meteo')

t.Position =  t.Position - [0,0.05,0];

t.Interpreter = 'latex';

clear t;

end

% annotate spokes

if max(theta)-min(theta)>180,

t = text(contourD(end).*1.3.*cosd(position),...

contourD(end).*1.3.*sind(position),...

[labelR],'verticalalignment','bottom',...

'horizontalAlignment', 'right',...

'handlevisibility','off','parent',cax);

else

t = text(contourD(end).*0.6.*cosd(position),...

contourD(end).*0.6.*sind(position),...

[labelR],'verticalalignment','bottom',...

'horizontalAlignment', 'right',...

'handlevisibility','off','parent',cax);

end

t.Interpreter = 'latex';

if min(round(theta))==90 && strcmpi(typeRose,'meteo'),

t.Position =  t.Position + [0,0.05,0];

clear t;

elseif max(round(theta))==90 && strcmpi(typeRose,'meteo'),

t.Position =  t.Position + [0,0.05,0];

clear t;

end

%                 if min(round(abs(90-theta)))<5 && strcmpi(typeRose,'default'),

%                     t.Position =  t.Position - [0,0.12,0];

%                     t.Interpreter = 'latex';

%                     clear t;

%                 end

end

end

function [rNorm] = getRnorm(Rscale,Origin,R,Rrange)

if strcmpi(Rscale,'linear')||strcmpi(Rscale,'lin')

rNorm = R-R(1)+Origin;

rNorm = (rNorm)/max(rNorm)*max(R/Rrange);

elseif strcmpi(Rscale,'log')||strcmpi(Rscale,'logarithmic')

if rMin<=0

error(' The radial vector cannot be lower or equal to 0 if the logarithmic scale is used');

end

rNorm = log10(R); %normalized radius [0,1]

rNorm =rNorm-rNorm(1);

rNorm = (rNorm)/max(rNorm)*max(R/Rrange);

else

error('''Rscale'' must be ''log'' or ''linear'' ');

end

end

end

• 2
点赞
• 6
收藏
觉得还不错? 一键收藏
• 打赏
• 0
评论
05-20
05-06 3万+
06-13 1万+
09-02
04-19 1万+

### “相关推荐”对你有帮助么？

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

¥2 ¥4 ¥6 ¥10 ¥20

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、C币套餐、付费专栏及课程。