オープンデータ解析例１
高岡市の公開しているオープンデータを使ってデータ解析をしてみます。
今回は、平成27年7月地区別世帯数及び人口集計表の「全体」シートを使って、高岡市を２つのグループに分類してみます。
必要なユーティリティの準備
データを簡単に処理するために、ggplot2等のライブラリーを読み込みます。
# 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を引数にセットします。
# 高岡市のオープンデータから平成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 頁 地区\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の欄に地区コード+地区名でセットされていました。
# カラム名を付け直し、最初の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 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とすることで不要なデータを取り除きます。
# 最初の番号（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 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に渡しています。
# d2をRのデータフレームd1に変換する
# 日本語が上手く渡せない
d2 = d[d.index > 1][['male', 'female', 'population', 'household']]
d2 = d2[d2.index <36]
PandaDf2RDf(d2, "d2")
|
|
# 分布を可視化
g_pairs = preGraph("pairs.pdf")
r('pairs(d2, main="Basic Scatter Plot Matrix")')
postGraph(g_pairs, fac=0.6)
|
|
世帯数（household）と人口（population）の関係を散布図でみてみましょう。
世帯数と人口は、概ね直線上にのっており、地区によるおおきな変化はみられません。
# 世帯数（household）と人口（population）の関係をみる
df = d1[d1.index <36]
ggplot(df, aes(x='household', y='population')) + geom_point()
<repr(<ggplot.ggplot.ggplot at 0x7f6ff54ba6d0>) failed: KeyError: 0.0> <repr(<ggplot.ggplot.ggplot at 0x7f6ff54ba6d0>) failed: KeyError: 0.0> |
GgSave('fig_1.png', fac=0.8)
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
グループ分け
高岡市の地区を2つのグループに分けてみましょう。 グループ分けに使用するデータは、男性、女性の人数と世帯数のデータとしました。
分析に使用する手法は、混合正規分布（GMM）を使用します。 ここでは、sklearnのmixtureパッケージのGMMを使っています。 n_components=2で2つのグループに分けることを指定しています。
# 混合正規分布
from sklearn import mixture
# 分類器の生成
classifier = mixture.GMM(n_components=2, covariance_type='full')
|
|
# 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) 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の値で散布図の点を色分けしてプロットしてみます。
# 分類結果をpredにセット
df['pred'] = classifier.predict(data)
# 世帯数（household）と人口（population）の関係をみる
ggplot(df, aes(x='population', y='household', colour='pred')) + geom_point()
<repr(<ggplot.ggplot.ggplot at 0x7f6ff5009810>) failed: KeyError: 0.0> <repr(<ggplot.ggplot.ggplot at 0x7f6ff5009810>) failed: KeyError: 0.0> |
GgSave('fig_2.png', fac=0.8)
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
世帯数と分類の結果
データを世帯数でソートして表示してみます。 predがグループ分けの結果で、平米地区が世帯数が中田地区よりも少ないのに グループ１に入っており、興味深い結果となっています。
# 世帯数でソート
# 平米地区に注目
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 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 |
男性と女性の数の関係
同様に、男性と女性の関係をプロットしてみます。 世帯数よりも男性と女性の関係の方が、はっきりした直線の関係がみられます。
# 男性（male）と女性（female）の関係をみる
ggplot(df, aes(x='male', y='female', colour='pred')) + geom_point()
<repr(<ggplot.ggplot.ggplot at 0x7f6ff0a7b710>) failed: KeyError: 0.0> <repr(<ggplot.ggplot.ggplot at 0x7f6ff0a7b710>) failed: KeyError: 0.0> |
GgSave('fig_3.png', fac=0.8)
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
ペア・プロットにグループで色分け
平米地区がグループ１に入っている理由をみるために、 ペア・プロットにグループで色分けしてみます。
# グループ分けの結果で色分けしてペア・プロットを表示
d2 = df[['male', 'female', 'population', 'household', 'pred']]
PandaDf2RDf(d2, "d2")
|
|
# 分布を可視化
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）を使用しました。
小学校区は、A27-10_16-g_SchoolDistrict.shpに定義され、小学校名はA27-10_16-g_SchoolDistrict.dbfに定義されています。 国土地理院のデータはShift-JISコードであり、これをUTF-8に変換する必要があります。 以下のサイトを参考にUTF-8に変換しました。
地区と小学校の関係
高岡市の地区と小学校区の関係をExcelシートに整理しました。市町村合併や小学校の統合によって、 分からなかったので、高岡市の防災マップを元に人に聞きながら作ったので、間違ってかも知れませんが、 このExcelファイルを修正すれば、地図はすぐに書き直せます。
# 高岡市地区の小学校区を読み込む
schoolDistrict = pd.read_excel(DATA + "schoolDistrict.xls")
# 予測結果に対する色と小学校をd2にセット
schoolDistrict.head()
code schoolDsitrict 0 01 平米地区 平米小学校 1 02 定塚地区 定塚小学校 2 03 下関地区 下関小学校 3 04 博労地区 博労小学校 4 05 横田地区 横田小学校 code schoolDsitrict 0 01 平米地区 平米小学校 1 02 定塚地区 定塚小学校 2 03 下関地区 下関小学校 3 04 博労地区 博労小学校 4 05 横田地区 横田小学校 |
つぎにcodeをキーにして、高岡市のクラス分け結果（df）と小学校区をマージします。
# 小学校区と地区をマージ
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 横田小学校 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ファイルとして出力します。
# 小学校区のクラス分け
schoolClass = df[['schoolDsitrict', 'pred', 'population']]. drop_duplicates()
schoolClass.head()
schoolDsitrict pred population 0 平米小学校 1 3380 1 定塚小学校 1 9395 2 下関小学校 1 9068 3 博労小学校 1 11150 4 横田小学校 0 5909 schoolDsitrict pred population 0 平米小学校 1 3380 1 定塚小学校 1 9395 2 下関小学校 1 9068 3 博労小学校 1 11150 4 横田小学校 0 5909 |
# 日本語を含むため、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=="高岡市立",]の条件抽出を使います。
# maptoolsをロード
# r("install.packages('maptools')")
r("library(maptools)")
[1] "maptools" "sp" "jsonlite" "ggplot2" "stats" "graphics" "grDevices" "utils" [9] "datasets" "methods" "base" [1] "maptools" "sp" "jsonlite" "ggplot2" "stats" "graphics" "grDevices" "utils" [9] "datasets" "methods" "base" |
# 富山県の小学校区の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: あさひ野小学校 さみさと小学校 井口小学校 井波小学校 宇奈月小学校 ... 立山北部小学校 [1] 成美小学校 戸出東部小学校 戸出東部小学校 博労小学校 戸出西部小学校 伏木小学校 [7] 万葉小学校 福岡小学校 太田小学校 国吉小学校 牧野小学校 千鳥丘小学校 [13] 西広谷小学校 古府小学校 能町小学校 西条小学校 野村小学校 川原小学校 [19] 石堤小学校 定塚小学校 横田小学校 平米小学校 東五位小学校 木津小学校 [25] 南条小学校 二塚小学校 中田小学校 下関小学校 下関小学校 成美小学校 122 Levels: あさひ野小学校 さみさと小学校 井口小学校 井波小学校 宇奈月小学校 ... 立山北部小学校 |
高岡市の小学校区のクラス分け結果のCSVファイルを読み込み、takaoka@data$A27_007に含まれる小学校のクラス結果predをセットします。
色分けを0=赤、１=青とし、Rのインデックスが１はじまりなので、pred+1としてtakaoka@data$A27_007に含まれる小学校の色をcolにセットします。
# shcoolClassを読み込む
r('shcoolClass <- read.csv("%s")' % schoolClass_file)
# クラス分けを変数predにセット
r('idx <- match(takaoka@data$A27_007 , shcoolClass$schoolDsitrict)')
r('pred <- shcoolClass$pred[idx]')
# 小学校区毎の人口をJINKOにセット
r('takaoka@data$JINKO <- shcoolClass$population[idx]')
# クラスの色を配列にセット
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" [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の小学校区のクラス分けを地図にプロットします。
赤で囲まれた青の地区は、「川原地区」です。
graph = preGraph("schoolDistrict.pdf")
r('plot(takaoka, col=col)')
postGraph(graph, fac=0.8)
|
|
# ここからは製作途中
#
r('library(sp)')
# 欠損値を許さない naに0をセット
r('takaoka@data$JINKO[is.na(takaoka@data$JINKO)] <- 0')
# 分位点で色分け
r('brks <- round(quantile(takaoka@data$JINKO, 0:10/10))')
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
0 2633 3359 4529 5838 6880 8602 9453 11307 13637 17439
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
0 2633 3359 4529 5838 6880 8602 9453 11307 13637 17439
|
# Sageではエラーになる
#graph = preGraph("JINKO.pdf")
#r('spplot(takaoka, zcol="JINKO", col.regions=terrain.colors(length(brks)-1), at=brks)')
#postGraph(graph, fac=0.8)
|
|
# plotGoogleMapsをロードする。CentOSへのインストールには、gdal-devel proj-devel proj-epsg proj-nad が必要
r("library(plotGoogleMaps)")
[1] "plotGoogleMaps" "spacetime" "maptools" "sp" "jsonlite" [6] "ggplot2" "stats" "graphics" "grDevices" "utils" [11] "datasets" "methods" "base" [1] "plotGoogleMaps" "spacetime" "maptools" "sp" "jsonlite" [6] "ggplot2" "stats" "graphics" "grDevices" "utils" [11] "datasets" "methods" "base" |
# 高岡の人口分布をtakaoka_jinko.htmlでGoogle Mapと一緒に表示
html_file = DATA + 'takaoka_jinko.html'
r('proj4string(takaoka) <- CRS("+proj=longlat +ellps=WGS84")')
r('m <- plotGoogleMaps(takaoka, zcol="JINKO", layerName="population", filename="%s", colPalette=terrain.colors(length(brks)-1), at=brks, strokeColor=1, mapTypeId="ROADMAP", openMap=FALSE)' % html_file)
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
0 2633 3359 4529 5838 6880 8602 9453 11307 13637 17439
0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
0 2633 3359 4529 5838 6880 8602 9453 11307 13637 17439
|
データ替えて同じ計算をする
最初のデータの読み込みを以下のように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」を実行すれば、 別のデータに同じような分析を実行することができるのです。
|
|