Pythonで世界銀行のオープンデータを使ってバブルチャートを作成する

世界銀行が公開しているオープンデータ（World Bank Open Data）があって、データ取得用のAPIも用意されている。Pythonでこのオープンデータを取得してバブルチャートを作成してみる。

環境

WSL2(Ubuntu20.04)とJupyter Lab。
$ lsb_release -dr Description: Ubuntu 20.04.3 LTS Release: 20.04 $ python3 -V Python 3.8.10

APIで世界銀行のオープンデータを取得

オープンデータの取得にはWBGAPIを使う。WBGAPIでは世界銀行のAPIを利用してPythonでオープンデータにアクセスできる。インストールはpipでできる。

$ pip3 install wbgapi

世界銀行のオープンデータには人口や所得、消費電力などいろいろなデータがある。WBGAPIで「power」という語を含むデータを検索してみる。

import wbgapi as wb

display(wb.series.info(q='power'))

以下のような結果が得られる。データを取得するときにはこの結果にあるidを使う。

次にオープンデータに含まれる国や地域を確認してみる。オープンデータには国別データの他にいくつかの国をまとめた地域ごとのデータもある。

import wbgapi as wb

display(wb.economy.info())

次のような国IDや地域などの情報が確認できる。

先の結果にあるregionの情報は以下のようにして確認できる。

import wbgapi as wb

display(wb.region.info())

データを取得するときにはデータの種類、国または地域、年を指定できる。以下では2019年のアメリカの人口を取得できる。

import wbgapi as wb

display(wb.data.DataFrame(['SP.POP.TOTL'], ['USA'], [2019]))

オープンデータをもとにバブルチャートを作成する

WBGAPIで取得したデータをもとにバブルチャートを作成してみる。まずは必要なPythonライブラリのインストール。

$ pip3 install pandas matplotlib seaborn

続いて4種類のデータを取得してバブルチャートを作成する。すべての国だと数が多すぎて見づらくなるので人口上位50にしぼってある。使用した4種類のデータは以下の通り。

• SP.POP.TOTL: 人口
• NY.GDP.PCAP.CD: ひとりあたりGDP
• SP.DYN.LE00.IN: 出生時平均寿命
• BN.GSR.FCTY.CD: 一次純所得

%matplotlib inline

import pandas as pd
import wbgapi as wb
import matplotlib.pyplot as plt
import seaborn as sns

# 地域を除いた国のリストを作成
# regionが空白になっていないidのみにする 
economy_info = wb.economy.info().items
countries = [ei['id'] for ei in economy_info if ei['region'] != '']

# 2019年の国別の以下のデータを取得する
# SP.POP.TOTL: 人口
# NY.GDP.PCAP.CD: ひとりあたりGDP
# SP.DYN.LE00.IN: 出生時平均寿命
# BN.GSR.FCTY.CD: 一次純所得
df=wb.data.DataFrame(
    ['SP.POP.TOTL', 'NY.GDP.PCAP.CD', 'SP.DYN.LE00.IN', 'BN.GSR.FCTY.CD'],
    countries,
    [2019]
)

# データがない国を省いて人口上位50にしぼる
df = df.dropna().sort_values('SP.POP.TOTL', ascending=False).head(50)

plt.figure(figsize=(18, 12))

# バブルチャート作成
# NY.GDP.PCAP.CD: Y軸
# SP.DYN.LE00.IN: X軸
# SP.POP.TOTL: 円の大きさ
# BN.GSR.FCTY.CD: 色
ax = sns.scatterplot(
    data=df, 
    x='NY.GDP.PCAP.CD', 
    y='SP.DYN.LE00.IN', 
    size='SP.POP.TOTL', 
    c=df['BN.GSR.FCTY.CD'], 
    cmap='autumn_r', 
    alpha=0.9, 
    edgecolors='grey', 
    linewidth=1, 
    legend=False, 
    sizes=(20, 5000)
)

# 各円に国IDを表示する
def put_label(x, y, labels, ax):
    a = pd.concat({'x': x, 'y': y, 'label': labels}, axis=1)
    for i, point in a.iterrows():
        ax.text(point['x']+.02, point['y'], str(point['label']))

put_label(df['NY.GDP.PCAP.CD'], df['SP.DYN.LE00.IN'], df.index.to_series(), plt.gca())

plt.show()

次のようなバブルチャートが作成された。

PythonのMatplotlibで矢印つきの注釈を表示する

Pythonのmatplotlibでラベル付き散布図を作成するでは散布図の各要素にラベル（注釈）を表示させたが、Matplotlibでは矢印つきの注釈をグラフ上に表示できる。今回は折れ線グラフ上に矢印つきの注釈を表示させてみる。

環境

Debian Buster（Dockerコンテナ）とJupyter Lab。さらに日本語フォント（Takao）もインストール済み。
# cat /etc/issue Debian GNU/Linux 10 \n \l # python -V Python 3.7.6 # pip list | grep -e pandas -e matplotlib -e jupyter jupyter 1.0.0 jupyter-client 5.3.4 jupyter-console 6.1.0 jupyter-core 4.6.1 jupyterlab 1.2.6 jupyterlab-server 1.0.6 matplotlib 3.1.3 pandas 1.0.1 # apt list --installed | grep fonts-takao fonts-takao-gothic/stable,now 00303.01-2 all [installed,automatic] fonts-takao-mincho/stable,now 00303.01-2 all [installed,automatic] fonts-takao/stable,now 00303.01-2 all [installed]

グラフで表示するデータ

まずはグラフ表示するデータを用意する。DTW（動的時間伸縮法）で時系列データ間の距離を求めるで作成したcsvを使う。元データは、厚生労働省のホームページ地域ごとの感染状況等の公表についてからダウンロードできる「３．帰国者・接触者相談センターへの相談件数の推移（都道府県別・各日）・帰国者・接触者外来の受診者数の推移・うちPCR検査実施件数の推移（都道府県別・各日） 」のPDFで、これをcsvに変換したもの。全国と各都道府県の日ごとの数値が含まれる。以下はExcelで開いた状態。

注釈のデータ

今回は主なニュースを注釈としてグラフ上に表示させる。グラフで表示するデータが日ごとの数値なので、日にちをキーとするPythonの辞書型でデータを用意する。

from datetime import datetime as dt
# グラフ上に表示するイベント
events = {
      dt(2020, 5, 2): '国内死者500人超え',
      dt(2020, 5, 4): '緊急事態宣言を5月31日まで延長',
      dt(2020, 5, 14): '39県で緊急事態宣言解除を決定',
      dt(2020, 5, 21): '京都、大阪、兵庫の緊急事態宣言解除を決定',
      dt(2020, 5, 25): '北海道および首都圏の1都3県緊急事態宣言解除を決定',    
}

グラフ用データの読み込み

PDFから抽出したデータを保存したcsv（000643758.csv）を読み込んで、このcsvに含まれる全国の「帰国者・接触者相談センター(全相談件数)」を使う。期間は５月のひと月分。日ごとのデータなので、日をインデックスとするPandasのDataFrameにする。

%matplotlib inline
 
from datetime import datetime as dt
import pandas as pd
import matplotlib.pyplot as plt

# 日本語フォントの指定
plt.rcParams['font.family'] = 'TakaoGothic'
 
DT_FROM = '2020-05-01'
DT_TO = '2020-05-30'

def load_csv(csvpath):
    df = pd.read_csv(csvpath, index_col=0, header=[0,1], parse_dates=True, encoding='shift-jis')
  
    # 全国の「帰国者・接触者相談センター(全相談件数)」だけを選択
    df_con = df.loc[:, [('全国', '帰国者・接触者相談センター(全相談件数)')]]
     
    # カラムのレベルを全国のみにする
    df_con.columns = df_con.columns.droplevel(level=1)
 
    return df_con[DT_FROM:DT_TO]
 
df = load_csv('000643758.csv')
df.plot()
plt.show()

プロットすると以下のようになる。

グラフ上に矢印付きの注釈を表示

はじめに折れ線グラフを表示し、さらにその上に辞書型の注釈データを使って矢印つき注釈を表示する。

%matplotlib inline
 
from datetime import datetime as dt
import pandas as pd
import matplotlib.pyplot as plt

# 日本語フォントの指定
plt.rcParams['font.family'] = 'TakaoGothic'
 
DT_FROM = '2020-05-01'
DT_TO = '2020-05-30'

# グラフ上に表示するイベント
events = {
      dt(2020, 5, 2): '国内死者500人超え',
      dt(2020, 5, 4): '緊急事態宣言を5月31日まで延長',
      dt(2020, 5, 14): '39県で緊急事態宣言解除を決定',
      dt(2020, 5, 21): '京都、大阪、兵庫の緊急事態宣言解除を決定',
      dt(2020, 5, 25): '北海道および首都圏の1都3県緊急事態宣言解除を決定',    
}

def load_csv(csvpath):
    df = pd.read_csv(csvpath, index_col=0, header=[0,1], parse_dates=True, encoding='shift-jis')
  
    # 全国の「帰国者・接触者相談センター(全相談件数)」だけを選択
    df_con = df.loc[:, [('全国', '帰国者・接触者相談センター(全相談件数)')]]
     
    # カラムのレベルを全国のみにする
    df_con.columns = df_con.columns.droplevel(level=1)
 
    return df_con[DT_FROM:DT_TO]
 
df = load_csv('000643758.csv')
df.plot()

# y軸の最大値
mx = df.max()

# イベントを注釈としてグラフに追加
for d, annot in events.items():
    # イベント日に縦の点線を描く
    plt.axvline(x=d, color='black', linestyle='--', alpha=0.5)

    # プロットのy軸上の位置
    if d in df.index:
        y = df.loc[d]
    else:
        continue

    # プロット値の上部に注釈と矢印を表示
    plt.annotate(annot, xy=(d, y+mx*0.1), 
                 xytext=(d, y+mx*0.2), 
                 arrowprops=dict(color='red', headwidth=4, width=2, headlength=4),
                 horizontalalignment='left', verticalalignment='top')

plt.legend()
plt.show()

以下のようにグラフ上に矢印つきの注釈が表示された。

DTW（動的時間伸縮法）で時系列データ間の距離を求める

以前にPythonで時系列クラスタリングをするで、DTW（動的時間伸縮法）による時系列データのクラスタリングをやってみた。このときはPythonパッケージのtslearnを使ってクラスタリングをしたが、単にDTWによる時系列データ間の類似度（距離）を求めたいだけならfastdtwも使える。ただ、fastdtwではDTWの近似アルゴリズムを使用しているので、tslearnとfastdtwの結果を比べてみることにした。

環境

Debian Buster（Dockerコンテナ）とJupyter Lab。
# cat /etc/issue Debian GNU/Linux 10 \n \l # python -V Python 3.7.6 # pip list | grep -e pandas -e matplotlib -e scikit-learn -e jupyter jupyter 1.0.0 jupyter-client 5.3.4 jupyter-console 6.1.0 jupyter-core 4.6.1 jupyterlab 1.2.6 jupyterlab-server 1.0.6 matplotlib 3.1.3 pandas 1.0.1 scikit-learn 0.22.1

時系列データ

PythonでPDFの表からデータを抽出する（その２） で使用した帰国者・接触者相談センター(全相談件数)を使う。これは厚生労働省のホームページ地域ごとの感染状況等の公表についてからダウンロードできる「３．帰国者・接触者相談センターへの相談件数の推移（都道府県別・各日）・帰国者・接触者外来の受診者数の推移・うちPCR検査実施件数の推移（都道府県別・各日） 」のPDFに含まれるデータ。最近EXCELファイルが追加されたのでわざわざPDFからデータを抽出する必要はないのだが、ここではPythonでPDFの表からデータを抽出する（その２） の方法で作成したcsvを使う。都道府県ごとのデータがあるので、都道府県データ間の類似度をDTWで比べてみる。

PDFからcsvにしたデータをExcelで開くと次のような感じになる。

ライブラリなどのインストール

DTWを計算するためにtslearnとfastdtwをインストールする。
# pip install tslearn fastdtw

インストールしたライブラリのバージョン。
# pip list | grep -e tslearn -e fastdtw fastdtw 0.3.4 tslearn 0.4.1

さらに、matplotlibのグラフで日本語表示できるように日本語フォントをインストールしておく。
# apt install fonts-takao

データの準備

PDFから抽出したデータを保存したcsv（000643758.csv）を読み込んで、このcsvに含まれる都道府県ごとの「帰国者・接触者相談センター(全相談件数)」を使う。期間は５月のひと月分。都道府県により絶対数に差があるので、0-1の範囲にデータを変換する。なお、欠損値があるとDTWによる計算でエラーになるので京都、滋賀、香川のデータは除外。

%matplotlib inline

import pandas as pd
from sklearn import preprocessing

import matplotlib.pyplot as plt

# 日本語フォントの指定
plt.rcParams['font.family'] = 'TakaoGothic'

DT_FROM = '2020-05-01'
DT_TO = '2020-05-30'

def normalize(df):
    cols = df.columns
    idx = df.index
    x = df.values
    min_max_scaler = preprocessing.MinMaxScaler()
    x_scaled = min_max_scaler.fit_transform(x)
    return pd.DataFrame(x_scaled, index=idx, columns=cols)

def _load(csvpath):
    df = pd.read_csv(csvpath, index_col=0, header=[0,1], parse_dates=True, encoding='shift-jis')
 
    # 全国を除く「帰国者・接触者相談センター(全相談件数)」だけを選択
    df_con = df.filter(like='帰国者・接触者相談センター(全相談件数)', axis=1).drop(columns=('全国', '帰国者・接触者相談センター(全相談件数)'))
    
    # カラム名を都道府県のみにする
    df_con.columns = df_con.columns.droplevel(level=1)

    # 0-1の範囲に正規化    
    df_con = normalize(df_con[DT_FROM:DT_TO])

    # 欠損値(NaN)があるとDTWの計算ができないので除外
    # 京都、滋賀、香川のデータを除外
    df_con = df_con.dropna(axis='columns')    

    #print(df_con.info())
    #display(df_con.head())
    ax = df_con.plot(legend=None)
    plt.show()
    
    return df_con

df = _load('000643758.csv')

都道府県ごとのデータのプロット。

DTWの計算

tslearnとfastdtwで都道府県データ間の距離をDTWで計算し、似ている／似ていない都道府県の組み合わせを求めてみる。

%matplotlib inline

from itertools import combinations
import pandas as pd
from sklearn import preprocessing

from fastdtw import fastdtw
from tslearn.metrics import dtw

import matplotlib.pyplot as plt

# 日本語フォントの指定
plt.rcParams['font.family'] = 'TakaoGothic'

DT_FROM = '2020-05-01'
DT_TO = '2020-05-30'

def calc_dtw(l_sr):
    dict_dtw = {}
    for sr1, sr2 in combinations(l_sr, 2):
        distance, path = fastdtw(sr1.to_numpy(), sr2.to_numpy())
        dict_dtw[(sr1.name, sr2.name)] = (dtw(sr1.to_numpy(), sr2.to_numpy()), distance)
    
    return dict_dtw

df = _load('000643758.csv')

# Seriesのリスト
l_sr = [item for l, item in df.items()]

# 都道府県間のDTWを計算
dict_dtw = calc_dtw(l_sr)

print('\n似ている都道府県トップ10（DTW）')
display(sorted(dict_dtw.items(), key=lambda x:x[1][0])[:10])
print('\n似ている都道府県トップ10（fastDTW）')
display(sorted(dict_dtw.items(), key=lambda x:x[1][1])[:10])
print('\n似ていない都道府県トップ10（DTW）')
display(sorted(dict_dtw.items(), key=lambda x:x[1][0], reverse=True)[:10])
print('\n似ていない都道府県トップ10（fastDTW）')
display(sorted(dict_dtw.items(), key=lambda x:x[1][1], reverse=True)[:10])

似ている都道府県の組み合わせトップ10は以下の通り。tslearnとfastdtwで違いはあるものの、最も距離が近い組み合わせの神奈川と富山はtslearnとfastdtwで同じ。

似ていない都道府県の組み合わせトップ10は以下の通り。こちらもtslearnとfastdtwで違いはあるものの、最も距離が遠い組み合わせの岐阜と高知はtslearnとfastdtwで同じ。

最も距離が短い神奈川と富山のプロット。ほぼ一致している。

最も距離が遠い岐阜と高知のプロット。確かにけっこう違う。

tslearnとfastdtwそれぞれで距離が短いトップ5を比較。以下のコードを上記コードのあとに追加する。

（2021.12.16 グラフ作成部のコードを追加）

tslearn_top5 = sorted(dict_dtw.items(), key=lambda x:x[1][0])[:5]
fastdtw_top5 = sorted(dict_dtw.items(), key=lambda x:x[1][1])[:5]

for i, ((prefs1, values1), (prefs2, values2)) in enumerate(zip(tslearn_top5, fastdtw_top5)):
    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 3))
    df.loc[:, prefs1].plot(ax=axes[0])
    df.loc[:, prefs2].plot(ax=axes[1])
    plt.show()   

左がtslearnで右がfastdtw。どちらが良いのかは判断が難しい・・

ipywidgetsとJupyter Notebookで経済センサスのコロプレス図を作成する

PythonとPlotlyでインタラクティブなコロプレス図を作成するではPlotlyでインタラクティブなコロプレス図を作成したが、ipywidgetsを使うとJupyter Notebook上にインタラクティブなウィジェットを作成できる。今回はe-Statからダウンロードした経済センサスのデータを使って、ドロップダウンウィジェットで項目を選択して表示項目を変えられるコロプレス図を作成してみる。

環境

Windows10(1903)のWSL(Ubuntu 18.04)とJupyter Notebookを使用。
$ lsb_release -dr Description: Ubuntu 18.04.3 LTS Release: 18.04 $ python3 -V Python 3.6.9 $ pip3 show pandas Name: pandas Version: 0.25.3 ... $ jupyter notebook --version 6.0.3

必要なライブラリのインストール

GeoPandas、Plotly、ipywidgetsをインストールする。それぞれpipでインストールできる。
$ pip3 install geopandas plotly ipywidgets

バージョンはそれぞれ以下の通り。
$ pip3 list | grep -e geopandas -e plotly -e ipywidgets geopandas (0.6.2) ipywidgets (7.5.1) plotly (4.5.0)

さらに以下のコマンドも実行しておく。
$ jupyter nbextension enable --py widgetsnbextension

使用するデータ

e-Statにある平成28年経済センサス（活動調査）のデータを使う。e-Statにある平成28年経済センサス（活動調査）には分類や集計単位が異なるいくつかのデータがあるが、使用するのは以下のデータ。

• 統計名：経済センサス‐活動調査 平成28年経済センサス－活動調査 事業所に関する集計 産業別集計 サービス関連産業Ｂに関する集計
• 表番号：7
• 表題：サービス関連産業Ｂ（細分類），単独・本所・支所（３区分）別民営事業所数，従業者数，売上（収入）金額及び収入を得た相手先別収入額―全国，都道府県

平成28年経済センサス（活動調査）については平成28年経済センサス‐活動調査を参照。ラーメン店やホテルなど様々な産業の事業所数や収入などのデータがある。経済センサスの産業はたくさんあるが、平成28年経済センサス‐活動調査 産業分類一覧に一覧がある。

Pythonとe-StatのAPIで統計データを取得するのようにAPIでもデータ取得できるが、今回はcsvをダウンロードする。e-Statでは項目を選択したり、表示形式を変更したりしてダウンロードできるが、今回は変更はせずにcsvダウンロード時に「注釈を表示する」のチェックを外した以外はデフォルトのまま。データ量が多いので2ファイルになる。

ダウンロードしたcsvの先頭の数行には統計名などデータの情報があり、その後にデータが部が続く。データは都道府県などの「地域」や「H28_産業分類」ごとのデータで、「H28_単独・本所・支所」と「表章項目」の2種類のカラムがある。

2つのcsvを取り込んで、最終的に「地域」「H28_産業分類」「H28_単独・本所・支所」の3つをインデックスとするマルチインデックスのDataFrameにする。

import pandas as pd

def load_sensus7(path_to_csv):
    print('*** Loading {} ***'.format(path_to_csv))

    # 都道府県をインデックスとして読み込む
    # 最終header行より前の行はスキップされる
    df = pd.read_csv(path_to_csv, dtype=object, header=[10, 12], index_col=1, skipinitialspace=True, encoding='shift-jis')

    # 地域コード、時間軸コード、時間軸、表章項目カラムを削除
    df.drop(columns=df.columns[:4], inplace=True)
    df.drop(columns=df.columns[1], inplace=True)

    # 産業分類をインデックスに追加してマルチインデックスにする
    df.set_index(df.columns[0], append=True, inplace=True)
    df.index.set_names(['都道府県', '産業分類'], inplace=True)
    
    # set_indexを使うとマルチカラムインデックスがシングルインデックスカラムになってしまうので再設定
    df.columns = pd.MultiIndex.from_tuples(df.columns, names=['単独・本所・支所', '項目'])
    
    # 「単独・本所・支所」カラムをマルチインデックスに追加
    # stackを実行すると順番が変わるのでreindexで元の順番を維持させる
    column_l = [lv1 for lv0, lv1 in df.columns]
    column_sorted = sorted(set(column_l), key=column_l.index)
    index_sorted = []
    base_l = [base for base, item in df.columns]
    base_sorted = sorted(set(base_l), key=base_l.index)
    for pref, industry in df.index:
        for base in base_sorted:
            index_sorted.append((pref, industry, base))
            
    df = df.stack(level=0).reindex(index_sorted).reindex(columns=column_sorted)


    # 取り込んだデータのカンマを除去して数値に変換
    df = df.applymap(lambda x: pd.to_numeric(x.replace(',', ''), errors='coerce'))

    # NaNを0とする
    # to_numericだとfloatになってしまうのでここでintにする
    df = df.fillna(0).astype('int64')

    # 「全国」の行を削除
    return df.iloc[~(df.index.get_level_values(0) == '全国')]

def main(path_to_csv_l):
    df = pd.DataFrame()
    for path_to_csv in path_to_csv_l:
        df = df.append(load_sensus7(path_to_csv))

    display(df.head())

if __name__ == '__main__':
    # ダウンロードしたcsvを指定
    main(['FEH_00200553_200243012738.csv', 'FEH_00200553_200243012807.csv'])

以下のようなデータになる。

さらに、PythonとPlotlyでインタラクティブなコロプレス図を作成するで作成した、簡略化したGeoJSON形式の都道府県境界図（jpnmap_simplified.geojson）も使う。

ウィジェットとコロプレス図の作成

コロプレス図の作成はPythonとPlotlyでインタラクティブなコロプレス図を作成すると同様だが、さらにウィジェットを使って表示項目などをインタラクティブに変更できるようにする。

ドロップボックスウィジェットとしてコロプレス図背景地図の「背景地図」、事業所数や従業員数などの「項目」、ラーメン店などの「産業分類」、事業所の本社・支店数などの「単独・本所・支所」の4つを作成する。なお、コロプレス図に使用するデータはすべて人口1000人あたりの数値とする。コードは以下の通り（経済センサスのcsvを読み込む関数は省略）。

import json

import pandas as pd
import geopandas as gpd
import plotly.express as px
import ipywidgets as widgets

# 都道府県境界図のパス
OUTLINEMAP_PATH = 'jpnmap_simplified.geojson'

# コロプレス図の中心
CENTER = {'lat': 36, 'lon': 140}

def load_outlinemap():
    outline = gpd.read_file(OUTLINEMAP_PATH, driver='GeoJSON')

    # GeoJSON形式で保存するとGeoDataFrameのindexが保存されないので再設定する
    # to_jsonのときにindexの値がIDとなるため
    outline.set_index('ID', drop=False, inplace=True)
    print(outline.info())
    display(outline.head())
    
    return outline

def main(path_to_csv_l):

    def plot_choropleth(background, colname, **indices):
        
        # 「産業」「単独・本所・支所」で対象データをしぼる
        prefs_df = df.copy()
        for idx_name, idx_value in indices.items():
            prefs_df = prefs_df.iloc[prefs_df.index.get_level_values(idx_name)==idx_value]           
        
        # マルチインデックスの都道府県以外のindexを削除する（outlineとの計算でエラーになるため）
        prefs_df.index = prefs_df.index.droplevel(prefs_df.index.names[1:])

        # outlineとprefs_dfのindexをIDに合わせる
        prefs_df.set_index('ID', drop=False, inplace=True)
        
        # 人口1000人当たりの数値にする
        prefs_df['color'] = prefs_df[colname] * 1000 / outline['JINKO']
    
        # color: 表示対象の項目
        # locations: geojsonのIDに対応する項目
        # geojson: geojsonをdictに変換したもの。IDがふられている
        fig = px.choropleth_mapbox(prefs_df, geojson=json.loads(outline['geometry'].to_json()), locations='ID', color='color',
                                title = '',
                                hover_name = 'PREF_NAME',
                                color_continuous_scale='OrRd', 
                                mapbox_style=background,
                                zoom=3.5, 
                                center = CENTER,
                                opacity=0.5,
                                labels={colname:colname}
                                )
        fig.update_layout(margin={'r':10,'t':30,'l':10,'b':10, 'pad':5})
        fig.show()

    
    # 都道府県境界図の読み込み
    outline = load_outlinemap()

    df = pd.DataFrame()
    for path_to_csv in path_to_csv_l:
        df = df.append(load_sensus7(path_to_csv))
    
    # 経済センサスのデータと都道府県境界図のデータを、都道府県名をキーにしてマージ
    # dfにgeometryに対応するIDを設定
    idx = pd.MultiIndex.from_tuples(df.index.values, names=df.index.names)
    df = df.merge(outline[['PREF_NAME', 'ID']], left_on=df.index.names[0], right_on='PREF_NAME', how='left')
    df.index = idx
    
    print(df.info())
    display(df.head())

    # 背景地図のウィジェット
    background = widgets.Dropdown(
        options=['white-bg', 'carto-positron', 'carto-darkmatter', 'stamen-terrain'],
        value='white-bg',
        description='背景地図:',
    )

    # 項目用ウィジット
    # PREF_NAMEとIDを対象外にする
    item = widgets.Dropdown(
        options=df.columns[:-2],
        value=df.columns[0],
        description='項目:',
    )

    # index項目用Dropdown作成
    dropdown_indices = {}
    level_n = df.index.nlevels
    for level in range(1, level_n):
        index_l = list(df.index.get_level_values(level).values)
        index_l_u = sorted(set(index_l), key=index_l.index)
        dropdown = widgets.Dropdown(
            options=index_l_u,
            value=index_l_u[0],
            description=df.index.names[level] + ':',
        )
        dropdown_indices[df.index.names[level]] = dropdown

    widgets.interact(plot_choropleth, background=background, colname=item, **dropdown_indices)
        
if __name__ == '__main__':
    # ダウンロードしたcsvを指定
    main(['FEH_00200553_200243012738.csv', 'FEH_00200553_200243012807.csv'])

次のようなコロプレス図とその左上にドロップボックスのウィジェットが作成された。

Jupyter Notebookで表示範囲が狭く見にくい場合は、メニューでCell＞All Output＞Toggle Scrollingとすると表示範囲が広くなる。

背景地図と産業分類を変更してみる。

全体のコードはこちら。

PythonとPlotlyでインタラクティブなコロプレス図を作成する

Pythonでシェープファイルから都道府県境界図を作成するではGeoPandasで都道府県境界図のシェープファイルを作成し、さらにコロプレス図をプロットしてみた。今回はPlotlyを利用してインタラクティブなコロプレス図を作成してみる。Plotlyを使うとPlotlyでStack OverflowのDeveloper Surveyを可視化するのように、マウスオーバーでテキストを表示したり、拡大・縮小などを行える図を作成できる。

環境

Windows10(1903)のWSL(Ubuntu 18.04)とJupyter Notebookを使用。
$ lsb_release -dr Description: Ubuntu 18.04.3 LTS Release: 18.04 $ python3 -V Python 3.6.9 $ jupyter notebook --version 6.0.3

必要なライブラリのインストール

GeoPandasとPlotlyをインストールする。両方ともpipでインストールできる。
$ pip3 install geopandas plotly

バージョンはそれぞれ以下の通り。
$ pip3 list | grep -e geopandas -e plotly geopandas (0.6.2) plotly (4.5.0)

GeoJSON形式の境界図の作成

Plotlyでコロプレス図を作成するにはGeoJSON形式のデータが必要になるので、シェープファイルではなくGeoJSON形式の都道府県境界図を用意する。GeoPandasではGeoJSON形式でも保存できる。Pythonでシェープファイルから都道府県境界図を作成するで保存した都道府県境界図のシェープファイルjpnmap.shpを読み込んで、GeoJSON形式で保存し直す。

import geopandas as gpd

jpn = gpd.read_file('jpnmap.shp', encoding = 'shift-jis')

# コロプレス図で使用する人口密度の計算
# AREA（面積）の単位は平方メートルなので平方キロメートルに直して計算
jpn['DENSITY'] = jpn['JINKO'] / (jpn['AREA'] / (10**6))

# コロプレス図を作成するときにgeometryとの対応をとるためのID
jpn['ID'] = jpn.index.values

# GeoJSON形式で保存
jpn.to_file('jpnmap.geojson', driver='GeoJSON')

GeoJSON形式ファイルの軽量化

GeoJSON形式の境界図ができたのでコロプレス図の作成に取りかかりたいところだが、ここで問題がある。作成したGeoJSON形式のファイルは200MBほどあり、これを読み込んでJupyter Notebookでコロプレス図を表示しようとしたらブラウザが落ちてしまった。データサイズが大きすぎてメモリ消費が増えすぎたことが原因のようなので、GeoJSON形式境界図の軽量化を行う。

軽量化にはTopoJSONを使用する。TopoJSONはGeoJSONの拡張形式で、GeoJSONからTopoJSONへの変換などいくつかのコマンドラインツールが用意されている。それらを使ってジオメトリのデータを間引く簡略化を行い、ファイルを軽量化する。

まずはTopoJSONのインストール。TopoJSONのインストールにはnpmが必要なので先にインストールする。
$ sudo apt install npm

続いてTopoJSONをインストール。
$ sudo npm install -g topojson

以下のコマンドでGeoJSON形式の境界図データをTopoJSON形式に変換して簡略化する。geo2topoでTopoJSON形式に変換、toposimplifyで簡略化。toposimplifyのオプションPについてはtoposimplifyにドキュメントがあるが、0から1の範囲の値を指定する。この値が小さいほど簡略化の程度が大きくなる。オプションfを指定すると、簡略化の閾値よりも小さな分離されたリング（ポリゴンの境界線）を取り除く。jpnはオブジェクト名で、任意の名前でOK。
$ geo2topo jpn=jpnmap.geojson | toposimplify -P 0.05 -f > simplified.topojson

簡略化したらGeoJSON形式へ再変換する。
$ topo2geo jpn=jpnmap_simplified.geojson < simplified.topojson

結果、200MBを超えていたファイルサイズが約10MBになった。

Plotlyでインタラクティブなコロプレス図を作成

PlotlyにExpressというAPIがあり、そのAPIのchoropleth_mapboxというメソッドを使うとGeoJSON形式のデータを読み込んでコロプレス図を作成できる。以下は簡略化したGeoJSON形式の境界図jpnmap_simplified.geojsonを使ってコロプレス図を作成するPythonコード。

import json
import geopandas as gpd
import plotly.express as px

# コロプレス図の中心
CENTER = {'lat': 36, 'lon': 140}

jpn = gpd.read_file('jpnmap_simplified.geojson')

# GeoJSON形式で保存したときにGeoDataFrameのindexが保存されず0からのindexがふられるのでindexを再設定する
# 今回はたまたまもとのindexが0からはじまる数字なので再設定しなくても動作する
jpn.set_index('ID', drop=False, inplace=True)
 
print(jpn.info())
display(jpn.head())

# color: 表示対象の項目
# locations: GeoJSONのIDに対応する項目
# geojson: GeoJSONをdictに変換したもの。このデータではID(0-46)がふられている
fig = px.choropleth_mapbox(jpn, geojson=json.loads(jpn['geometry'].to_json()), locations='ID', color='DENSITY',
      title = '人口密度',
      hover_name = 'PREF_NAME',
      color_continuous_scale="OrRd", 
      mapbox_style='white-bg',
      zoom=4, 
      center=CENTER,
      opacity=0.5,
      labels={'DENSITY':'人口密度'}
      )
fig.update_layout(margin={'r':10,'t':40,'l':10,'b':10,'pad':5})
fig.show()

次のようなコロプレス図が作成された。マウスオーバーするとchoropleth_mapboxのhover_name、locations、labelsで指定した項目が表示される。

PythonでRのスクリプトを実行する

Raspberry Piに最新のR(3.5.2)をインストールするでRaspberry PiにRをインストールしたので、PythonからRを使ってみる。

PythonからRを使う方法としてはPypeRなどがあるらしいが、Rで検定を実行してその結果を表示するだけなら、Pythonのsubprocessモジュールでも十分。Pythonのsubprocessでコマンドを実行するのようにPythonからコマンドを実行できるので、RscriptをPythonから実行する。

環境

Raspberry PiとRaspbian Stretch with Desktop。
$ lsb_release -a No LSB modules are available. Distributor ID: Raspbian Description: Raspbian GNU/Linux 9.6 (stretch) Release: 9.6 Codename: stretch $ python3 -V Python 3.5.3 $ R --version R version 3.5.2 (2018-12-20) -- "Eggshell Igloo" Copyright (C) 2018 The R Foundation for Statistical Computing Platform: armv7l-unknown-linux-gnueabihf (32-bit)

RスクリプトをPythonで実行する

Rスクリプトをコマンドラインで実行するにはRscriptコマンドを使う。このコマンドは引数を受け取れるので、PythonでRスクリプトファイルを作成して、引数を含めてsubprocessモジュールで実行する。

以下のコードでは、平均が異なる正規分布に従う値を10ずつ生成してt検定を行うRスクリプトファイルを作成し、Rscriptで実行した結果を表示する。

※Rスクリプトがエラーのときにメッセージを出力するようにコード変更（2019.2.11）

結果は以下の通り。
g1= [1] 0.5362420 0.9076155 0.8669404 0.7854252 -0.0402832 0.7736883 [7] 0.1422077 -0.4394270 -1.5471381 0.2457207 Min. 1st Qu. Median Mean 3rd Qu. Max. -1.54714 0.00534 0.39098 0.22310 0.78249 0.90762 g2= [1] 1.178631 1.688101 2.748449 3.125478 1.366282 1.800813 2.291191 2.578756 [9] 2.612918 2.660056 Min. 1st Qu. Median Mean 3rd Qu. Max. 1.179 1.716 2.435 2.205 2.648 3.125 Two Sample t-test data: g1 and g2 t = -6.2329, df = 18, p-value = 7.019e-06 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -2.650029 -1.313908 sample estimates: mean of x mean of y 0.2230992 2.2050676

関連リンク

• Raspberry Piに最新のR(3.5.2)をインストールする

階層的クラスタリングでユークリッド距離とコサイン類似度の結果を比べてみる

階層的クラスタリングでは距離や類似度をもとにクラスターを作っていくわけだけど、距離や類似度にはいくつか種類があって、分析に適したものを選択することになる。距離としてはユークリッド距離が使われる例がたくさんある一方で、テキスト分析ではコサイン類似度がよく使われるらしい。じゃぁ実際のところ、この2つはどう違うのだろうかというのが疑問で、実際に同じデータをPythonのScipyで階層的クラスタリングして比べてみた。

環境

Bash on Ubuntu on WindowsとJupyter Notebook。
$ cat /etc/issue Ubuntu 16.04.3 LTS \n \l $ python3 -V Python 3.5.2 $ pip3 show scikit-learn pandas matplotlib scipy Name: scikit-learn Version: 0.19.1 ... Name: pandas Version: 0.22.0 ... Name: matplotlib Version: 2.1.0 ... Name: scipy Version: 1.0.0 ... $ jupyter --version 4.4.0

ユークリッド距離とコサイン類似度

ユークリッド距離は、一般的に言うところの2点間の距離なので感覚的に理解できる。今回はScipyで階層的クラスタリングをするけれど、デフォルトはユークリッド距離になっている。一方のコサイン類似度は２つのベクトルによる角度をθとしたときのcosθで、θが0度のときは1，90度のときは0になる。つまり角度が小さいほど類似度が1に近づいて高くなる。コサイン類似度（高校数学の美しい物語）の説明がわかりやすい。

データの準備

階層的クラスタリングするためのデータを準備する。データはscikit-learnのirisデータセットを使う。まずはデータセットをPandasのDataframeにして正規化する。

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.preprocessing import MinMaxScaler

%matplotlib inline
import matplotlib.pyplot as plt

def iris_df():
    # scikit-learnのirisデータセット読み込み
    iris = load_iris()
    # irisデータセットをPandasのDetaFrameに変換
    df = pd.DataFrame(data= np.c_[iris['data'], iris['target']], columns= iris['feature_names'] + ['target'])

    # 正規化
    items = df.columns
    index = df.index
    df = pd.DataFrame(MinMaxScaler().fit_transform(df))
    df.columns = items #カラム名を再設定
    df.index = index # インデックス名を再設定
    
    return df

if __name__ == '__main__':
    df = iris_df()
    print(df)

以下のようなDataFrameになる。

Pythonで階層的クラスタリング

Scipyのlinkageで階層的クラスタリングを行い、dendrogramでデンドログラム（樹形図）を作成する。まずはユークリッド距離から。比較しやすいように使うデータはpetal length (cm)とpetal width (cm) の2次元のみにする。クラスタリング手法はウォード法。

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.preprocessing import MinMaxScaler
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage

%matplotlib inline
import matplotlib.pyplot as plt

def draw_dendrogram(Z, labels):
    # デンドログラムの作成
    plt.figure(figsize=(14, 5))
    plt.ylabel('Distance')
    dendrogram(
        Z,
        leaf_rotation=90.,
        leaf_font_size=8.,
        labels=labels
    )
    plt.show()

if __name__ == '__main__':
    df = iris_df()
    print(df)

    # petal length (cm)とpetal width (cm) で階層的クラスタリング
    Z = linkage(pdist(df[df.columns[2:4]].as_matrix(), metric='euclidean'), method='ward')
    draw_dendrogram(Z, df.index.values)

次はユークリッド距離の代わりにコサイン類似度で階層的クラスタリングをする。

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.preprocessing import MinMaxScaler
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage

%matplotlib inline
import matplotlib.pyplot as plt

def draw_dendrogram(Z, labels):
    # デンドログラムの作成
    plt.figure(figsize=(14, 5))
    plt.ylabel('Distance')
    dendrogram(
        Z,
        leaf_rotation=90.,
        leaf_font_size=8.,
        labels=labels
    )
    plt.show()

if __name__ == '__main__':
    df = iris_df()
    print(df)

    # petal length (cm)とpetal width (cm) で階層的クラスタリング
    Z = linkage(pdist(df[df.columns[2:4]].as_matrix(), metric='cosine'), method='ward')

    draw_dendrogram(Z, df.index.values)

階層的クラスタリングの結果を比較してみる

デンドログラムだとわかりずらいので散布図で比較してみる。

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.preprocessing import MinMaxScaler
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

%matplotlib inline
import matplotlib.pyplot as plt

def draw_scatter(df):
    # targetごとに色分けしてプロットする
    groups = df.groupby('fcluster')
    for status, group in groups:
        plt.plot(group[df.columns[2]], group[df.columns[3]], marker='o', linestyle='', ms=4)
    plt.grid(True)
    plt.show()

if __name__ == '__main__':
    df = iris_df()
    print(df)

    metric_t = ('euclidean', 'cosine') # ユークリッド距離、コサイン類似度
    th_t = (4.0, 0.6) # 分割する距離の閾値
    for metric, th in zip(metric_t, th_t):
        # petal length (cm)とpetal width (cm) で階層的クラスタリング
        Z = linkage(pdist(df[df.columns[2:4]].as_matrix(), metric=metric), method='ward')
        draw_dendrogram(Z, df.index.values)

        # 距離をもとにフラットクラスターを作成
        df['fcluster'] = fcluster(Z, th, criterion='distance')

        draw_scatter(df)

ユークリッド距離を使ったときの結果。

コサイン類似度を使ったときの結果。

コサイン類似度のときはベクトルの角度で分割されているのがわかる。

全体コードはこちら。

Pythonで分類木を視覚化する

都道府県別の飲み物の購入額をk-means法でクラスタリングするでk-means法によるクラスタリングで分類したが、どういう基準で分類されたのかわからない。そこでこのクラスタリングの結果をもとに決定木の一種である分類木を作成して、分類基準やその過程を確認してみる。scikit-learnで分類木による分類を行い、Graphvizで結果を視覚化する。

環境

Raspbian JessieとJupyter Notebook。
$ cat /etc/os-release PRETTY_NAME="Raspbian GNU/Linux 8 (jessie)" NAME="Raspbian GNU/Linux" VERSION_ID="8" VERSION="8 (jessie)" ... $ python3 -V Python 3.4.2 $ pip3 show scikit-learn Name: scikit-learn Version: 0.19.0 ... $ pip3 show pandas Name: pandas Version: 0.20.2 ... $ pip3 show jupyter Name: jupyter Version: 1.0.0 ...

1．データの取得

使用するデータは都道府県別の飲み物の購入額をk-means法でクラスタリングするで使用したe-Statの2009年の全国消費実態調査データの都道府県ごとの「緑茶」「紅茶」「コーヒー」の消費額と、k-means法によるクラスタリングの結果。e-StatはAPIでデータ取得できるので、Pythonとe-StatのAPIで統計データを取得すると同じコード（estat.py）で2009年の全国消費実態調査データの都道府県ごとの「緑茶」「紅茶」「コーヒー」の消費額を取得し、k-means法によるクラスタリングを行う。

import pandas as pd
from sklearn.cluster import KMeans

import estat

def calc_rate(df):
    # データが文字列なので数値に変換
    df = df.apply(pd.to_numeric, errors='coerce')

    df['sum'] = df[df.columns[0]] + df[df.columns[1]] + df[df.columns[2]]
    df['緑茶R'] = df[df.columns[0]] / df['sum']
    df['紅茶R'] = df[df.columns[1]] / df['sum']
    df['コーヒーR'] = df[df.columns[2]] / df['sum']
 
    # 「緑茶」「紅茶」「コーヒー」の、都道府県ごとの購入金額の割合
    return df[df.columns[4:]]

if __name__ == '__main__':
    # e-Stat APIでデータ取得
    res = estat.get_estat_json()
    df = estat.to_dataframe(res)
 
    df = calc_rate(df)
    # 5つにクラスタリング
    df['cluster'] = KMeans(n_clusters=5, random_state=0).fit_predict(df.as_matrix())
    print(df)

結果は以下のようになる。

このデータを使って分類木を作成する。

2．Pythonで分類木を可視化するための準備

データの分類はscikit-learnでできるが、図の作成にはGraphvizが必要。pipでGraphvizバイナリとPythonパッケージがインストールできる。
$ sudo pip3 install graphviz

3. 分類木による分類と可視化

分類木による分類はsklearn.tree.DecisionTreeClassifierクラスでできる。さらに分類結果をJupyter Notebookで可視化するコードは以下の通り。

from IPython.display import Image
from IPython.display import display

import pandas as pd
from sklearn.cluster import KMeans
from sklearn import tree
import graphviz 

import estat

def classification_tree(df):
    X = df[df.columns[:-1]]
    y = df[df.columns[-1]]

    # 最大深さ3の分類木を作成
    clf = tree.DecisionTreeClassifier(max_depth=3)    
    clf.fit(X, y)

    """
    分類木の可視化
    """
    # 分類木をDOT言語で出力
    dot_data = tree.export_graphviz(clf, out_file = None, feature_names = df.columns[:-1],
                           class_names = df.columns[-1], filled = True, rounded = True,
                           special_characters=True)
    graph = graphviz.Source(dot_data) 
    graph.format = 'png' # pngで保存
    display(Image(graph.render('ct_drink'))) # 保存したpngを表示
    
if __name__ == '__main__':
    # e-Stat APIでデータ取得
    res = estat.get_estat_json()
    df = estat.to_dataframe(res)
 
    df = calc_rate(df)
    # 5つにクラスタリング
    df['cluster'] = KMeans(n_clusters=5, random_state=0).fit_predict(df.as_matrix())
    print(df)

    classification_tree(df)

4．結果

分類木を可視化した結果は以下の通り。こうしてみると紅茶は緑茶、コーヒーに比べて割合が少なかったせいか、分類の判断に使われていない。

5．標準化してクラスタリングした結果で分類木を作成

クラスタリングする際にデータを標準化することで紅茶が分類基準に使われるようになるか試してみる。

以下のようにsklearn.preprocessing.StandardScalerクラスで標準化する。

    df = calc_rate(df)

上記コード部に以下のように標準化のコードを追加。

from sklearn.preprocessing import StandardScaler
...

    df = calc_rate(df)

    items = df.columns
    index = df.index
    df = pd.DataFrame(StandardScaler().fit_transform(df))
    df.columns = items #カラム名を再設定
    df.index = index # インデックス名を再設定

結果は以下の通り。クラスタリングの結果も変わっているが、今度はコーヒーが分類に使われていない。緑茶や紅茶の差の影響がより大きいということなのかもしれない。

標準化すると値がマイナスになって正しくクラスタリングされないと思われるので正規化に変更。5．を差し替え。（2018/7/31）

5．正規化してクラスタリングした結果で分類木を作成

クラスタリングする際にデータを正規化することで紅茶が分類基準に使われるようになるか試してみる。

以下のようにsklearn.preprocessing.MinMaxScalerクラスで正規化する。

    df = calc_rate(df)

上記コード部に以下のように正規化のコードを追加。

from sklearn.preprocessing import MinMaxScaler
...

    df = calc_rate(df)

    items = df.columns
    index = df.index
    df = pd.DataFrame(MinMaxScaler().fit_transform(df))
    df.columns = items #カラム名を再設定
    df.index = index # インデックス名を再設定

結果は以下の通り。クラスタリングの結果も変わっているが、今度は紅茶も分類に使われるようになった。

Pythonとe-StatのAPIで統計データを取得する

都道府県別の飲み物の購入額をk-means法でクラスタリングするではe-Statからcsvをダウンロードして利用したが、e-StatにはAPIがあって、さまざまな統計データをプログラムで取得できる。Python3でe-StatのAPIでデータを取得して、最終的にPandasのDataFrameに変換する手順をまとめておく。

環境

Raspbian JessieとPython3.4.2。
$ cat /etc/os-release PRETTY_NAME="Raspbian GNU/Linux 8 (jessie)" NAME="Raspbian GNU/Linux" VERSION_ID="8" VERSION="8 (jessie)" ... $ python3 -V Python 3.4.2 $ pip3 show pandas Name: pandas Version: 0.20.2 ...

1．アプリケーションIDの取得

APIを使用するにはアプリケーションIDが必要で、それを取得するにはまずe-Statの「新規登録」でユーザー登録をする。ユーザー登録が完了したらログインして、「マイページ」ページの「API機能（アプリケーションID発行）」へ移動する。名称やURLを入力する必要があるが、利用ガイドにあるように、公開サイトで利用するのでなければURLは「http://localhost」でOK。名称も後で変更できるので、ちょっと試すだけなら、とりあえず深く考えずにつけても問題なさそう。

「発行」ボタンを押すとappIdの欄にアプリケーションIDが表示されるので、これをどこかにコピーしておく。

2．リクエストURLの取得

APIを使用するためのリクエストURLを取得する。URLは各データのページで取得するデータを表示したうえで「API」ボタンを押すことで得られる。

以下のようなポップアップウィンドウにAPI用のURLが表示される。

今回は都道府県別の飲み物の購入額をk-means法でクラスタリングすると同じデータを取得する。2009地域は全国と各都道府県、2009品目分類は「緑茶」「紅茶」「コーヒー」のみで、あとは合計値が選択されるようにしておく。以下は選択した状態の画面。

「表示を更新」を押すと選択したデータが表示される。この画面で「API」ボタンを押して表示されるURLをコピーしておく。

3．APIによるデータの取得

APIで取得できるデータ形式はXML、JSON、JSONP、CSVの4種類。取得するデータ形式に合わせてコピーしておいたAPIのURLを変更し、さらにアプリケーションIDを設定してURLを使う。API仕様の詳細はAPI仕様を参照。

JSONで取得するPython3のコードは以下の通り。

import urllib.request
import json

# e-Stat APIのappId
appId = 'コピーしておいたappIdをここにペーストする'
# APIのリクエストURL
# コピーしておいたURLのhttpをhttpsにしておく
url = 'https://api.e-stat.go.jp/rest/2.1/app/getStatsData?cdCat03=47001&cdCat04=48001&cdArea=00000%2C01000%2C02000%2C03000%2C04000%2C05000%2C06000%2C07000%2C08000%2C09000%2C10000%2C11000%2C12000%2C13000%2C14000%2C15000%2C16000%2C17000%2C18000%2C19000%2C20000%2C21000%2C22000%2C23000%2C24000%2C25000%2C26000%2C27000%2C28000%2C29000%2C30000%2C31000%2C32000%2C33000%2C34000%2C35000%2C36000%2C37000%2C38000%2C39000%2C40000%2C41000%2C42000%2C43000%2C44000%2C45000%2C46000%2C47000&cdCat02=H1000%2CH1010%2CH1050&cdCat01=07001&appId=&lang=J&statsDataId=0003018185&metaGetFlg=Y&cntGetFlg=N&sectionHeaderFlg=1'

def url_json(url):
    """
    URLをJSON取得用に変更する
    「/apps/」のあとに「json/」を追加
    """
    return url.replace('/app/', '/app/json/')

def set_appid(url):
    """
    URLにappIdを設定する
    「appId=」のあとにappIdを追加
    """
    return url.replace('appId=', 'appId=' + appId)

def get_estat_json():
    request_url = set_appid(url_json(url))

    with urllib.request.urlopen(request_url) as f:
        res = json.loads(f.read().decode('utf8'))

    return res

def main():
    res = get_estat_json()
    # 取得したデータ（JSON）を表示
    print(json.dumps(res, indent=2, sort_keys=True))

if __name__ == '__main__':
    main()

実行すると以下のようにJSONでデータが得られる。
{ "GET_STATS_DATA": { "PARAMETER": { "CNT_GET_FLG": "N", "DATA_FORMAT": "J", "LANG": "J", "METAGET_FLG": "Y", "NARROWING_COND": { "CODE_AREA_SELECT": "00000,01000,02000,03000,04000,05000,06000,07000,08000,09000,10000,11000,12000,13000,14000,15000,16000,17000,18000,19000,20000,21000,22000,23000,24000,25000,26000,27000,28000,29000,30000,31000,32000,33000,34000,35000,36000,37000,38000,39000,40000,41000,42000,43000,44000,45000,46000,47000", "CODE_CAT01_SELECT": "07001", "CODE_CAT02_SELECT": "H1000,H1010,H1050", "CODE_CAT03_SELECT": 47001, "CODE_CAT04_SELECT": 48001 }, "SECTION_HEADER_FLG": 1, "START_POSITION": 1, "STATS_DATA_ID": "0003018185" }, ...

4．JSONをPandasのDataFrameに変換する

後の処理をやりやすくするためにJSONをPandasのDataFrameに変換する。まずはJSONをpandas.io.json.json_normalizeでDataFrameに変換し、pandas.DataFrame.pivotで都道府県をインデックス、品目をカラムにする。

import urllib.request
import json
import pandas as pd
from pandas.io.json import json_normalize

# e-Stat APIのappId
appId = 'コピーしておいたappIdをここにペーストする'
# APIのリクエストURL
# コピーしておいたURLのhttpをhttpsにしておく
url = 'https://api.e-stat.go.jp/rest/2.1/app/getStatsData?cdCat03=47001&cdCat04=48001&cdArea=00000%2C01000%2C02000%2C03000%2C04000%2C05000%2C06000%2C07000%2C08000%2C09000%2C10000%2C11000%2C12000%2C13000%2C14000%2C15000%2C16000%2C17000%2C18000%2C19000%2C20000%2C21000%2C22000%2C23000%2C24000%2C25000%2C26000%2C27000%2C28000%2C29000%2C30000%2C31000%2C32000%2C33000%2C34000%2C35000%2C36000%2C37000%2C38000%2C39000%2C40000%2C41000%2C42000%2C43000%2C44000%2C45000%2C46000%2C47000&cdCat02=H1000%2CH1010%2CH1050&cdCat01=07001&appId=&lang=J&statsDataId=0003018185&metaGetFlg=Y&cntGetFlg=N&sectionHeaderFlg=1'

def to_dataframe(res):
    item = [] # 品目のリスト
    pref = [] # 都道府県のリスト
    for obj in res['GET_STATS_DATA']['STATISTICAL_DATA']['CLASS_INF']['CLASS_OBJ']:
        # 分類などの出力
        if isinstance(obj['CLASS'], dict): # dict
            print(obj['@id'], obj['@name'], obj['CLASS']['@name'])
        else: # list
            id = obj['@id']
            for cat in obj['CLASS']:
                print(id, obj['@name'], cat['@name'], cat['@code'])
                if id == 'cat02': # 2009品目分類
                    item.append(cat['@name'])
                elif id == 'area': # 2009地域
                    pref.append(cat['@name'])

    # JSONをDataFrameに変換
    df = json_normalize(res['GET_STATS_DATA']['STATISTICAL_DATA']['DATA_INF']['VALUE'])
    # 必要なカラムのみ抽出
    df = df.loc[:, ['$', '@area', '@cat02']]

    # 都道府県をインデックス、品目（緑茶、紅茶、コーヒー）をカラムにする
    g = df.pivot(index='@area', columns='@cat02')
    g.index = pref
    g.columns = item
    g = g.iloc[1:48] # 「全国」の行を省く

    return g

def main():
    res = get_estat_json()
    # 取得したデータ（JSON）を表示
    print(json.dumps(res, indent=2, sort_keys=True))

    df = to_dataframe(res)
    print(df)

if __name__ == '__main__':
    main()

５．結果

4のスクリプトを実行すると以下の結果が得られる。
緑茶 紅茶 コーヒー 北海道 373 66 473 青森県 232 33 402 岩手県 422 48 399 宮城県 370 38 419 秋田県 412 62 427 山形県 339 37 317 福島県 309 49 329 茨城県 419 60 383 栃木県 429 47 431 ...

都道府県別の飲み物の購入額をk-means法でクラスタリングするでダウンロードしたcsvからDataFrameに変換したデータと比べると、同じ数値になっている。
  緑茶【円】 紅茶【円】 コーヒー【円】 2009地域 北海道 373 66 473 青森県 232 33 402 岩手県 422 48 399 宮城県 370 38 419 秋田県 412 62 427 山形県 339 37 317 福島県 309 49 329 茨城県 419 60 383 栃木県 429 47 431 ...

k-means法の最適なクラスター数を選択する

都道府県別の飲み物の購入額をk-means法でクラスタリングするでは、Python3とscikit-learnを使ってk-means法によるクラスタリングを行ったが、このときにクラスター数は思いつきの5にした。この5という数字には何の根拠もないので、実際はベターなクラスター数がありうる。k-means法のクラスター数を決める方法としてシルエット分析とエルボー法があるので、この2つの方法で最適なクラスター数を求めてみる。参照は以下のサイト。

• k-meansの最適なクラスター数を調べる方法
• Scikit Learn - K-Means - Elbow - criterion
• Selecting the number of clusters with silhouette analysis on KMeans clustering

環境とクラスタリングするデータ

都道府県別の飲み物の購入額をk-means法でクラスタリングすると同じで、以下の環境とe-Statの全国消費実態調査 平成21年全国消費実態調査 全国 購入先・購入地域編の［二人以上の世帯］　【品目別1世帯当たり1か月間の支出】　地域，購入地域，購入先別データセット。

Bash on Ubuntu on Windows
$ cat /etc/issue Ubuntu 16.04.3 LTS \n \l $ python3 -V Python 3.5.2 $ pip3 show scikit-learn pandas matplotlib Name: scikit-learn Version: 0.19.1 ... Name: pandas Version: 0.22.0 ... Name: matplotlib Version: 2.1.0 ...

Jupyter Notebook
$ jupyter --version 4.4.0

シルエット分析

シルエット分析では各点（ここでは都道府県）のシルエット値を計算する。シルエット値とは、MathWorksのサイトシルエット プロットによると「他のクラスターの点と比べて、その点が自身のクラスター内の他の点にどれくらい相似しているかを示す尺度」で、-1から1の範囲の値。シルエット値が1に近く、かつシルエット値をプロットしたシルエット図でクラスター間の幅の差が最も少ないクラスター数が最適となる。

エルボー法

クラスター数を変えてクラスタリングしたときの各SSE(クラスター内誤差の平方和)をプロットしたエルボー図で、ひじのように曲がっているところのクラスター数が最適というもの。scikit-learnのKMeansクラスではfitメソッドを適用したあとに得られるinertia_ (Sum of squared distances of samples to their closest cluster center)がSSEに当たるらしい。

シルエット分析とエルボー法のPython3コード

データをPnadasのDataFframeにするまでは都道府県別の飲み物の購入額をk-means法でクラスタリングするとコードは同じ。今回はクラスター数2～8の中からシルエット分析とエルボー法で最適なクラスター数を選ぶ。シルエット分析のコードはSelecting the number of clusters with silhouette analysis on KMeans clusteringのコードをほぼそのまま使えるが、長くなるのでクラスター数ごとの平均シルエット値を算出するコードのみ掲載しておく。エルボー図はinertia_をプロットしているだけ。

import codecs
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples, silhouette_score

%matplotlib inline
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt

def select_n_cluster(X):
    # クラスター数2～8を比べる
    range_n_clusters = [i for i in range(2, 9)]    

    sse = []
    for n_clusters in range_n_clusters:
        clusterer = KMeans(init='k-means++', n_clusters=n_clusters, random_state=0)
        cluster_labels = clusterer.fit_predict(X)
        kmeans = clusterer.fit(X)
        
        # SSE（クラスター内誤差の平方和）
        sse.append(kmeans.inertia_) 
        
        # シルエット値（-1～1）の平均
        silhouette_avg = silhouette_score(X, cluster_labels)
        print('For n_clusters =', n_clusters,
              'The average silhouette_score is :', silhouette_avg)    

    # エルボー図のプロット
    plt.plot(range_n_clusters, sse, marker='o')
    plt.xlabel('Number of clusters')
    plt.ylabel('SSE')
    plt.show()

if __name__ == '__main__':
    # csvのデータを読み込んでPandasのDataFrameに
    df = import_estat_csv()

    # 「緑茶」「紅茶」「コーヒー」の合計に対する割合を計算
    df = calc_rate(df)

    # 平均シルエット値の算出とエルボー図のプロット
    select_n_cluster(df.as_matrix())

最適なクラスター数

クラスター数ごとの平均シルエット値では、クラスター数5のときが1に最も近い。
For n_clusters = 2 The average silhouette_score is : 0.6545332827467689 For n_clusters = 3 The average silhouette_score is : 0.7336089740755304 For n_clusters = 4 The average silhouette_score is : 0.8274996031662673 For n_clusters = 5 The average silhouette_score is : 0.9524996088455784 For n_clusters = 6 The average silhouette_score is : 0.8685375857328894 For n_clusters = 7 The average silhouette_score is : 0.7165802427106107 For n_clusters = 8 The average silhouette_score is : 0.6332673269192498

別途シルエット図もプロットしてみたが、クラスター数5のときがクラスターごとのプロット幅が最も均等になった。以下はクラスター数5のときのシルエット図（左）と、「緑茶」「紅茶」の散布図（右）。

エルボー図ではクラスター数5のところがひじの箇所になるので、5が最適クラスター数となる。

シルエット分析でも、エルボー法でも、たまたま選んだクラスター数5が最適という結果になった。

都道府県別の飲み物の購入額をk-means法でクラスタリングする

e-Statには全国消費実態調査データというものがあって、さまざまなものにどのくらいの金額が使われているかがわかる。Python3とそのライブラリのscikit-learnを使って、このデータをk-means法でクラスタリングする。

全国消費実態調査データには集計方法などが違ういくつかのデータセットがあるが、今回使うのは全国消費実態調査 平成21年全国消費実態調査 全国 購入先・購入地域編の［二人以上の世帯］　【品目別1世帯当たり1か月間の支出】　地域，購入地域，購入先別データセットにある「緑茶」「紅茶」「コーヒー」の都道府県別購入金額。

環境

Bash on Ubuntu on Windows
$ cat /etc/issue Ubuntu 16.04.3 LTS \n \l $ python3 -V Python 3.5.2 $ pip3 show scikit-learn pandas matplotlib Name: scikit-learn Version: 0.19.1 ... Name: pandas Version: 0.22.0 ... Name: matplotlib Version: 2.1.0 ...

Jupyter Notebook
$ jupyter --version 4.4.0

データの取得

e-Statからデータを取得する。取得するのは全国消費実態調査 平成21年全国消費実態調査 全国 購入先・購入地域編の［二人以上の世帯］　【品目別1世帯当たり1か月間の支出】　地域，購入地域，購入先別データセット。デフォルトだと項目数が多すぎるので、表示項目選択で項目を絞る。

「2009品目分類」では「緑茶」「紅茶」「コーヒー」を選択し、「2009-47購入地域」と「2009-48購入先」は合計のみ、「2009地域」は全国と都道府県を選択する。表示を更新して、csvをダウンロード。このときに注釈は使わないので「注釈を表示しない」のチェックを外しておく。

csvはsurvey_drink.csvというファイル名で保存しておく。

csvからデータの抽出

ダウンロードしたcsvからPython3で必要なデータを抽出してPandasのDataFrameにする。

import codecs
import pandas as pd

def import_estat_csv():
    with codecs.open('survey_drink.csv', 'r', 'shift_jis', 'ignore') as f:
        # 都道府県をインデックスとして読み込む
        df = pd.read_csv(f, index_col=7, skiprows=10, skipinitialspace=True)

    # 「緑茶」「紅茶」「コーヒー」列のみを抽出
    df = df[df.columns[12:]]

    # 都道府県の行のみ抽出
    return df.iloc[1:48]

if __name__ == '__main__':
    df = import_estat_csv()

    # DataFrameの確認
    print(df)

DataFrameの中身は以下のようになる。 緑茶【円】 紅茶【円】 コーヒー【円】 2009地域 北海道 373 66 473 青森県 232 33 402 岩手県 422 48 399 宮城県 370 38 419 秋田県 412 62 427 山形県 339 37 317 福島県 309 49 329 茨城県 419 60 383 栃木県 429 47 431 ...

クラスタリング用にデータの計算

「緑茶」「紅茶」「コーヒー」の消費金額は都道府県ごとに合計金額がまちまち。このままクラスタリングすると合計金額が多い・少ないかが影響してクラスタリングされてしまうかもしれないので、都道府県ごとに「緑茶」「紅茶」「コーヒー」の合計に対する割合を計算してクラスタリングする。以下はcsvを読み込んで作成したDataFrameで「緑茶」「紅茶」「コーヒー」の合計に対する割合を計算するコード。

import codecs
import pandas as pd

def calc_rate(df):
    df['sum'] = df[df.columns[0]] + df[df.columns[1]] + df[df.columns[2]]
    df['緑茶'] = df[df.columns[0]] / df['sum']
    df['紅茶'] = df[df.columns[1]] / df['sum']
    df['コーヒー'] = df[df.columns[2]] / df['sum']

    # 「緑茶」「紅茶」「コーヒー」の、都道府県ごとの購入金額の割合
    return df[df.columns[4:]]

if __name__ == '__main__':
    df = import_estat_csv()

    df = calc_rate(df)

    print(df)

以下のように各飲み物の都道府県ごとの割合が計算される。
緑茶 紅茶 コーヒー 2009地域 北海道 0.408991 0.072368 0.518640 青森県 0.347826 0.049475 0.602699 岩手県 0.485616 0.055236 0.459148 宮城県 0.447400 0.045949 0.506651 秋田県 0.457270 0.068812 0.473918 山形県 0.489177 0.053391 0.457431 福島県 0.449782 0.071325 0.478894 茨城県 0.486079 0.069606 0.444316 栃木県 0.472988 0.051819 0.475193 ...

k-means法によるクラスタリング

scikit-learnのKMeansクラスを使ってk-means法によりクラスタリングする。k-means法は初期値を設定する必要があるが、KMeansではデフォルトで初期値設定が不要のk-means++法が使われる。また、k-means法ではクラスター数をはじめに決める必要がある。ここではとりあえず５にする。

import codecs
import pandas as pd
from sklearn.cluster import KMeans

def clustering(df):
    # 5つにクラスタリング
    df['cluster'] = KMeans(n_clusters=5, random_state=0).fit_predict(df.as_matrix())

    # クラスタリング結果の確認
    groups = df.groupby('cluster')
    for c, g in groups:
        print('Cluster', c)
        print(g.index.values)

    return groups

if __name__ == '__main__':
    # csvのデータを読み込んでPandasのDataFrameに
    df = import_estat_csv()

    # 「緑茶」「紅茶」「コーヒー」の合計に対する割合を計算
    df = calc_rate(df)

    # クラスタリングと結果表示
    clustering(df)

以下はクラスタリングの結果。緑茶の生産量が多い静岡県、鹿児島県が同じクラスターにに属している。また、東京とその周辺、大阪とその周辺がそれぞれ同じクラスターに属しているので、関東と関西の違いがあるのかもしれない。
  Cluster 0 ['石川県' '奈良県' '和歌山県' '鳥取県' '岡山県' '広島県' '徳島県' '愛媛県' '高知県' '沖縄県'] Cluster 1 ['静岡県' '佐賀県' '長崎県' '鹿児島県'] Cluster 2 ['岩手県' '秋田県' '山形県' '福島県' '茨城県' '栃木県' '埼玉県' '千葉県' '東京都' '神奈川県' '新潟県' '長野県' '熊本県' '宮崎県'] Cluster 3 ['北海道' '宮城県' '群馬県' '福井県' '三重県' '滋賀県' '島根県' '福岡県' '大分県'] Cluster 4 ['青森県' '富山県' '山梨県' '岐阜県' '愛知県' '京都府' '大阪府' '兵庫県' '山口県' '香川県']

クラスタリング結果のグラフ化

文字だけ見てもどういう結果になったか分かりにくいので、クラスターごとの平均値をmatplotlibでグラフ化する。

import codecs
import pandas as pd
from sklearn.cluster import KMeans

%matplotlib inline
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt

def draw_bar(groups):
    means = groups.mean()
    errors = groups.std()

    fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(1,1,1)
    means.plot.bar(color=['green', 'orange', 'brown'], yerr=errors, ax=ax)
    ax.set_xticklabels(means.index.values, rotation=0)
    plt.show()

if __name__ == '__main__':
    # csvのデータを読み込んでPandasのDataFrameに
    df = import_estat_csv()

    # 「緑茶」「紅茶」「コーヒー」の合計に対する割合を計算
    df = calc_rate(df)

    # クラスタリングと結果表示
    groups = clustering(df)

    # クラスターごとの平均値の棒グラフ
    draw_bar(groups)

クラスターごとの平均購入額割合

静岡県が属するクラスター１の緑茶の購入額割合が多いのは納得。一方、クラスター４に属する京都府は意外と緑茶の購入額割合が多くない。また、首都圏が属するクラスター２は緑茶とコーヒーの購入額割合が拮抗していており、関西圏が属するクラスター４はコーヒーの購入額割合が緑茶に比べてかなり多い。首都圏と関西圏では傾向に違いがあることがわかる。