转载请注明本文链接
及文章作者：slandarer
MATLAB 桑基图绘制函数大更新啦，本次更新主要是根据fileexchange上大家提出和反馈的建议增添的新功能，我们首先看一下该工具的介绍，基本用法，然后再看看新功能叭：
工具介绍
SSankey 是本人编写的 MATLAB 桑基图绘制工具，可以实现高度自定义的桑基图，例如：
该工具可在以下fileexchange链接或者文末提供的gitee仓库下载：
https://www.mathworks.com/matlabcentral/fileexchange/128679-sankey-plot
基本用法
基本用法非常简单，例如：
links={'a1','A',1.2;'a2','A',1;'a1','B',.6;'a3','A',1; 'a3','C',0.5;
'b1','B',.4; 'b2','B',1;'b3','B',1; 'c1','C',1;
'c2','C',1; 'c3','C',1;'A','AA',2; 'A','BB',1.2;
'B','BB',1.5; 'B','AA',1.5; 'C','BB',2.3; 'C','AA',1.2};
% 创建桑基图对象(Create a Sankey diagram object)
SK=SSankey(links(:,1),links(:,2),links(:,3));
% 开始绘图(Start drawing)
SK.draw()
更多用法请见以下文章：
更新内容
前两条更新意见是由Kamuran Turksoy提出的：The figures look great but its difficult to use the tool if you have multiple layers with many data points. You need to define the links one by one. It would be great if it accepted data in table/matrix format and created the links automatically.
1 邻接矩阵创建桑基图
支持使用邻接矩阵创建桑基图,不预先设定节点名称的话，节点将被命名为node-n： （注意，目前不支持两个节点互相有数据流动的情况，即流动必须是单向的）
adjMat=[0,0,0,1,2,1,0,0,0,0;
0,0,0,1,2,3,0,0,0,0;
0,0,0,2,0,1,0,0,0,0;
0,0,0,0,0,0,1,4,0,0;
0,0,0,0,0,0,2,1,0,0;
0,0,0,0,0,0,0,3,0,0;
0,0,0,0,0,0,0,0,1,5;
0,0,0,0,0,0,0,0,2,3;
0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0];
nodeList=compose('C%d',1:10);
% 创建桑基图对象(Create a Sankey diagram object)
SK=SSankey([],[],[],'NodeList',nodeList,'AdjMat',adjMat);
% method 1
% SK=SSankey([],[],[],'AdjMat',adjMat);
% method 2
% SK=SSankey([],[],[],'NodeList',nodeList,'AdjMat',adjMat)
% method 3
% SK=SSankey([],[],[]);
% SK.AdjMat=adjMat;
% 开始绘图(Start drawing)
SK.draw()
2 自由调整节点层级
可以通过调整Layer属性来实现，Layer是一个数组，显示每一个节点属于哪一层，比如第i个数值是n就意味着第i个节点在n层（注意，暂不支持相互之间有流动的节点在同一层级，即层级内暂不支持数据流动）
adjMat=[0,0,0,1,2,1,0,0,0,0;
0,0,0,1,2,3,0,0,0,0;
0,0,0,2,0,1,0,0,0,0;
0,0,0,0,0,0,1,4,0,0;
0,0,0,0,0,0,2,1,0,0;
0,0,0,0,0,0,0,3,0,0;
0,0,0,0,0,0,0,0,1,5;
0,0,0,0,0,0,0,0,2,3;
0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0];
% nodeList={'C1','C2',C3',...'C10'}
nodeList=compose('C%d',1:10);
layer=[1,1,2,4,4,3,6,6,7,7];
% 创建桑基图对象(Create a Sankey diagram object)
SK=SSankey([],[],[],'NodeList',nodeList,'AdjMat',adjMat,'Layer',layer);
% SK.Layer = layer;
% 开始绘图(Start drawing)
SK.draw()
SK.moveBlockY(3,-10)
SK.moveBlockY(6,-10)
第三条更新意见是由Paul Maximilian Röhrig提出的：Is it possible to move blocks in the X direction as well as the Y direction? Many thanks for your help.
3 移动节点x坐标
使用以下方式移动节点x坐标：
obj.moveBlockX(n,dx)
links={'a1','A',1.2;'a2','A',1;'a1','B',.6;'a3','A',1; 'a3','C',0.5;
'b1','B',.4; 'b2','B',1;'b3','B',1; 'c1','C',1;
'c2','C',1; 'c3','C',1;'A','AA',2; 'A','BB',1.2;
'B','BB',1.5; 'B','AA',1.5; 'C','BB',2.3; 'C','AA',1.2};
% 创建桑基图对象(Create a Sankey diagram object)
SK=SSankey(links(:,1),links(:,2),links(:,3));
% 开始绘图(Start drawing)
SK.draw()
for i=1:9
SK.moveBlockX(i,i/10)
end
for i=10:12
SK.moveBlockX(i,(i-9)*3/10-.3)
end
新版本工具函数完整代码
classdef SSankey < handle
% Copyright (c) 2023, Zhaoxu Liu / slandarer
% =========================================================================
% @author : slandarer
% 公众号 : slandarer随笔
% 知乎 : slandarer
% =========================================================================
% # update 2.0.0(2024-02-04)
% see natureSankeyDemo1.m
%
% + 层向右对齐(Align layers to the right)
% try : obj.LayerOrder='reverse';
%
% + 单独调整每层间隙大小(Adjust the Sep size of each layer separately)
% try : obj.Sep=[.2,.06,.05,.07,.07,.08,.15];
% =========================================================================
% # update 3.0.0(2024-04-15)
% see sankeyDemo9.m sankeyDemo10.m sankeyDemo11.m
%
% + 通过邻接矩阵创建桑基图(Creating a Sankey diagram through adjacency matrix)
% method 1 :
% SK=SSankey([],[],[],'AdjMat',adjMat);
% method 2 :
% SK=SSankey([],[],[],'NodeList',nodeList,'AdjMat',adjMat)
% method 3 :
% SK=SSankey([],[],[]);
% SK.AdjMat=adjMat;
%
% try :
% adjMat=zeros(10,10);
% layerNum=[3,3,2,2];
% layerInd=cumsum([0,layerNum]);
% for i=1:length(layerInd)-2
% adjMat(layerInd(i)+1:layerInd(i+1),layerInd(i+1)+1:layerInd(i+2))=randi([1,6],[layerNum([i,i+1])]);
% end
% disp(adjMat)
% SK=SSankey([],[],[],'NodeList',nodeList,'AdjMat',adjMat);
% SK.draw()
%
% + 每层情况可被设置(Each layer state can be set)
% try : obj.Layer = [1,1,1, 2,2,2, 3,3, 4,4,...];
%
% + 每个节点可在x方向上位移(Each node can be displaced in the x-direction)
% try : obj.moveBlockX(n,dx)
properties
Source;Target;Value;
SourceInd;TargetInd;
Layer;LayerPos;LayerOrder='normal';
AdjMat;BoolMat;
RenderingMethod='interp' % 'left'/'right'/'interp'/'map'/'simple'
LabelLocation='left' % 'left'/'right'/'top'/'center'/'bottom'
Align='center' % 'up'/'down'/'center'
BlockScale=0.05; % BlockScale>0 ! !
Sep=0.05; % Sep>=0 ! !
NodeList={};
ColorList=[[65,140,240;252,180,65;224,64,10;5,100,146;191,191,191;26,59,105;255,227,130;18,156,221;
202,107,75;0,92,219;243,210,136;80,99,129;241,185,168;224,131,10;120,147,190]./255;
[127,91,93;187,128,110;197,173,143;59,71,111;104,95,126;76,103,86;112,112,124;
72,39,24;197,119,106;160,126,88;238,208,146]./255];
BlockHdl;LinkHdl;LabelHdl;ax;Parent;
BN;LN;VN;TotalLen;SepLen;
arginList={'RenderingMethod','LabelLocation','BlockScale','Layer'...
'Sep','Align','ColorList','Parent','NodeList','AdjMat'}
end
% 构造函数 =================================================================
methods
function obj=SSankey(varargin)
% 获取基本数据 -------------------------------------------------
if isa(varargin{1},'matlab.graphics.axis.Axes')
obj.ax=varargin{1};varargin(1)=[];
else
end
obj.Source=varargin{1};
obj.Target=varargin{2};
obj.Value=varargin{3};
varargin(1:3)=[];
% 获取其他信息 -------------------------------------------------
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
if isempty(obj.ax)&&(~isempty(obj.Parent)),obj.ax=obj.Parent;end
if isempty(obj.ax),obj.ax=gca;end
obj.ax.NextPlot='add';
% 基本数据预处理 -----------------------------------------------
if isempty(obj.NodeList)
if isempty(obj.Source)
if ~isempty(obj.AdjMat)
obj.NodeList=compose('node%d',1:size(obj.AdjMat,1));
end
else
obj.NodeList=[obj.Source;obj.Target];
obj.NodeList=unique(obj.NodeList,'stable');
end
end
obj.BN=length(obj.NodeList);
if length(obj.NodeList)>size(obj.ColorList,1)
obj.ColorList=[obj.ColorList;rand(length(obj.NodeList),3).*.7];
end
% obj.VN=length(obj.Value);
% 坐标区域基础设置 ---------------------------------------------
obj.ax.YDir='reverse';
obj.ax.XColor='none';
obj.ax.YColor='none';
end
% 绘图函数 =================================================================
function draw(obj)
% 生成整体邻接矩阵 ---------------------------------------------
if isempty(obj.AdjMat)
obj.AdjMat=zeros(obj.BN,obj.BN);
for i=1:length(obj.Source)
obj.SourceInd(i)=find(strcmp(obj.Source{i},obj.NodeList));
obj.TargetInd(i)=find(strcmp(obj.Target{i},obj.NodeList));
obj.AdjMat(obj.SourceInd(i),obj.TargetInd(i))=obj.Value{i};
end
end
if isempty(obj.NodeList)
obj.NodeList=compose('node%d',1:size(obj.AdjMat,1));
end;obj.BN=length(obj.NodeList);help SSankey
obj.BoolMat=abs(obj.AdjMat)>0;
if any(any(obj.BoolMat+obj.BoolMat.'==2))
warning('Currently, bidirectional flow sankey diagram plotting is not supported.')
end
obj.VN=sum(sum(obj.BoolMat));
% 计算每个对象位于的层、每层方块长度、每个方块位置 ----------------
if isempty(obj.Layer)
if strcmp(obj.LayerOrder,'normal')
obj.Layer=zeros(obj.BN,1);
obj.Layer(sum(obj.BoolMat,1)==0)=1;
startMat=diag(obj.Layer);
for i=1:(obj.BN-1)
tLayer=(sum(startMat*obj.BoolMat^i,1)>0).*(i+1);
obj.Layer=max([obj.Layer,tLayer'],[],2);
end
else
obj.Layer=zeros(obj.BN,1);
obj.Layer(sum(obj.BoolMat,2)==0)=-1;
startMat=diag(obj.Layer);
for i=1:(obj.BN-1)
tLayer=(sum(startMat*(obj.BoolMat.')^i,1)<0).*(-i-1);
obj.Layer=min([obj.Layer,tLayer'],[],2);
end
obj.Layer=obj.Layer-min(obj.Layer)+1;
end
end
obj.Layer=obj.Layer(:);
obj.LN=max(obj.Layer);
obj.TotalLen=max([sum(obj.AdjMat,1).',sum(obj.AdjMat,2)],[],2);
obj.SepLen=max(obj.TotalLen).*obj.Sep;
obj.LayerPos=zeros(obj.BN,4);
for i=1:obj.LN
tBlockInd=find(obj.Layer==i);
tBlockLen=[0;cumsum(obj.TotalLen(tBlockInd))];
tY1=tBlockLen(1:end-1)+(0:length(tBlockInd)-1).'.*obj.SepLen(min(i,length(obj.Sep)));
tY2=tBlockLen(2:end)+(0:length(tBlockInd)-1).'.*obj.SepLen(min(i,length(obj.Sep)));
obj.LayerPos(tBlockInd,3)=tY1;
obj.LayerPos(tBlockInd,4)=tY2;
% for j=1:length(tY2)
% plot([i,i],[tY1(j),tY2(j)],'LineWidth',2)
% end
end
obj.LayerPos(:,1)=obj.Layer;
obj.LayerPos(:,2)=obj.Layer+obj.BlockScale;
% 根据对齐方式调整Y坐标 -----------------------------------------
tMinY=min(obj.LayerPos(:,3));
tMaxY=max(obj.LayerPos(:,4));
for i=1:obj.LN
tBlockInd=find(obj.Layer==i);
tBlockPos3=obj.LayerPos(tBlockInd,3);
tBlockPos4=obj.LayerPos(tBlockInd,4);
switch obj.Align
case 'up'
case 'down'
obj.LayerPos(tBlockInd,3)=obj.LayerPos(tBlockInd,3)+tMaxY-max(tBlockPos4);
obj.LayerPos(tBlockInd,4)=obj.LayerPos(tBlockInd,4)+tMaxY-max(tBlockPos4);
case 'center'
obj.LayerPos(tBlockInd,3)=obj.LayerPos(tBlockInd,3)+...
min(tBlockPos3)/2-max(tBlockPos4)/2+tMinY/2-tMaxY/2;
obj.LayerPos(tBlockInd,4)=obj.LayerPos(tBlockInd,4)+...
min(tBlockPos3)/2-max(tBlockPos4)/2+tMinY/2-tMaxY/2;
end
end
% 绘制连接 -----------------------------------------------------
if isempty(obj.SourceInd)
[obj.SourceInd,obj.TargetInd]=find(obj.AdjMat~=0);
end
for i=1:obj.VN
tSource=obj.SourceInd(i);
tTarget=obj.TargetInd(i);
tS1=sum(obj.AdjMat(tSource,1:(tTarget-1)))+obj.LayerPos(tSource,3);
tS2=sum(obj.AdjMat(tSource,1:tTarget))+obj.LayerPos(tSource,3);
tT1=sum(obj.AdjMat(1:(tSource-1),tTarget))+obj.LayerPos(tTarget,3);
tT2=sum(obj.AdjMat(1:tSource,tTarget))+obj.LayerPos(tTarget,3);
if isempty(tS1),tS1=0;end
if isempty(tT1),tT1=0;end
tX=[obj.LayerPos(tSource,1),obj.LayerPos(tSource,2),obj.LayerPos(tTarget,1),obj.LayerPos(tTarget,2)];
if abs(tX(1)-tX(3))<eps
warning('Currently, flow between the same layer is not supported.')
break
end
qX=linspace(obj.LayerPos(tSource,1),obj.LayerPos(tTarget,2),200);qT=linspace(0,1,50);
qY1=interp1(tX,[tS1,tS1,tT1,tT1],qX,'pchip');
qY2=interp1(tX,[tS2,tS2,tT2,tT2],qX,'pchip');
XX=repmat(qX,[50,1]);YY=qY1.*(qT'.*0+1)+(qY2-qY1).*(qT');
MeshC=ones(50,200,3);
switch obj.RenderingMethod
case 'left'
MeshC(:,:,1)=MeshC(:,:,1).*obj.ColorList(tSource,1);
MeshC(:,:,2)=MeshC(:,:,2).*obj.ColorList(tSource,2);
MeshC(:,:,3)=MeshC(:,:,3).*obj.ColorList(tSource,3);
case 'right'
MeshC(:,:,1)=MeshC(:,:,1).*obj.ColorList(tTarget,1);
MeshC(:,:,2)=MeshC(:,:,2).*obj.ColorList(tTarget,2);
MeshC(:,:,3)=MeshC(:,:,3).*obj.ColorList(tTarget,3);
case 'interp'
MeshC(:,:,1)=repmat(linspace(obj.ColorList(tSource,1),obj.ColorList(tTarget,1),200),[50,1]);
MeshC(:,:,2)=repmat(linspace(obj.ColorList(tSource,2),obj.ColorList(tTarget,2),200),[50,1]);
MeshC(:,:,3)=repmat(linspace(obj.ColorList(tSource,3),obj.ColorList(tTarget,3),200),[50,1]);
case 'map'
MeshC=MeshC(:,:,1).*obj.Value{i};
case 'simple'
MeshC(:,:,1)=MeshC(:,:,1).*.6;
MeshC(:,:,2)=MeshC(:,:,2).*.6;
MeshC(:,:,3)=MeshC(:,:,3).*.6;
end
obj.LinkHdl(i)=surf(obj.ax,XX,YY,XX.*0,'EdgeColor','none','FaceAlpha',.3,'CData',MeshC);
end
% 绘制方块 -----------------------------------------------------
for i=1:obj.BN
obj.BlockHdl(i)=fill(obj.ax,obj.LayerPos(i,[1,2,2,1]),...
obj.LayerPos(i,[3,3,4,4]),obj.ColorList(i,:),'EdgeColor','none');
end
% 绘制文本 -----------------------------------------------------
for i=1:obj.BN
switch obj.LabelLocation
case 'right'
obj.LabelHdl(i)=text(obj.ax,obj.LayerPos(i,2),mean(obj.LayerPos(i,[3,4])),...
[' ',obj.NodeList{i}],'FontSize',15,'FontName','Times New Roman','HorizontalAlignment','left');
case 'left'
obj.LabelHdl(i)=text(obj.ax,obj.LayerPos(i,1),mean(obj.LayerPos(i,[3,4])),...
[obj.NodeList{i},' '],'FontSize',15,'FontName','Times New Roman','HorizontalAlignment','right');
case 'top'
obj.LabelHdl(i)=text(obj.ax,mean(obj.LayerPos(i,[1,2])),obj.LayerPos(i,3),...
obj.NodeList{i},'FontSize',15,'FontName','Times New Roman','HorizontalAlignment','center','VerticalAlignment','bottom');
case 'center'
obj.LabelHdl(i)=text(obj.ax,mean(obj.LayerPos(i,[1,2])),mean(obj.LayerPos(i,[3,4])),...
obj.NodeList{i},'FontSize',15,'FontName','Times New Roman','HorizontalAlignment','center');
case 'bottom'
obj.LabelHdl(i)=text(obj.ax,mean(obj.LayerPos(i,[1,2])),obj.LayerPos(i,4),...
obj.NodeList{i},'FontSize',15,'FontName','Times New Roman','HorizontalAlignment','center','VerticalAlignment','top');
end
end
% -------------------------------------------------------------
axis tight;
end
% =========================================================================
function setBlock(obj,n,varargin)
set(obj.BlockHdl(n),varargin{:})
end
function setLink(obj,n,varargin)
set(obj.LinkHdl(n),varargin{:})
end
function setLabel(obj,n,varargin)
set(obj.LabelHdl(n),varargin{:})
end
function moveBlockX(obj,n,dx)
obj.LayerPos(n,[1,2])=obj.LayerPos(n,[1,2])+dx;
set(obj.BlockHdl(n),'XData',obj.LayerPos(n,[1,2,2,1]));
switch obj.LabelLocation
case 'right',set(obj.LabelHdl(n),'Position',[obj.LayerPos(n,2),mean(obj.LayerPos(n,[3,4]))]);
case 'left',set(obj.LabelHdl(n),'Position',[obj.LayerPos(n,1),mean(obj.LayerPos(n,[3,4]))]);
case 'top',set(obj.LabelHdl(n),'Position',[mean(obj.LayerPos(n,[1,2])),obj.LayerPos(n,3)]);
case 'center',set(obj.LabelHdl(n),'Position',[mean(obj.LayerPos(n,[1,2])),mean(obj.LayerPos(n,[3,4]))]);
case 'bottom',set(obj.LabelHdl(n),'Position',[mean(obj.LayerPos(n,[1,2])),obj.LayerPos(n,4)]);
end
for i=1:obj.VN
tSource=obj.SourceInd(i);
tTarget=obj.TargetInd(i);
if tSource==n||tTarget==n
tS1=sum(obj.AdjMat(tSource,1:(tTarget-1)))+obj.LayerPos(tSource,3);
tS2=sum(obj.AdjMat(tSource,1:tTarget))+obj.LayerPos(tSource,3);
tT1=sum(obj.AdjMat(1:(tSource-1),tTarget))+obj.LayerPos(tTarget,3);
tT2=sum(obj.AdjMat(1:tSource,tTarget))+obj.LayerPos(tTarget,3);
if isempty(tS1),tS1=0;end
if isempty(tT1),tT1=0;end
tX=[obj.LayerPos(tSource,1),obj.LayerPos(tSource,2),obj.LayerPos(tTarget,1),obj.LayerPos(tTarget,2)];
qX=linspace(obj.LayerPos(tSource,1),obj.LayerPos(tTarget,2),200);qT=linspace(0,1,50);
qY1=interp1(tX,[tS1,tS1,tT1,tT1],qX,'pchip');
qY2=interp1(tX,[tS2,tS2,tT2,tT2],qX,'pchip');
YY=qY1.*(qT'.*0+1)+(qY2-qY1).*(qT');
set(obj.LinkHdl(i),'YData',YY,'XData',qX);
end
end
end
function moveBlockY(obj,n,dy)
obj.LayerPos(n,[3,4])=obj.LayerPos(n,[3,4])-dy;
set(obj.BlockHdl(n),'YData',obj.LayerPos(n,[3,3,4,4]));
switch obj.LabelLocation
case 'right',set(obj.LabelHdl(n),'Position',[obj.LayerPos(n,2),mean(obj.LayerPos(n,[3,4]))]);
case 'left',set(obj.LabelHdl(n),'Position',[obj.LayerPos(n,1),mean(obj.LayerPos(n,[3,4]))]);
case 'top',set(obj.LabelHdl(n),'Position',[mean(obj.LayerPos(n,[1,2])),obj.LayerPos(n,3)]);
case 'center',set(obj.LabelHdl(n),'Position',[mean(obj.LayerPos(n,[1,2])),mean(obj.LayerPos(n,[3,4]))]);
case 'bottom',set(obj.LabelHdl(n),'Position',[mean(obj.LayerPos(n,[1,2])),obj.LayerPos(n,4)]);
end
for i=1:obj.VN
tSource=obj.SourceInd(i);
tTarget=obj.TargetInd(i);
if tSource==n||tTarget==n
tS1=sum(obj.AdjMat(tSource,1:(tTarget-1)))+obj.LayerPos(tSource,3);
tS2=sum(obj.AdjMat(tSource,1:tTarget))+obj.LayerPos(tSource,3);
tT1=sum(obj.AdjMat(1:(tSource-1),tTarget))+obj.LayerPos(tTarget,3);
tT2=sum(obj.AdjMat(1:tSource,tTarget))+obj.LayerPos(tTarget,3);
if isempty(tS1),tS1=0;end
if isempty(tT1),tT1=0;end
tX=[obj.LayerPos(tSource,1),obj.LayerPos(tSource,2),obj.LayerPos(tTarget,1),obj.LayerPos(tTarget,2)];
qX=linspace(obj.LayerPos(tSource,1),obj.LayerPos(tTarget,2),200);qT=linspace(0,1,50);
qY1=interp1(tX,[tS1,tS1,tT1,tT1],qX,'pchip');
qY2=interp1(tX,[tS2,tS2,tT2,tT2],qX,'pchip');
YY=qY1.*(qT'.*0+1)+(qY2-qY1).*(qT');
set(obj.LinkHdl(i),'YData',YY);
end
end
end
end
% Copyright (c) 2023, Zhaoxu Liu / slandarer
% =========================================================================
% @author : slandarer
% 公众号 : slandarer随笔
% 知乎 : slandarer
% -------------------------------------------------------------------------
end
完
以上已经是完整代码，未经允许本代码请勿作商业用途，引用的话可以引用我file exchange上的链接，可使用如下格式：
Zhaoxu Liu / slandarer (2023). sankey plot (https://www.mathworks.com/matlabcentral/fileexchange/128679-sankey-plot), MATLAB Central File Exchange. 检索来源 2023/4/28.
若转载请保留以上file exchange链接及本文链接！！！！!
该工具可通过上述fileexchange链接获取，或者通过以下gitee仓库下载：
https://gitee.com/slandarer/matlab-sankey-diagram
