Hiroshi TAKEMOTO (take@pwv.co.jp)

オープンデータ解析例１

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

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

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

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

{{{id=1| # 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を引数にセットします。

{{{id=2| # 高岡市のオープンデータから平成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の欄に地区コード+地区名でセットされていました。

{{{id=4| # カラム名を付け直し、最初の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とすることで不要なデータを取り除きます。

{{{id=7| # 最初の番号（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に渡しています。

{{{id=28| # d2をRのデータフレームd1に変換する # 日本語が上手く渡せない d2 = d[d.index > 1][['male', 'female', 'population', 'household']] d2 = d2[d2.index <36] PandaDf2RDf(d2, "d2") /// }}} {{{id=27| # 分布を可視化 g_pairs = preGraph("pairs.pdf") r('pairs(d2, main="Basic Scatter Plot Matrix")') postGraph(g_pairs, fac=0.6) /// }}}

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

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

{{{id=9| # 世帯数（household）と人口（population）の関係をみる df = d1[d1.index <36] ggplot(df, aes(x='household', y='population')) + geom_point() /// ) failed: KeyError: 0.0> }}} {{{id=10| GgSave('fig_1.png', fac=0.8) /// Saving 11.0 x 8.0 in image. }}}

グループ分け

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

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

{{{id=12| # 混合正規分布 from sklearn import mixture # 分類器の生成 classifier = mixture.GMM(n_components=2, covariance_type='full') /// }}} {{{id=13| # 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の値で散布図の点を色分けしてプロットしてみます。

{{{id=14| # 分類結果をpredにセット df['pred'] = classifier.predict(data) # 世帯数（household）と人口（population）の関係をみる ggplot(df, aes(x='population', y='household', colour='pred')) + geom_point() /// ) failed: KeyError: 0.0> }}} {{{id=15| GgSave('fig_2.png', fac=0.8) /// Saving 11.0 x 8.0 in image. }}}

世帯数と分類の結果

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

{{{id=25| # 世帯数でソート # 平米地区に注目 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 }}}

男性と女性の数の関係

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

{{{id=16| # 男性（male）と女性（female）の関係をみる ggplot(df, aes(x='male', y='female', colour='pred')) + geom_point() /// ) failed: KeyError: 0.0> }}} {{{id=20| GgSave('fig_3.png', fac=0.8) /// Saving 11.0 x 8.0 in image. }}}

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

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

{{{id=30| # グループ分けの結果で色分けしてペア・プロットを表示 d2 = df[['male', 'female', 'population', 'household', 'pred']] PandaDf2RDf(d2, "d2") /// }}} {{{id=31| # 分布を可視化 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ファイルを修正すれば、地図はすぐに書き直せます。

{{{id=39| # 高岡市地区の小学校区を読み込む schoolDistrict = pd.read_excel(DATA + "schoolDistrict.xls") # 予測結果に対する色と小学校をd2にセット schoolDistrict.head() /// code schoolDsitrict 0 01 平米地区 平米小学校 1 02 定塚地区 定塚小学校 2 03 下関地区 下関小学校 3 04 博労地区 博労小学校 4 05 横田地区 横田小学校 }}}

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

{{{id=62| # 小学校区と地区をマージ 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ファイルとして出力します。

{{{id=63| # 小学校区のクラス分け 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 }}} {{{id=64| # 日本語を含むため、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=="高岡市立",]の条件抽出を使います。

{{{id=40| # maptoolsをロード # r("install.packages('maptools')") r("library(maptools)") /// [1] "maptools" "sp" "jsonlite" "ggplot2" "stats" "graphics" "grDevices" "utils" [9] "datasets" "methods" "base" }}} {{{id=42| # 富山県の小学校区の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にセットします。

{{{id=65| # 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" }}}

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

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

{{{id=67| graph = preGraph("schoolDistrict.pdf") r('plot(takaoka, col=col)') postGraph(graph, fac=0.8) /// }}} {{{id=69| # ここからは製作途中 # 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 }}} {{{id=71| # Sageではエラーになる #graph = preGraph("JINKO.pdf") #r('spplot(takaoka, zcol="JINKO", col.regions=terrain.colors(length(brks)-1), at=brks)') #postGraph(graph, fac=0.8) /// }}} {{{id=75| # 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" }}} {{{id=77| # 高岡の人口分布を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 }}} {{{id=78| html('高岡の人口分布') /// 高岡の人口分布 }}}

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

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

{{{id=79| /// }}}