【MATLAB】地図の色付けをやってみた

# はじめに


昨今、赤色や青色のアメリカ地図をよくニュースで見かけますが、MATLABでそのような地図をどう作るのでしょうか。政治がらみの記事にはしたくないので、私の過去を振り返るネタで書きたいと思います。

以前、アメリカの本社でアプリケーションエンジニアをしていた頃、出張で全国を飛び回っていました。実際どの州を訪問したか可視化してみることにしました。

![image_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/194780/e5726e5f-2499-25b8-8ca6-ab9cd36ca32e.png)


今回はこのアニメーションの作成方法について解説します。


# Mapping Toolbox


地図上に色付きのパッチを表示するために Mapping Toolbox を活用しました。主に使った機能はこちらの３つ。



   -  地図データ（州の境界線の地図座標）をシェープファイルから読み込む `shaperead` 関数 
   -  アメリカのMap axesを作成する `usamap` 関数（世界の場合 `worldmap` 関数を使う） 
   -  色付きパッチを表示する `geoshow` 関数 



あとは、出張の訪問履歴から州情報を抜き出し、対応する州を色付けし、ループでグラフィックスをアップデートします。


# 出張データ

```matlab:Code
visitData = readtable("VisitedStates.xlsx","Sheet","Sheet2","TextType","string");
```



こんな感じでExcelデータを用意しました。



```matlab:Code
visitData(1:10,:)
```

| |Date|State|
|:--:|:--:|:--:|
|1|2012/01/18|"Massachusetts"|
|2|2012/02/15|"Tennessee"|
|3|2012/03/05|"Pennsylvania"|
|4|2012/03/07|"Pennsylvania"|
|5|2012/03/09|"Canada"|
|6|2012/03/21|"Pennsylvania"|
|7|2012/03/27|"Virginia"|
|8|2012/04/19|"Ohio"|
|9|2012/05/02|"Illinois"|
|10|2012/05/21|"Washington"|



後々、処理しやすいように、日付を月始めに調整します。



```matlab:Code
visitData.Date = dateshift(visitData.Date,'start','month');
visitData(1:10,:)
```

| |Date|State|
|:--:|:--:|:--:|
|1|2012/01/01|"Massachusetts"|
|2|2012/02/01|"Tennessee"|
|3|2012/03/01|"Pennsylvania"|
|4|2012/03/01|"Pennsylvania"|
|5|2012/03/01|"Canada"|
|6|2012/03/01|"Pennsylvania"|
|7|2012/03/01|"Virginia"|
|8|2012/04/01|"Ohio"|
|9|2012/05/01|"Illinois"|
|10|2012/05/01|"Washington"|

# 地図データの準備


地図データの描画はほぼ Mapping Toolbox で行いました。




まずはアメリカの州の境界線データをMapping Toolboxに含まれているシェープファイルから読み取ります。usastatelo.shp にはアメリカの50州+Washington D.C.の境界線の座標データ（緯度、経度）が含まれています。



```matlab:Code
states = shaperead("usastatelo","UseGeoCoords",true);

% 州名を抽出
names = string({states.Name});

% アラスカとハワイは別軸に表示されるため、それぞれのインデックスを取得
indexAlaska = names == "Alaska";
indexHawaii = names == "Hawaii";
```



次に、使う色を指定します。まだ訪問していない州、訪問済みの州、今回訪問した州、の３種類の色をしていします。



```matlab:Code
haveBeenColor = [.5 1 .5];
justHaveBeenColor = [0 .5 0];
haveNotBeenColor = [.75 .75 .75];
```



まずアメリカの map axes を用意するのですが、Toolboxにある `usamap` 関数を使います。



```matlab:Code
hFig = figure("Units","pixels","Position",[100 100 800 600],...
  "Color","white");
hAx = usamap("all");
```

![figure_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/194780/f06cb5df-b7d6-597f-216f-83b2945afcec.png)


# 地図の描画


次に `geoshow` 関数で州の境界線を描きます。最初の軸にはハワイとアラスカ以外の州を、2つ目の軸にはアラスカを、3つ目の軸にはハワイを描きます。



```matlab:Code
hMain = geoshow(hAx(1),states(~indexHawaii & ~indexAlaska),"FaceColor",haveNotBeenColor);
hAlaska = geoshow(hAx(2),states(indexAlaska),"FaceColor",haveNotBeenColor);
hHawaii = geoshow(hAx(3),states(indexHawaii),"FaceColor",haveNotBeenColor);

hTitle = title(hAx(1),"","FontSize",16);

% 訪問していない州のリスト
hStateList = axes(hFig,"Units","normalized","Position",[0.8 0 0.2 1],"Visible","off");
hStateListText = text(0,1,["訪問していない州（" + length(names) + "）：",names], ...
                      "VerticalAlignment","top");
```


![figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/194780/fd16a5d0-6a57-89ab-5b54-778b74f39ecf.png)




hMainを見てみると、オブジェクトの Children に49のパッチオブジェクトがあるのが分かります。



```matlab:Code
hMain
```


```text:Output
hMain = 
  Group のプロパティ:

    Children: [49x1 Patch]
     Visible: on
     HitTest: on

  すべてのプロパティ を表示
```



つまり49州のパッチオブジェクトです。




map axes の枠を消しましょう。



```matlab:Code
for k = 1:length(hAx)
  setm(hAx(k),"Frame","off","Grid","off",...
    "ParallelLabel","off","MeridianLabel","off")
end
```

![figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/194780/e03589fc-f0fb-b642-674a-5ffd375bdd88.png)


# 訪問した州（訪問してない州）


準備完了です！あとは、`for` ループで月ごとに地図をアップデートします。




試しに、ある月で見てましょう。例えば2010年10月時点での訪問データを可視化してみます。



```matlab:Code
cutoffDate = datetime(2010,10,1);
```

### 1. 出張データから、ある日付（月）までに訪問した州名を抽出。また、その月に訪問した州名も抽出
```matlab:Code
goneToStates = unique(visitData.State(visitData.Date < cutoffDate));
justGoneToStates = unique(visitData.State(visitData.Date == cutoffDate));
```

### 2. 抽出した州名のインデックスを取得
```matlab:Code
haveBeenIDX = contains(names, goneToStates);
justHaveBeenIDX = contains(names, justGoneToStates);
haveNotBeenIDX = ~(haveBeenIDX | justHaveBeenIDX);
```

### 3. 州に適切な色で色付け
```matlab:Code
% 訪問していない州
set(hMain.Children(haveNotBeenIDX(~indexHawaii & ~indexAlaska)),...
    "FaceColor",haveNotBeenColor)
set(hAlaska.Children(haveNotBeenIDX(indexAlaska)),"FaceColor",haveNotBeenColor)
set(hHawaii.Children(haveNotBeenIDX(indexHawaii)),"FaceColor",haveNotBeenColor)

% 訪問した州
set(hMain.Children(haveBeenIDX(~indexHawaii & ~indexAlaska)),...
    "FaceColor",haveBeenColor)
set(hAlaska.Children(haveBeenIDX(indexAlaska)),"FaceColor",haveBeenColor)
set(hHawaii.Children(haveBeenIDX(indexHawaii)),"FaceColor",haveBeenColor)

% この月に訪問した州
set(hMain.Children(justHaveBeenIDX(~indexHawaii & ~indexAlaska)),...
    "FaceColor",justHaveBeenColor)
set(hAlaska.Children(justHaveBeenIDX(indexAlaska)),"FaceColor",justHaveBeenColor)
set(hHawaii.Children(justHaveBeenIDX(indexHawaii)),"FaceColor",justHaveBeenColor)
```

### 4. タイトルと訪問していない州リストの更新
```matlab:Code
hTitle.String = char(cutoffDate,"uuuu年MM月");
hStateListText.String = ["訪問していない州（" + nnz(haveNotBeenIDX) + ...
                         "）：",names(haveNotBeenIDX)];
```

![figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/194780/aac41c57-99d7-902f-b71c-83f78ef72764.png)


# 完成版


全てをくっつけるとこんな感じになります。



```matlab:Code(Display)
% 出張データ読み込み
visitData = readtable("VisitedStates.xlsx","Sheet","Sheet2","TextType","string");
visitData.Date = dateshift(visitData.Date,'start','month');

% アメリカの州の境界線データ読み込み
states = shaperead("usastatelo","UseGeoCoords",true);

% 州名を抽出
names = string({states.Name});

% アラスカとハワイは別軸に表示されるため、それぞれのインデックスを取得
indexAlaska = names == "Alaska";
indexHawaii = names == "Hawaii";

% 色の定義
haveBeenColor = [.5 1 .5];
justHaveBeenColor = [0 .5 0];
haveNotBeenColor = [.75 .75 .75];

% Map Axesの準備
hFig = figure("Units","pixels","Position",[100 100 800 600],...
  "Color","white");
hAx = usamap("all");

% 州の境界線
hMain = geoshow(hAx(1),states(~indexHawaii & ~indexAlaska),"FaceColor",haveNotBeenColor);
hAlaska = geoshow(hAx(2),states(indexAlaska),"FaceColor",haveNotBeenColor);
hHawaii = geoshow(hAx(3),states(indexHawaii),"FaceColor",haveNotBeenColor);

hTitle = title(hAx(1),"","FontSize",16);

% 訪問していない州のリスト
hStateList = axes(hFig,"Units","normalized","Position",[0.8 0 0.2 1],"Visible","off");
hStateListText = text(0,1,["訪問していない州（"+length(names)+"）：",names],"VerticalAlignment","top");

% Map Axesの枠削除
for k = 1:length(hAx)
  setm(hAx(k),"Frame","off","Grid","off",...
    "ParallelLabel","off","MeridianLabel","off")
end

% 月ベクトルの定義
monthlist = datetime(2006,5,1):calmonths(1):datetime(2013,12,1);

for idx = 1:length(monthlist)
    cutoffDate = monthlist(idx);
    
    % 出張データから、ある日付（月）までに訪問した州名を抽出。また、その月に訪問した州名も抽出
    goneToStates = unique(visitData.State(visitData.Date < cutoffDate));      % この月までに訪問した州名
    justGoneToStates = unique(visitData.State(visitData.Date == cutoffDate)); % この月に訪問した州名
    
    % 抽出した州名のインデックスを取得
    haveBeenIDX = contains(names, goneToStates);
    justHaveBeenIDX = contains(names, justGoneToStates);
    haveNotBeenIDX = ~(haveBeenIDX | justHaveBeenIDX);
    
    % 州に適切な色で色付け
    % 訪問していない州
    set(hMain.Children(haveNotBeenIDX(~indexHawaii & ~indexAlaska)),"FaceColor",haveNotBeenColor)
    set(hAlaska.Children(haveNotBeenIDX(indexAlaska)),"FaceColor",haveNotBeenColor)
    set(hHawaii.Children(haveNotBeenIDX(indexHawaii)),"FaceColor",haveNotBeenColor)
    
    % 訪問した州
    set(hMain.Children(haveBeenIDX(~indexHawaii & ~indexAlaska)),"FaceColor",haveBeenColor)
    set(hAlaska.Children(haveBeenIDX(indexAlaska)),"FaceColor",haveBeenColor)
    set(hHawaii.Children(haveBeenIDX(indexHawaii)),"FaceColor",haveBeenColor)
    
    % この月に訪問した州
    set(hMain.Children(justHaveBeenIDX(~indexHawaii & ~indexAlaska)),"FaceColor",justHaveBeenColor)
    set(hAlaska.Children(justHaveBeenIDX(indexAlaska)),"FaceColor",justHaveBeenColor)
    set(hHawaii.Children(justHaveBeenIDX(indexHawaii)),"FaceColor",justHaveBeenColor)
    
    % タイトルと訪問していない州リストの更新
    hTitle.String = char(cutoffDate,"uuuu年MM月");
    hStateListText.String = ["訪問していない州（" + nnz(haveNotBeenIDX) + "）：",names(haveNotBeenIDX)];
    
    pause(0.1)
end
```



最後にもう一度アニメーションを見てみましょう。


![image_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/194780/0741d82f-a5c1-edf1-4a05-bbc9310ae54a.png)

これで、データさえあればアメリカの大統領選挙の結果も簡単に可視化できますね！