クラスタリングしたデータの地図へのマッピング

添付させて頂いたテキストはそれぞれ地球全体の地表の一年間の平均気温(ave) 平均変化度(var)です。中のデータは地球全体のマップ2.5度間隔で 72 x 144 array　となっています。 そのデータを以下のプログラムでクラスタリングしたのですが、その後そのクラスタリングデータをその分類された色のままマップ上にマッピングしたいのですが何か方法はあるでしょうか。アドバイス頂けると幸いです 以下クラスタリングのプログラムです

xfilename =('TS_ave.txt');
x = load(xfilename);
yfilename =('TS_var.txt');
y = load(yfilename);
a=15;
rng default;
X = [x;y];
X = reshape(X,10368,2);
NX = zscore(X);
opts = statset('Display','final');
[idx,C] = kmeans(NX,3,'Distance','cityblock','Replicates',5,'Options',opts);
plot(NX(idx==1,1),NX(idx==1,2),'r.','MarkerSize',12);
hold on;
plot(NX(idx==2,1),NX(idx==2,2),'b.','MarkerSize',12);
plot(NX(idx==3,1),NX(idx==3,2),'g.','MarkerSize',12);
%plot(NX(idx==4,1),NX(idx==4,2),'k.','MarkerSize',12);
plot(C(:,1),C(:,2),'kx','MarkerSize',15,'LineWidth',3);
legend('Cluster 1','Cluster 2','Cluster 3','Centroids','Location','NW');
title 'Cluster Assignments and Centroids';
xlabel('average temperature');
ylabel('temperature variations');
hold off;

以下データのマッピングプログラム

filename =('TS_ave.txt'); 
%filename =('TS_var.txt');
x = load(filename);
x = reshape(x,72,144);
rv = [0.4 90 180]; 
ax = worldmap('World');
%worldmap('world');
geoshow(x, rv, 'displaytype', 'texturemap');
%C = load('coast');
%plotm(C.lat, C.long, 'k');
setm(ax, 'Origin', [0 0 0])
land = shaperead('landareas', 'UseGeoCoords', true);
geoshow(ax, land, 'FaceColor', [0.5 0.7 0.5], 'FaceAlpha',[0]);

Answer 1

[michio]

クラスタリング結果 idx をインデックス画像として取扱う方法はいかがでしょうか。

 geoshow(idx,cmap,rv);

kmeans 関数の出力の１つ idx にはクラスタリング後のクラス番号が出力されています。cmap で番号に対応する色を直接指定することができます。

   cmap = [1,0,0; 0,0,1; 0,1,0];

n行目が数値nの色に対応します。 すなわち、 idx を画像として扱い、要素が１となっている部分は赤色、２となっている部分は青、３となっている部分は緑と表示させることになります。下のサンプルコードでは cmap で色の対応を定義しています。

• [1,0,0] => 赤
• [0,0,1] => 青
• [0,1,0] => 緑　

rng default;
% データの読み込み
xfilename =('TS_ave.txt');
yfilename =('TS_var.txt');
x = load(xfilename);
y = load(yfilename);

% データの正規化
X = [x,y];
NX = zscore(X);
opts = statset('Display','final');

% クラスタリング
[idx,C] = kmeans(NX,3,'Distance','cityblock','Replicates',5,'Options',opts);
idx = reshape(idx,72,144);
rv = [0.4 90 180]; 
ax = worldmap('World');
cmap = [1,0,0; 0,0,1; 0,1,0]; % 色指定
% 1 for [1,0,0] = red;
% 2 for [0,0,1] = blue;
% 3 for [0,1,0] = green
geoshow(idx,cmap,rv); % 色分けして表示

setm(ax, 'Origin', [0 0 0])
land = shaperead('landareas', 'UseGeoCoords', true);
geoshow(ax, land, 'FaceColor', [0.5 0.7 0.5], 'FaceAlpha',[0]);