sage/オープンデータ解析例１

FrontPage

2016/02/16からのアクセス回数 3231

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/61/

また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

オープンデータ解析例１ †

高岡市の公開しているオープンデータを使ってデータ解析をしてみます。

今回は、平成27年7月地区別世帯数及び人口集計表の「全体」シートを使って、高岡市を２つのグループに分類してみます。

必要なユーティリティの準備 †

データを簡単に処理するために、ggplot2等のライブラリーを読み込みます。

sageへの入力:

# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を追加+GgSaveを追加（2015/07/12）
load(DATA + 'RUtil.py')

平成27年7月地区別世帯数及び人口集計表 †

それでは、オープンデータとして公開されている高岡市の「平成27年7月地区別世帯数及び人口集計表」のデータを 読み込んでPandasのデータフレームに取り込みます。

平成27年7月地区別世帯数及び人口集計表のデータは、以下の様に入力されています。

ヘッダは、4,5行目にあり、カラムは、2列目から始まっています。

pd.read_excel関数に読み込むExcelファイルのあるURLを指定します。 ヘッダは、3（0始まり）とし、カラムは1（0始まり）にあることを指定するため、header=3, index_col=1を引数にセットします。

sageへの入力:

# 高岡市のオープンデータから平成27年7月の人口データ(h270731.xlsx)を取り込みます
d = pd.read_excel('http://www.city.takaoka.toyama.jp/joho/shise/opendata/documents/h270731.xlsx', header=3, index_col=1)
d.head()

      頁  地区\nコード  地区名  人　　口 Unnamed: 4 Unnamed: 5   世帯数
NaN NaN      NaN  NaN     男          女          計   NaN
 2  NaN  01 平米地区  NaN  1573       1807       3380  1591
 3  NaN  02 定塚地区  NaN  4428       4967       9395  3977
 4  NaN  03 下関地区  NaN  4439       4629       9068  3971
 5  NaN  04 博労地区  NaN  5317       5833      11150  4582

カラム名のセットと不要なカラムの削除 †

カラム名が日本語のままだと操作しにくいので、カラム名をd.columnsに英語名でセットします。 頁の後にNaNという変な値が入っているので、not_usedというカラムとしました。 地区名はセル結合されて、codeの欄に地区コード+地区名でセットされていました。

sageへの入力:

# カラム名を付け直し、最初の5個を表示
d.columns = ['not_used', 'code', 'region', 'male', 'female', 'population', 'household']
d.head()

     not_used     code  region  male female population  household
NaN       NaN      NaN     NaN     男      女          計        NaN
 2        NaN  01 平米地区     NaN  1573   1807       3380       1591
 3        NaN  02 定塚地区     NaN  4428   4967       9395       3977
 4        NaN  03 下関地区     NaN  4439   4629       9068       3971
 5        NaN  04 博労地区     NaN  5317   5833      11150       4582

博労地区と佐野地区は木津地区を除いた集計がされているので、indexがNaNとなっており、 d.index > 1とすることで不要なデータを取り除きます。

sageへの入力:

# 最初の番号（index）が1より大きなcode, population, householdを抽出します。
d1 = d[d.index > 1][['code', 'male', 'female', 'population', 'household']]
d1.head()

      code  male female population  household
2  01 平米地区  1573   1807       3380       1591
3  02 定塚地区  4428   4967       9395       3977
4  03 下関地区  4439   4629       9068       3971
5  04 博労地区  5317   5833      11150       4582
6  05 横田地区  2839   3070       5909       2286

データの可視化 †

各要素にどのような関係があるかを、Rのpairsプロットでみてみましょう。

SageのPandaで作られたデータフレームをRのデータフレームに変換します。 残念ながら、日本語のデータが正しくRに変換できないため、ここでは地域 codeを除いてRに渡しています。

sageへの入力:

# d2をRのデータフレームd1に変換する
# 日本語が上手く渡せない
d2 = d[d.index > 1][['male', 'female', 'population', 'household']]
d2 = d2[d2.index <36]
PandaDf2RDf(d2, "d2")

sageへの入力:

# 分布を可視化
g_pairs = preGraph("pairs.pdf")
r('pairs(d2, main="Basic Scatter Plot Matrix")')
postGraph(g_pairs, fac=0.6)

世帯数（household）と人口（population）の関係を散布図でみてみましょう。

世帯数と人口は、概ね直線上にのっており、地区によるおおきな変化はみられません。

sageへの入力:

# 世帯数（household）と人口（population）の関係をみる
df = d1[d1.index <36]
ggplot(df, aes(x='household', y='population')) + geom_point()

<repr(<ggplot.ggplot.ggplot at 0x7f7943c6a6d0>) failed: KeyError: 0.0>

sageへの入力:

GgSave('fig_1.png', fac=0.8)

Saving 11.0 x 8.0 in image.

グループ分け †

高岡市の地区を2つのグループに分けてみましょう。 グループ分けに使用するデータは、男性、女性の人数と世帯数のデータとしました。

分析に使用する手法は、混合正規分布（GMM）を使用します。 ここでは、sklearnのmixtureパッケージのGMMを使っています。 n_components=2で2つのグループに分けることを指定しています。

sageへの入力:

# 混合正規分布
from sklearn import mixture
# 分類器の生成
classifier = mixture.GMM(n_components=2, covariance_type='full')

sageへの入力:

# GMMを求める
# GMMについては、http://www15191ue.sakura.ne.jp:8000/home/pub/57/ を参照してください。
data = df[['male', 'female', 'household']].values
# 混合正規分布を求める
classifier.fit(data)

GMM(covariance_type='full', init_params='wmc', min_covar=0.001,
  n_components=2, n_init=1, n_iter=100, params='wmc', random_state=None,
  thresh=None, tol=0.001)

分析結果のプロット †

グループ分けされた結果は、classifier.predict関数で得ることができます。 この結果をデータフレームdfのpredにセットします。

predの値で散布図の点を色分けしてプロットしてみます。

sageへの入力:

# 分類結果をpredにセット
df['pred'] = classifier.predict(data)
# 世帯数（household）と人口（population）の関係をみる
ggplot(df, aes(x='population', y='household', colour='pred')) + geom_point()

<repr(<ggplot.ggplot.ggplot at 0x7f79437b9810>) failed: KeyError: 0.0>

sageへの入力:

GgSave('fig_2.png', fac=0.8)

Saving 11.0 x 8.0 in image.

世帯数と分類の結果 †

データを世帯数でソートして表示してみます。 predがグループ分けの結果で、平米地区が世帯数が中田地区よりも少ないのに グループ１に入っており、興味深い結果となっています。

sageへの入力:

# 世帯数でソート
# 平米地区に注目
df.sort('household')

           code  male female population  household  pred
33  30 福岡町五位山地区   204    217        421        170     0
18      16 小勢地区   488    498        986        277     0
34   31 福岡町赤丸地区   577    612       1189        414     0
21      19 石堤地区   613    619       1232        443     0
32  29 福岡町西五位地区   938    975       1913        613     0
31   28 福岡町大滝地区  1197   1240       2437        796     0
25      23 太田地区  1221   1387       2608        889     0
17      15 福田地区  1309   1347       2656        900     0
23      21 守山地区  1259   1380       2639        933     0
10      09 二上地区  1270   1366       2636        999     0
29   26 福岡町福岡地区  1477   1616       3093       1153     0
19      17 立野地区  1578   1698       3276       1181     0
22      20 国吉地区  1719   1844       3563       1192     0
15      13 二塚地区  1770   1894       3664       1218     0
30   27 福岡町山王地区  1958   1980       3938       1244     0
8       07 川原地区  1554   1719       3273       1336     0
2       01 平米地区  1573   1807       3380       1591     1
20     18 東五位地区  2354   2546       4900       1700     0
35       木津地区　※  2765   2966       5731       2135     0
28      25 中田地区  3025   3270       6295       2220     0
6       05 横田地区  2839   3070       5909       2286     0
16      14 佐野地区  3238   3530       6768       2399     0
7       06 西条地区  3427   3565       6992       2676     1
12      11 牧野地区  4678   4909       9587       3461     1
9       08 成美地区  3877   4414       8291       3512     1
4       03 下関地区  4439   4629       9068       3971     1
3       02 定塚地区  4428   4967       9395       3977     1
24      22 伏木地区  5368   5939      11307       4442     1
11      10 能町地区  5667   5874      11541       4449     1
5       04 博労地区  5317   5833      11150       4582     1
26      24 戸出地区  6716   6921      13637       4948     1
13      12 野村地区  8443   8996      17439       6907     1

男性と女性の数の関係 †

同様に、男性と女性の関係をプロットしてみます。 世帯数よりも男性と女性の関係の方が、はっきりした直線の関係がみられます。

sageへの入力:

# 男性（male）と女性（female）の関係をみる
ggplot(df, aes(x='male', y='female', colour='pred')) + geom_point()

<repr(<ggplot.ggplot.ggplot at 0x7f793f22b710>) failed: KeyError: 0.0>

sageへの入力:

GgSave('fig_3.png', fac=0.8)

Saving 11.0 x 8.0 in image.

ペア・プロットにグループで色分け †

平米地区がグループ１に入っている理由をみるために、 ペア・プロットにグループで色分けしてみます。

sageへの入力:

# グループ分けの結果で色分けしてペア・プロットを表示
d2 = df[['male', 'female', 'population', 'household', 'pred']]
PandaDf2RDf(d2, "d2")

sageへの入力:

# 分布を可視化
graph = preGraph("pairs2.pdf")
r('pairs(d2[, c(2, 3, 4, 5)], main="Takaoka Group", pch=21, bg=c("red", "blue")[factor(d2$pred)])')
postGraph(graph, fac=0.8)

クラス分けを地図に表示 †

国土地理院の小学校区のshapeファイル †

クラス分けの結果を地図にプロットするために、国土地理院が提供している富山県の小学校区のデータ（A27-10_16_GML）を使用しました。

• 国土交通省 GIS データのダウンロード

小学校区は、A27-10_16-g_SchoolDistrict.shpに定義され、小学校名はA27-10_16-g_SchoolDistrict.dbfに定義されています。 国土地理院のデータはShift-JISコードであり、これをUTF-8に変換する必要があります。

以下のサイトを参考にUTF-8に変換しました。

• Google Earth の使い方～その２　道路や鉄道、将来人口情報などの表示

地区と小学校の関係 †

高岡市の地区と小学校区の関係をExcelシートに整理しました。市町村合併や小学校の統合によって、 分からなかったので、高岡市の防災マップを元に人に聞きながら作ったので、間違ってかも知れませんが、 このExcelファイルを修正すれば、地図はすぐに書き直せます。

sageへの入力:

# 高岡市地区の小学校区を読み込む
schoolDistrict = pd.read_excel(DATA + "schoolDistrict.xls")
# 予測結果に対する色と小学校をd2にセット
schoolDistrict.head()

      code schoolDsitrict
0  01 平米地区          平米小学校
1  02 定塚地区          定塚小学校
2  03 下関地区          下関小学校
3  04 博労地区          博労小学校
4  05 横田地区          横田小学校

つぎにcodeをキーにして、高岡市のクラス分け結果（df）と小学校区（schoolDistrict）をマージします。

sageへの入力:

# 小学校区と地区をマージ
df.set_index('code')
schoolDistrict.set_index('code')
df = df.merge(schoolDistrict)
df.head()

      code  male female population  household  pred schoolDsitrict
0  01 平米地区  1573   1807       3380       1591     1          平米小学校
1  02 定塚地区  4428   4967       9395       3977     1          定塚小学校
2  03 下関地区  4439   4629       9068       3971     1          下関小学校
3  04 博労地区  5317   5833      11150       4582     1          博労小学校
4  05 横田地区  2839   3070       5909       2286     0          横田小学校

地図の表示には、Rのmaptoolsを使用するため、 小学校区（schoolDsitrict）とクラス分け結果（pred）をCSVファイルとして出力します。

sageへの入力:

# 小学校区のクラス分け
schoolClass = df[['schoolDsitrict', 'pred']]. drop_duplicates()
schoolClass.head()

  schoolDsitrict  pred
0          平米小学校     1
1          定塚小学校     1
2          下関小学校     1
3          博労小学校     1
4          横田小学校     0

sageへの入力:

# 日本語を含むため、csvファイルを使ってRへ渡す
schoolClass_file = DATA + 'schoolCalss.csv'
schoolClass.to_csv(schoolClass_file)

maptoolsで地図にプロット †

Rのmaptoolsをロードし、富山県の小学校区のshapeファイルA27-10_16-g_SchoolDistrict.shpをshapeに読み込みます。

shape変数から高岡市立の小学校区のみを抽出するため、shape[shape@data$A27_006=="高岡市立",]の条件抽出を使います。

sageへの入力:

# maptoolsをロード
# r("install.packages('maptools')")
r("library(maptools)")

 [1] "maptools"  "sp"        "jsonlite"  "ggplot2"   "stats"     "graphics"  "grDevices" "utils"   
 [9] "datasets"  "methods"   "base"    

sageへの入力:

# 富山県の小学校区のshapeファイルを読み込み、高岡市立の小学校のみを取り出す
shape_file = DATA + "A27-10_16-g_SchoolDistrict.shp"
r('shape <- readShapePoly("%s")' % shape_file)
r('takaoka <- shape[shape@data$A27_006=="高岡市立",]')
r('takaoka@data$A27_007')

 [1] 成美小学校     戸出東部小学校 戸出東部小学校 博労小学校     戸出西部小学校 伏木小学校   
 [7] 万葉小学校     福岡小学校     太田小学校     国吉小学校     牧野小学校     千鳥丘小学校 
[13] 西広谷小学校   古府小学校     能町小学校     西条小学校     野村小学校     川原小学校   
[19] 石堤小学校     定塚小学校     横田小学校     平米小学校     東五位小学校   木津小学校   
[25] 南条小学校     二塚小学校     中田小学校     下関小学校     下関小学校     成美小学校   
122 Levels: あさひ野小学校 さみさと小学校 井口小学校 井波小学校 宇奈月小学校 ... 立山北部小学校

高岡市の小学校区のクラス分け結果のCSVファイルを読み込み、takaoka@data$A27_007に含まれる小学校のクラス結果predをセットします。

色分けを0=赤、１=青とし、Rのインデックスが１はじまりなので、pred+1としてtakaoka@data$A27_007に含まれる小学校の色をcolにセットします。

sageへの入力:

# shcoolClassを読み込む
r('shcoolClass <- read.csv("%s")' % schoolClass_file)
# クラス分けを変数predにセット
r('pred <- shcoolClass$pred[match(takaoka@data$A27_007 , shcoolClass$schoolDsitrict)]')
# クラスの色を配列にセット
r('cols <- c("red", "blue")')
# クラス分けを色に変換（Rのインデックスが1始まりなのでpred+1とする）
r('col <- cols[pred+1]')

 [1] "blue" "blue" "blue" "blue" "blue" "blue" "red"  "red"  "red"  "red"  "blue" "red"  NA   
[14] "blue" "blue" "blue" "blue" "red"  "red"  "blue" "red"  "blue" "red"  "red"  "red"  "red"
[27] "red"  "blue" "blue" "blue"

準備ができたので、takaokaの小学校区のクラス分けを地図にプロットします。

赤で囲まれた青の地区は、「川原地区」です。

sageへの入力:

graph = preGraph("schoolDistrict.pdf")
r('plot(takaoka, col=col)')
postGraph(graph, fac=0.8)

データ替えて同じ計算をする †

最初のデータの読み込みを以下のように6月のデータに替えて、ActionメニューからEvaluate Allを選択すると、 平成27年6月地区別世帯数及び人口集計表に対して同じ計算を行い、その結果をグラフに表示します。

     d = pd.read_excel('http://www.city.takaoka.toyama.jp/joho/shise/opendata/documents/h270630.xlsx', header=3, index_col=1)

これが、Sageを使った「リプロデューシブルリサーチ」です。 Sageは、実行可能なドキュメントなので、入力のExcelファイルを替えて「Evaluate All」を実行すれば、 別のデータに同じような分析を実行することができるのです。

コメント †

選択肢	投票
おもしろかった	0
そうでもない	0
わかりずらい	0

皆様のご意見、ご希望をお待ちしております。

お名前:

(Input image string)