OpenData_Analysis_01

1318 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

オープンデータ解析例1

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

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

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

データを簡単に処理するために、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がグループ分けの結果で、平米地区が世帯数が中田地区よりも少ないのに グループ1に入っており、興味深い結果となっています。

# 世帯数でソート # 平米地区に注目 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.

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

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

# グループ分けの結果で色分けしてペア・プロットを表示 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=赤、1=青とし、Rのインデックスが1はじまりなので、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 
html('<a href="takaoka_jinko.html">高岡の人口分布</a>') 

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

最初のデータの読み込みを以下のように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」を実行すれば、 別のデータに同じような分析を実行することができるのです。