不改变数据形状的极坐标转换这应是原创,网上能找到的都是会改变数据形状的,比如正方形的xy图变为圆形的极坐标图。
而我的程序则是直接投影,即正方形的xy图变为正方形的数据区块投影在极坐标图上。
一共四个文件:zjfunc_pol2cart.m,zjfunc_cart2pol.m,demozj.m,polarPcolor.m。使用时放于同一个文件夹下,运行demozj.m即可得到下图。
极坐标(左)xy坐标(右)
原理
1.极坐标转xy坐标
利用interp2函数,即interp2(theta,R,v_p’,theta_now,rho_now)。首先对待转换的极坐标下的矩阵进行插值重构为xy坐标矩阵,将所需xy坐标对应的theta和rho计算出来存为极坐标下的theta_now和rho_now,通过查询theta_now和rho_now在插值重构后的xy坐标矩阵中对应位置的值获得xy坐标的结果。
2.xy坐标转极坐标
同样的,利用interp2函数,即interp2(p,q,v,x,y)。首先对待转换的xy坐标下的矩阵进行插值重构为极坐标矩阵,将所需极坐标对应的x和y计算出来存为xy坐标下的x和y,通过查询x和y在插值重构后的极坐标矩阵中对应位置的值获得极坐标的结果。其中涉及对不同象限数据的分类讨论,较xy坐标转极坐标略复杂。其中涉及对不同象限数据的分类讨论,较极坐标转xy坐标略复杂。
使用方法
修改demozj.m中的v矩阵为需要转换的矩阵。
局限性
由于时间问题,目前程序只能转换正方形数据,长宽不一致的数据需要适当修改后使用喔。
代码
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).
% - dimension: radial distance.
% * 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
% see also pcolor
%
%% InputParseer
p = inputParser();
p.CaseSensitive = false;
p.addOptional('Ncircles',5);
p.addOptional('autoOrigin','on');
p.addOptional('Nspokes',8);
p.addOptional('labelR','');
p.addOptional('RtickLabel',[]);
p.addOptional('colBar',1);
p.addOptional('Rscale','linear');
p.addOptional('colormap','parula');
p.addOptional('ncolor',[]);
p.addOptional('typeRose','meteo'); % 'meteo' or 'default'
p.addOptional('circlesPos',[]);
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)])
shading flat
set(cax,'dataaspectratio',[1 1 1]);axis off;
if ~ishold(cax);
% make a radial grid
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');
% plot graduations of angles
% 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
% radius tick label
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
% radial graduations
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
今天的文章MATLAB极坐标与xy坐标互相转换_不改变数据形状_极坐标变量v_p(theta,r)的平面图分享到此就结束了,感谢您的阅读,如果确实帮到您,您可以动动手指转发给其他人。
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/33202.html