MATLAB | 全网唯一,三元相映射图绘制(三元相分级统计地图)

示例

---

对上期的补充

本期文章名叫三元相映射图绘制，但不单单讲这一点，同时会对上一篇多元映射地图做出补充https://slandarer.blog.csdn.net/article/details/126702967：

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

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

代码中出现的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;
1
2
3
4
5
6
7
8
9
10
11
12
13

---

三元相映射图绘制

使用基本说明

三元相映射图绘制函数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;
1
2
3
4
5
6
7
8
9
10
11
12
13

图例可拖拽移动：

配色展示

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;
1
2
3
4
5
6
7
8
9
10
11
12
13

全部配色展示：

示例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 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

示例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)']; 
1
2
3
4
5
6
7
8
9
10
11
12
13

---

三元相映射图完整代码

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163

---

完

如果有用请留个赞叭~

未经允许本代码请勿作商业用途，引用的话可以引用我file exchange上的链接，可使用如下格式：

Zhaoxu Liu (2022). Ternary phase map 三元相映射图 (https://www.mathworks.com/matlabcentral/fileexchange/117195-ternary-phase-map), MATLAB Central File Exchange. 检索来源 2022/9/5.

若转载请保留以上file exchange链接及本文链接！！！

全部m文件及素材包获取:
链接：https://pan.baidu.com/s/1s_hr5RHeFsJ9RCv-pP7NlQ?pwd=slan
提取码：slan
–来自百度网盘超级会员V5的分享