三元相图绘图 python

示例

---

对上期的补充

本期文章名叫三元相映射图绘制，但不单单讲这一点，同时会对上一篇多元映射地图做出补充：

上篇中可以绘制这样的三变量映射地图：

但这只是对多边形填色，有人问我如果是网格数据或者像素级别的数据该如何填色？实际上只要输入变量是三个同大小的矩阵就行，因此像素填色完全可以做到，本篇先对此做出解答，举个实例(这里也是懒得搜集数据，直接将经纬度及海拔作为映射变量)：

代码中出现的multiVarMap3函数请见上篇文章：

% 导入数据并绘制贴图版地图
load korea5c
worldmap(korea5c,korea5cR);
koreaHdl=geoshow(korea5c,korea5cR,'DisplayType','texturemap');

A=koreaHdl.XData;A=A(1:end-1,1:end-1);
B=koreaHdl.YData;B=B(1:end-1,1:end-1);
C=koreaHdl.CData;


% % 调用工具函数生成图例和映射表
[CMapData,CMapHdl]=multiVarMap3(A,B,C,'colorList',1,'pieceNum',50);
koreaHdl.CData=CMapData;

---

三元相映射图绘制

使用基本说明

三元相映射图绘制函数multiVarMapTri基本使用方法和上期推送的multiVarMap3函数一模一样，

依旧是colorList代表配色，pieceNum代表颜色划分数量：

% 导入数据并绘制贴图版地图
load korea5c
worldmap(korea5c,korea5cR);
koreaHdl=geoshow(korea5c,korea5cR,'DisplayType','texturemap');

A=koreaHdl.XData;A=A(1:end-1,1:end-1);
B=koreaHdl.YData;B=B(1:end-1,1:end-1);
C=koreaHdl.CData;


% % 调用工具函数生成图例和映射表
[CMapData,CMapHdl]=multiVarMapTri(A,B,C,'colorList',6,'pieceNum',8);
koreaHdl.CData=CMapData;

图例可拖拽移动：

配色展示

colorList可以是基础配色编号1-6，或者是3x3大小的RGB颜色数组。

准备了六个基础配色，例如改用配色5：

% 导入数据并绘制贴图版地图
load korea5c
worldmap(korea5c,korea5cR);
koreaHdl=geoshow(korea5c,korea5cR,'DisplayType','texturemap');

A=koreaHdl.XData;A=A(1:end-1,1:end-1);
B=koreaHdl.YData;B=B(1:end-1,1:end-1);
C=koreaHdl.CData;


% % 调用工具函数生成图例和映射表
[CMapData,CMapHdl]=multiVarMapTri(A,B,C,'colorList',5,'pieceNum',8);
koreaHdl.CData=CMapData;

全部配色展示：

示例2多边形地图填色

axm=usamap("conus");

states=shaperead("usastatelo.shp",'UseGeoCoords',true);
% 俩州离太远画不开，不要
for i=length(states):-1:1
    if states(i).Name=="Alaska"||states(i).Name=="Hawaii"
        states(i)=[];
    end
end

faceColors=makesymbolspec('Polygon',{'INDEX',[1 numel(states)],'FaceColor',...
  polcmap(numel(states))});
usaMap=geoshow(states, 'DisplayType','polygon','SymbolSpec', faceColors);

% 不想去收集数据了，随便生成点随机数
A=rand(numel(states),1);
B=rand(numel(states),1);
C=rand(numel(states),1);

% 调用工具函数生成图例和映射表
[CMapData,CMapHdl]=multiVarMapTri(A,B,C,'colorList',2,'pieceNum',8);
% 循环修改每个州颜色
for i=1:length(usaMap.Children)
    usaMap.Children(i).FaceColor=[CMapData(i,1,1),CMapData(i,1,2),CMapData(i,1,3)];
end

示例3气泡图

x=1:20;
y=rand(1,20);
sz=rand(1,20);
BC=bubblechart(x,y,sz); 

% 不想去收集数据了，随便生成点随机数
A=rand(1,20);
B=rand(1,20);
C=rand(1,20);

% 调用工具函数生成图例和映射表
[CMapData,CMapHdl]=multiVarMapTri(A,B,C,'colorList',2,'pieceNum',8);
BC.CData=[CMapData(:,:,1)',CMapData(:,:,2)',CMapData(:,:,3)'];

---

三元相映射图完整代码

function [CMapData,CMapHdl]=multiVarMapTri(A,B,C,varargin) 
% @author : slandarer
% gzh  : slandarer随笔

help multiVarMapTri
% =========================================================================
% 基本属性提取
obj.arginList={'colorList','pieceNum','varName','title'};
obj.colorList=1;
obj.pieceNum=6;
obj.varName={'var1','var2','var3'};
obj.title='ColorMap';
for i=1:2:(length(varargin)-1)
    tid=ismember(obj.arginList,varargin{i});
    if any(tid)
        obj.(obj.arginList{tid})=varargin{i+1};
    end
end
% =========================================================================
% 基础配色表
baseAtla{1}=[1,0,0;0,1,0;0,0,1];
baseAtla{2}=[0,211,166;211,152,255;233,164,67]./255;
baseAtla{3}=[249,3,252;252,242,8;2,251,249]./255;
baseAtla{4}=[215,0,102;0,137,52;3,79,162]./255;
baseAtla{5}=[243,180,0;180,102,0;0,0,0]./255;
baseAtla{6}=[31,5,94;128,128,128;180,0,0]./255;
if numel(obj.colorList)==1,atla=baseAtla{obj.colorList};end
% =========================================================================
% 数据处理及图例色卡数据构建
% CMapData=getCMap3(A,B,C,atla,obj.pieceNum);
    function CMap=getCMap3(A,B,C,atla,n)
        nA=(A(:)-min(A(:)))./(max(A(:))-min(A(:)));
        nB=(B(:)-min(B(:)))./(max(B(:))-min(B(:)));
        nC=(C(:)-min(C(:)))./(max(C(:))-min(C(:)));

        nA=(ceil(nA.*n.*3)-1)./(n.*3-1);
        nB=(ceil(nB.*n.*3)-1)./(n.*3-1);
        nC=(ceil(nC.*n.*3)-1)./(n.*3-1);
        nA(nA<0)=0;nB(nB<0)=0;nC(nC<0)=0;
        
        ABC=[nA,nB,nC];ABC=ABC./sum(ABC,2);ABC(isnan(ABC))=1/3;
        colorAB=atla(1,:).*ABC(:,1)+...
                atla(2,:).*ABC(:,2)+...
                atla(3,:).*ABC(:,3);

        CMap=zeros([size(A),3]);
        CMap(:,:,1)=reshape(colorAB(:,1),size(A));
        CMap(:,:,2)=reshape(colorAB(:,2),size(A));
        CMap(:,:,3)=reshape(colorAB(:,3),size(A));
    end

% =========================================================================
% 坐标区域初定位及修饰
ax1=gca;fig=ax1.Parent;
ax1Pos=ax1.Position;
ax2=axes('Parent',fig,'Position',[ax1Pos(1)+6/9*ax1Pos(3),...
    ax1Pos(2)+6/9*ax1Pos(4),ax1Pos(3)*2.5/9,ax1Pos(4)*2.5/9]);
ax1.UserData=1;
ax2.UserData=2;
ax2.XLim=[-.5,.5];
ax2.YLim=[-sqrt(3)/6,sqrt(3)/3.*1.3];
ax2.XTick=[];
ax2.YTick=[];
ax2.XColor='none';
ax2.YColor='none';
ax2.Color='none';
ax2.DataAspectRatio=[1,1,1];
ax2.FontName='Arial';
ax2.FontSize=9;
ax2.Title.String=obj.title;
ax2.Toolbar.Visible='off';
hold(ax2,'on');view(2);
% =========================================================================
% 图例中三元相图绘制
tt=[pi/2,pi/2+2*pi/3,pi/2+4*pi/3,pi/2];
tX=cos(tt)./sqrt(3);
tY=sin(tt)./sqrt(3);
plot(ax2,tX,tY,'Color',[0,0,0],'LineWidth',1.5)
YList=linspace(tY(1),tY(2),obj.pieceNum+1);
X0List=linspace(tX(1),tX(2),obj.pieceNum+1);
X1List=linspace(tX(1),tX(3),obj.pieceNum+1);
VList=zeros(round((obj.pieceNum+1)*(obj.pieceNum+2)/2),2);
for i=1:length(YList)
    tXList=linspace(X0List(i),X1List(i),i)';
    tYList=YList(i).*(ones(i,1));
    tN=round(i*(i-1)/2);
    VList((tN+1):(tN+i),:)=[tXList,tYList];
end
RList=zeros(obj.pieceNum*obj.pieceNum,3);
FList=zeros(obj.pieceNum*obj.pieceNum,3);n=1;
for i=2:length(YList)
    NLList=(round((i-1)*(i-2)/2)+1):(round((i-1)*(i-2)/2)+i-1);
    NList=(round(i*(i-1)/2)+1):(round(i*(i-1)/2)+i);
    for j=1:(i-1)
        FList(n,:)=[NList(j),NList(j+1),NLList(j)];
        RList(n,:)=[j-1,i-1-j,obj.pieceNum+1-i];
        n=n+1;
    end
end
for i=2:(length(YList)-1) 
    NList=(round(i*(i-1)/2)+1):(round(i*(i-1)/2)+i);
    NNList=(round((i+1)*(i)/2)+1):(round((i+1)*(i)/2)+i+1);
    for j=1:(i-1)
        FList(n,:)=[NList(j),NList(j+1),NNList(j+1)];
        RList(n,:)=[j-1+1/3,i-1-j+1/3,obj.pieceNum+1-i-2/3];
        n=n+1;
    end
end
legendCMap=getCMap3(RList(:,1),RList(:,2),RList(:,3),atla,obj.pieceNum);
patch(ax2,'Faces',FList,'Vertices',VList,'CData',legendCMap,...
    'FaceColor','flat','EdgeColor','none')
text(ax2,cos(pi/2+2*pi/3).*0.65,sin(pi/2+2*pi/3).*0.8,obj.varName{1},'FontSize',11,'FontName','Arial','HorizontalAlignment','center')
text(ax2,cos(pi/2).*0.69,sin(pi/2).*0.69,obj.varName{2},'FontSize',11,'FontName','Arial','HorizontalAlignment','center')
text(ax2,cos(pi/2-2*pi/3).*0.68,sin(pi/2-2*pi/3).*0.8,obj.varName{3},'FontSize',11,'FontName','Arial','HorizontalAlignment','center')


text(ax2,cos(-pi/2).*0.4,sin(-pi/2).*0.4,'→','FontSize',15,'FontName','Arial','HorizontalAlignment','center')
text(ax2,cos(-pi/2+2*pi/3).*0.4,sin(-pi/2+2*pi/3).*0.4,'→','FontSize',15,'FontName','Arial','HorizontalAlignment','center','Rotation',120)
text(ax2,cos(-pi/2-2*pi/3).*0.4,sin(-pi/2-2*pi/3).*0.4,'→','FontSize',15,'FontName','Arial','HorizontalAlignment','center','Rotation',-120)

RList=RList./sum(RList,2);
XY1=[RList(:,1)-RList(:,3)./2,RList(:,3).*sqrt(3)./2];
AA=(A(:)-min(A(:)))./(max(A(:))-min(A(:)));
BB=(B(:)-min(B(:)))./(max(B(:))-min(B(:)));
CC=(C(:)-min(C(:)))./(max(C(:))-min(C(:)));
ABCO=[AA,BB,CC];ABCO=ABCO./sum(ABCO,2);ABCO(isnan(ABCO))=1/3;
XY2=[ABCO(:,1)-ABCO(:,3)./2,ABCO(:,3).*sqrt(3)./2];
DIS=pdist2(XY1,XY2);
[~,IND]=min(DIS,[],1);
CMapData(:,:,1)=reshape(legendCMap(IND,:,1),size(A));
CMapData(:,:,2)=reshape(legendCMap(IND,:,2),size(B));
CMapData(:,:,3)=reshape(legendCMap(IND,:,3),size(C));

% 图例拖动回调
set(fig,'WindowButtonDownFcn',@BtnDownFunc, ...
        'WindowButtonUpFcn',@BtnUpFunc, ...
        'WindowButtonMotionFcn',@BtnMotionFunc, ...
        'KeyPressFcn',@keyFunc)
selectedAx=false;
    % 选择图例时的回调
    function BtnDownFunc(~,~)
        taxes=get(fig,'currentAxes');
        if taxes.UserData==2
            selectedAx=true;
        end
    end
    % 松开图例时的回调
    function BtnUpFunc(~,~)
        selectedAx=false;
        set(fig,'currentAxes',ax1)
    end
    % 选择图例并拖动时的回调
    function BtnMotionFunc(~,~)
        if selectedAx
            tpos=get(fig,'CurrentPoint');
            figPos=fig.Position;
            normPos=tpos./figPos(3:4);
            ax2.Position(1:2)=-ax2.Position(3:4)./2+normPos;
        end
    end
set(fig,'currentAxes',ax1)
CMapHdl=ax2;
end

---