![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
[[FrontPage]] #contents 2016/02/16からのアクセス回数 &counter; ここで紹介したSageワークシートは、以下のURLからダウンロードできます。 http://www15191ue.sakura.ne.jp:8000/home/pub/61/ また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。 ** オープンデータ解析例１ [#k120014b] 高岡市の公開しているオープンデータを使ってデータ解析をしてみます。 今回は、平成27年7月地区別世帯数及び人口集計表の「全体」シートを使って、高岡市を２つのグループに分類してみます。 *** 必要なユーティリティの準備 [#rd844635] データを簡単に処理するために、ggplot2等のライブラリーを読み込みます。 sageへの入力: #pre{{ # RとPandasで必要なユーティリティ # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # RUtilにRとPandasのデータフレームを相互に変換する関数を追加+GgSaveを追加（2015/07/12） load(DATA + 'RUtil.py') }} ** 平成27年7月地区別世帯数及び人口集計表 [#y88d0032] それでは、オープンデータとして公開されている高岡市の「平成27年7月地区別世帯数及び人口集計表」のデータを 読み込んでPandasのデータフレームに取り込みます。 平成27年7月地区別世帯数及び人口集計表のデータは、以下の様に入力されています。 &ref(Excel_format.jpg); ヘッダは、4,5行目にあり、カラムは、2列目から始まっています。 pd.read_excel関数に読み込むExcelファイルのあるURLを指定します。 ヘッダは、3（0始まり）とし、カラムは1（0始まり）にあることを指定するため、header=3, index_col=1を引数にセットします。 sageへの入力: #pre{{ # 高岡市のオープンデータから平成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() }} #pre{{ 頁 地区\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 }} *** カラム名のセットと不要なカラムの削除 [#pde60d63] カラム名が日本語のままだと操作しにくいので、カラム名をd.columnsに英語名でセットします。 頁の後にNaNという変な値が入っているので、not_usedというカラムとしました。 地区名はセル結合されて、codeの欄に地区コード+地区名でセットされていました。 sageへの入力: #pre{{ # カラム名を付け直し、最初の5個を表示 d.columns = ['not_used', 'code', 'region', 'male', 'female', 'population', 'household'] d.head() }} #pre{{ 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への入力: #pre{{ # 最初の番号（index）が1より大きなcode, population, householdを抽出します。 d1 = d[d.index > 1][['code', 'male', 'female', 'population', 'household']] d1.head() }} #pre{{ 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 }} *** データの可視化 [#w55636ee] 各要素にどのような関係があるかを、Rのpairsプロットでみてみましょう。 SageのPandaで作られたデータフレームをRのデータフレームに変換します。 残念ながら、日本語のデータが正しくRに変換できないため、ここでは地域 codeを除いてRに渡しています。 sageへの入力: #pre{{ # d2をRのデータフレームd1に変換する # 日本語が上手く渡せない d2 = d[d.index > 1][['male', 'female', 'population', 'household']] d2 = d2[d2.index <36] PandaDf2RDf(d2, "d2") }} sageへの入力: #pre{{ # 分布を可視化 g_pairs = preGraph("pairs.pdf") r('pairs(d2, main="Basic Scatter Plot Matrix")') postGraph(g_pairs, fac=0.6) }} &ref(th_pairs.png); 世帯数（household）と人口（population）の関係を散布図でみてみましょう。 世帯数と人口は、概ね直線上にのっており、地区によるおおきな変化はみられません。 sageへの入力: #pre{{ # 世帯数（household）と人口（population）の関係をみる df = d1[d1.index <36] ggplot(df, aes(x='household', y='population')) + geom_point() }} #pre{{ <repr(<ggplot.ggplot.ggplot at 0x7f7943c6a6d0>) failed: KeyError: 0.0> }} sageへの入力: #pre{{ GgSave('fig_1.png', fac=0.8) }} #pre{{ Saving 11.0 x 8.0 in image. }} &ref(th_fig_1.png); ** グループ分け [#sefb380b] 高岡市の地区を2つのグループに分けてみましょう。 グループ分けに使用するデータは、男性、女性の人数と世帯数のデータとしました。 分析に使用する手法は、混合正規分布（GMM）を使用します。 ここでは、sklearnのmixtureパッケージのGMMを使っています。 n_components=2で2つのグループに分けることを指定しています。 sageへの入力: #pre{{ # 混合正規分布 from sklearn import mixture # 分類器の生成 classifier = mixture.GMM(n_components=2, covariance_type='full') }} sageへの入力: #pre{{ # GMMを求める # GMMについては、http://www15191ue.sakura.ne.jp:8000/home/pub/57/ を参照してください。 data = df[['male', 'female', 'household']].values # 混合正規分布を求める classifier.fit(data) }} #pre{{ 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) }} *** 分析結果のプロット [#oe1932f1] グループ分けされた結果は、classifier.predict関数で得ることができます。 この結果をデータフレームdfのpredにセットします。 predの値で散布図の点を色分けしてプロットしてみます。 sageへの入力: #pre{{ # 分類結果をpredにセット df['pred'] = classifier.predict(data) # 世帯数（household）と人口（population）の関係をみる ggplot(df, aes(x='population', y='household', colour='pred')) + geom_point() }} #pre{{ <repr(<ggplot.ggplot.ggplot at 0x7f79437b9810>) failed: KeyError: 0.0> }} sageへの入力: #pre{{ GgSave('fig_2.png', fac=0.8) }} #pre{{ Saving 11.0 x 8.0 in image. }} &ref(th_fig_2.png); *** 世帯数と分類の結果 [#lcb3ca4f] データを世帯数でソートして表示してみます。 predがグループ分けの結果で、平米地区が世帯数が中田地区よりも少ないのに グループ１に入っており、興味深い結果となっています。 sageへの入力: #pre{{ # 世帯数でソート # 平米地区に注目 df.sort('household') }} #pre{{ 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 }} *** 男性と女性の数の関係 [#te754baa] 同様に、男性と女性の関係をプロットしてみます。 世帯数よりも男性と女性の関係の方が、はっきりした直線の関係がみられます。 sageへの入力: #pre{{ # 男性（male）と女性（female）の関係をみる ggplot(df, aes(x='male', y='female', colour='pred')) + geom_point() }} #pre{{ <repr(<ggplot.ggplot.ggplot at 0x7f793f22b710>) failed: KeyError: 0.0> }} sageへの入力: #pre{{ GgSave('fig_3.png', fac=0.8) }} #pre{{ Saving 11.0 x 8.0 in image. }} &ref(th_fig_3.png); *** ペア・プロットにグループで色分け [#uc1a1dc1] 平米地区がグループ１に入っている理由をみるために、 ペア・プロットにグループで色分けしてみます。 sageへの入力: #pre{{ # グループ分けの結果で色分けしてペア・プロットを表示 d2 = df[['male', 'female', 'population', 'household', 'pred']] PandaDf2RDf(d2, "d2") }} sageへの入力: #pre{{ # 分布を可視化 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) }} &ref(th_pairs2.png); ** クラス分けを地図に表示 [#x95d57fc] *** 国土地理院の小学校区のshapeファイル [#l323446a] クラス分けの結果を地図にプロットするために、国土地理院が提供している富山県の小学校区のデータ（A27-10_16_GML）を使用しました。 - [[国土交通省 GIS データのダウンロード>http://nrb-www.mlit.go.jp/kokjo/inspect/landclassification/download/index.html]] 小学校区は、A27-10_16-g_SchoolDistrict.shpに定義され、小学校名はA27-10_16-g_SchoolDistrict.dbfに定義されています。 国土地理院のデータはShift-JISコードであり、これをUTF-8に変換する必要があります。 以下のサイトを参考にUTF-8に変換しました。 - [[Google Earth の使い方～その２ 道路や鉄道、将来人口情報などの表示>http://diy-ilands.com/2015/02/23/google-earth-usage-read-gis/]] *** 地区と小学校の関係 [#ae43b3b4] 高岡市の地区と小学校区の関係をExcelシートに整理しました。市町村合併や小学校の統合によって、 分からなかったので、高岡市の防災マップを元に人に聞きながら作ったので、間違ってかも知れませんが、 このExcelファイルを修正すれば、地図はすぐに書き直せます。 sageへの入力: #pre{{ # 高岡市地区の小学校区を読み込む schoolDistrict = pd.read_excel(DATA + "schoolDistrict.xls") # 予測結果に対する色と小学校をd2にセット schoolDistrict.head() }} #pre{{ code schoolDsitrict 0 01 平米地区 平米小学校 1 02 定塚地区 定塚小学校 2 03 下関地区 下関小学校 3 04 博労地区 博労小学校 4 05 横田地区 横田小学校 }} つぎにcodeをキーにして、高岡市のクラス分け結果（df）と小学校区（schoolDistrict）をマージします。 sageへの入力: #pre{{ # 小学校区と地区をマージ df.set_index('code') schoolDistrict.set_index('code') df = df.merge(schoolDistrict) df.head() }} #pre{{ 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への入力: #pre{{ # 小学校区のクラス分け schoolClass = df[['schoolDsitrict', 'pred']]. drop_duplicates() schoolClass.head() }} #pre{{ schoolDsitrict pred 0 平米小学校 1 1 定塚小学校 1 2 下関小学校 1 3 博労小学校 1 4 横田小学校 0 }} sageへの入力: #pre{{ # 日本語を含むため、csvファイルを使ってRへ渡す schoolClass_file = DATA + 'schoolCalss.csv' schoolClass.to_csv(schoolClass_file) }} *** maptoolsで地図にプロット [#q6aab298] Rのmaptoolsをロードし、富山県の小学校区のshapeファイルA27-10_16-g_SchoolDistrict.shpをshapeに読み込みます。 shape変数から高岡市立の小学校区のみを抽出するため、shape[shape@data$A27_006=="高岡市立",]の条件抽出を使います。 sageへの入力: #pre{{ # maptoolsをロード # r("install.packages('maptools')") r("library(maptools)") }} #pre{{ [1] "maptools" "sp" "jsonlite" "ggplot2" "stats" "graphics" "grDevices" "utils" [9] "datasets" "methods" "base" }} sageへの入力: #pre{{ # 富山県の小学校区の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') }} #pre{{ [1] 成美小学校 戸出東部小学校 戸出東部小学校 博労小学校 戸出西部小学校 伏木小学校 [7] 万葉小学校 福岡小学校 太田小学校 国吉小学校 牧野小学校 千鳥丘小学校 [13] 西広谷小学校 古府小学校 能町小学校 西条小学校 野村小学校 川原小学校 [19] 石堤小学校 定塚小学校 横田小学校 平米小学校 東五位小学校 木津小学校 [25] 南条小学校 二塚小学校 中田小学校 下関小学校 下関小学校 成美小学校 122 Levels: あさひ野小学校 さみさと小学校 井口小学校 井波小学校 宇奈月小学校 ... 立山北部小学校 }} 高岡市の小学校区のクラス分け結果のCSVファイルを読み込み、takaoka@data$A27_007に含まれる小学校のクラス結果predをセットします。 色分けを0=赤、１=青とし、Rのインデックスが１はじまりなので、pred+1としてtakaoka@data$A27_007に含まれる小学校の色をcolにセットします。 sageへの入力: #pre{{ # 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]') }} #pre{{ [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への入力: #pre{{ graph = preGraph("schoolDistrict.pdf") r('plot(takaoka, col=col)') postGraph(graph, fac=0.8) }} &ref(th_schoolDistrict.png); ** データ替えて同じ計算をする [#f3df4804] 最初のデータの読み込みを以下のように6月のデータに替えて、ActionメニューからEvaluate Allを選択すると、 平成27年6月地区別世帯数及び人口集計表に対して同じ計算を行い、その結果をグラフに表示します。 #pre{{ 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」を実行すれば、 別のデータに同じような分析を実行することができるのです。 ** コメント [#ve498b7a] #vote(おもしろかった,そうでもない,わかりずらい) 皆様のご意見、ご希望をお待ちしております。 #comment_kcaptcha
