注記 2014/11/8 シェイプファイルの利用元を ESRI から国土数値情報に変更し, 若干修正
概要
筆者は空間情報データにも空間統計学にも素人なので注意してください.
導入
前回に引き続き, 都知事選のデータを題材にRでのデータを視覚化する実例を紹介する. 今回は, 地理空間情報を加工した地図データを重ね合わせることに挑戦してみる. 地図データは株式会社 ESRIジャパンが無償公開している日本国内の市区町村界データ (http://www.esrij.com/products/data/japan-shp/) を使用する. 使用条件第2条に「ArcGIS以外のソフトで使用するな」とあるが, ShapeFileライブラリによると非営利目的の範囲なら R でもデータを加工編集してもいいらしいのでありがたく使わせていただく*1.
GIS データの操作については, Rによる空間データの統計分析, 統計データの視覚化を参考にした.
シェイプファイルの読み込み
ダウンロードしたデータには .shp, shx, .dbf, ... といくつもファイルがあるが, すべて必要で, 同じディレクトリに入れておく. 以降これらをまとめてシェイプファイルと呼ぶ.
# R 3.1.2 library(maptools) # ver. 0.8-30 library(ggplot2) # ver. 1.0.0 tokyo.map <- readShapeSpatial("GIS/mlit/tokyo/N03-13_13_130401.shp") # パスは各自任意で変更すること. # シェイプファイルに含まれる属性データは以下の通り. # N03_001: 都道府県名 # N03_002: 支庁・振興局名 (北海道のみ) # N03_003: 郡・政令市名 # N03_004: 市区町村名 # N03_007: 行政区域コード # fortify で data.frame に変換した後の結合キーにするため, rowname を変数に代入 # 参考: http://www.r-bloggers.com/using-r-working-with-geospatial-data-and-ggplot2/ tokyo.map$id <- rownames(tokyo.map@data) # 市区町村ごとの選挙情報データとの結合のため市区町村の変数を作成する: 区名だけ N03_004に書かれているので, # 市区町村名が欠損なら群・政令市名を市区町村名とする. tokyo.map$SIKUCHOSON <- ifelse(is.na(tokyo.map$N03_004),yes = as.character(tokyo.map$N03_003),no=as.character(tokyo.map$N03_004)) # 不要になった変数削除 tokyo.map$N03_001 <- NULL; tokyo.map$N03_002 <- NULL; tokyo.map$N03_003 <- NULL tokyo.map$N03_004 <- NULL; tokyo.map$N03_007 <- NULL
読み込まれたmap オブジェクトには, 市区町村界の形状を記録したポリゴンデータや, 対応する都道府県・市区町村名など含まれている. そのため, 都道府県名・市区町村名をたどることで, 東京都以外のデータを切り落とすことができる.
一方, 逆に結合する場合はやや面倒になる. 今回は結合はしないので, 統計データの視覚化で紹介されている市区町村界データを結合することで市区町村界を消した福岡県のデータを作成する方法などを参照されたい.
ただし, シェイプファイルはシフトJIS (cp932) なので, ウィンドウズ以外では文字化けする可能性がある. そこで, 自分は 上記の読み込みを行う前に ShapeFileライブラリの全国市区町村界ShapeFileの県別分解プログラム(新）のコードを流用して文字コードを変換した. そのため, 元データではローマ字になっている map$SIKUCHOSON も漢字表記に変わっていることに注意する.
このデータを地図としてプロットするには, そのまま plot(tokyo.map) でできる. これをもとにコロプレス地図 (色分け地図) を作成したり, その他の付随情報を追加編集する場合 たとえば classInt パッケージや spdep パッケージが必要となる. この辺は先ほども挙げたRによる空間データの統計分析, 統計データの視覚化などが詳しい*2. しかし, これらのパッケージで提供される関数は若干複雑であり, またすでに ggplot2 の解説をしてきたこともあり, 今回も意地でも ggplot2 を用いることとする.
コロプレス図の描画
ggplot2 でデータをプロットするには, データフレーム型でデータを渡さなければならない. しかし, ポリゴンデータはデータフレーム型ではないため, 変換が必要となる. Maps with ggplot2 ではまさにその方法が提供されているが, 今回自分が試したところ, 一部の境界線が正常に出力されなかった*3. 加えてポリゴンデータと属性情報のみの変換であり, 各市区町村のラベルを置く座標を記録した labpt を取得できないことから, 今回は上記のサイトの関数は使わなかった.
実は ggplot2 にある fortify() 関数でポリゴンデータの変換が可能である*4. 以下のコードでは, さらに選挙に関する情報 df をシェイプファイルの属性データ tokyo.map@data に結合させ, fortify() で tokyo.map@polygons からポリゴンデータを取り出したものに, 市区町村ごとの属性情報を再度結合し, グラフ作成用データを作成した. df は東京選挙管理委員会等の公表データをもとに, 筆者が編集したデータフレームで,
area vote.ratio retired.ratio mean.income age masuzoe hosokawa 1 あきる野市 40.09986 59.90014 548.1566 45.36686 19.64500 7.401318 2 稲城市 47.97954 52.02046 600.8897 41.74498 21.61101 9.541792 3 羽村市 38.84584 61.15416 546.1129 43.47425 18.55585 7.156338 4 奥多摩町 44.37660 55.62340 NA 57.36676 26.27536 6.815048 5 葛飾区 40.84311 59.15689 514.6153 45.01322 19.30639 6.711528 6 江戸川区 39.34888 60.65112 522.6092 42.35804 19.10781 6.343519
というふうに市区町村ごとの数値を記録してある. tokyo.map$SIKUCHOSON と df$area に市区町村名が入っているので, この2変数で2つのデータをひもづけることができる.
この後, fortify() を使って ggplot() が読み込めるデータフレームに変換するのだが, 現状では fortify() 関数はシェイプファイルのポリゴンデータ (つまり市区町村の形状のデータ) 部分しか変換してくれない. そこで, データフレーム化したポリゴンデータに, 再度属性データを結合することになる. 今度の結合は, 市区町村名ではなく, 予め作成しておいた id を使う.
# 属性データに選挙情報を結合 # 市区町村名をキーとして左結合 tokyo.map <- merge(x=tokyo.map, y=df, by.x="SIKUCHOSON",by.y="area", all.x=T) # 島嶼部除外 islands <- c("大島町", "利島村", "新島村", "神津島村", "三宅村", "御蔵島村", "八丈町", "青ヶ島村", "小笠原村", "所属未定地") tokyo.map <- tokyo.map[!tokyo.map$SIKUCHOSON %in% islands,] # ポリゴンデータをデータフレームに変換 tokyo.gg.df <- fortify(tokyo.map) # データフレーム化したポリゴンデータに属性データを左結合 tokyo.gg.df <- merge(x=tokyo.gg.df, y=tokyo.map@data, by="id", all.x = T)
以上で必要なデータを作成できた. しかしその前にあともうひとつ, グラフを見やすくするため, 市区町村名を表示させるためのデータも作成しておこう. シェイプファイルには labpt という変数が保持されており, これが市区町村ごとに表示させるラベルの表示位置座標として使用できる. そこで, 下記のコードで座標と市区町村名を持つデータフレームを作成した.
# 区名ラベル位置座標 location <- lapply(tokyo.map@polygons, function(x) slot(x,"labpt") ) location <- as.data.frame(matrix(unlist(location), ncol=2, byrow=T, dimnames=list(NULL,c("x","y"))) ) location$id <-unlist(lapply(tokyo.map@polygons, function(x) slot(x,"ID") ) ) # ラベル位置と市区町村名結合 label.gg.df <- merge(x=location, y=unique(tokyo.gg.df[,c("SIKUCHOSON","id")]), by.x="id", by.y="id", all.x = T) label.gg.df <- label.gg.df[!duplicated(label.gg.df$SIKUCHOSON),] # 重複削除 (雑)
ただし, このやり方だと, 飛び地がある場合に表示位置がずれる可能性がある. 今回はあまり影響がなかったのでこの方法で通したが, これでずれる場合は最後の1行を工夫する必要がある. 手動で直すという手も無いわけでもない.
ここまですれば後は今までと本質的に同じである. なお, ここでは ggsave() 関数を用いて画像を自動で保存している. この関数は保存時のサイズをインチ, センチ, ミリメートルでしか指定できない関数である.
注記 2014/03/01 市区町村の境界線がなく, 見づらかったので修正. 以降すべての geom_polygon() 関数に color="white" を追加しました
# 脚注の文章作成 note <- data.frame(x=139, y=35.55, text1="行政区域は国土数値情報より利用\n 平均年齢は住民基本台帳による数字\n http://www.toukei.metro.tokyo.jp/juukiy/2013/jy13000001.htm より2013年1月時点のもの", text2="行政区域は国土数値情報より利用\n 世帯年収は2010年住宅・土地統計調査確報より推計.\n データえっせい『東京の年収地図』\n(http://tmaita77.blogspot.jp/2012/12/blog-post.html) を参考に", text3="行政区域は国土数値情報より利用\n 投票率・棄権率は東京都選挙管理委員会の発表\n (http://www.senkyo.metro.tokyo.jp/h26chijisokuho/index.html) より") credit1 <- geom_text(aes(x,y, label=text1, hjust=0),data=note, size=3.5) credit2 <- geom_text(aes(x,y, label=text2, hjust=0),data=note, size=3.5) credit3 <- geom_text(aes(x,y, label=text3, hjust=0),data=note, size=3.5) g.theme <-theme(axis.title=element_blank(), axis.text=element_blank(), axis.ticks=element_blank()) ggplot(tokyo.gg.df) + geom_polygon(aes(long,lat, group=group, fill=age), color="white") + labs(title="市区町村別平均年齢") + g.theme + scale_fill_continuous(name="平均年齢") + geom_text(data =label.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + credit1 + coord_equal() ggsave(filename = "tokyo_avg_age.png", width = 10.24, height = 6.4, unit ="in" ) # 幅はインチ, cm, mm でしか指定できない? ggplot(tokyo.gg.df) + geom_polygon(aes(long,lat, group=group, fill=mean.income), color="white") + labs(title="市区町村別世帯年収分布") + g.theme + scale_fill_continuous(name="平均世帯年収\n(万円)") + geom_text(data =label.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + credit2 + coord_equal() ggsave(filename = "tokyo_avg_income.png", width = 10.24, height = 6.4, unit ="in" ) ggplot(tokyo.gg.df) + geom_polygon(aes(long,lat, group=group, fill=vote.ratio), color="white") + labs(title="2014年都知事選\n市区町村別投票率") + g.theme + scale_fill_continuous(name="投票率\n(%)") + geom_text(data =label.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + credit3 + coord_equal() ggsave(filename = "tokyo_vote.png", width = 10.24, height = 6.4, unit ="in" ) ggplot(tokyo.gg.df) + geom_polygon(aes(long,lat, group=group, fill=retired.ratio), color="white") + labs(title="2014年度都知事選\n市区町村別棄権率") + g.theme + scale_fill_continuous(name="棄権率\n(%)") + geom_text(data =label.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + credit3 + coord_equal() ggsave(filename = "tokyo_retire.png", width = 10.24, height = 6.4, unit ="in" )
ポリゴンデータの描画には, geom_polygon() 関数を使うことになる. ただし, 境界内の塗りつぶしが不要なら, geom_path() でも良い. group=group とすることでで閉じた図形を定義し, それぞれの fill に色分けしたい属性情報を与えている. なお, group=id とすると, 地区町村で一つの図形と認識されるため, お台場など, 飛び地となっている図形が正常に表示されないので注意*5. 座標軸は緯度・経度で与えられているが, 今回は表示する意味がないため, 非表示にしている. この設定は先に g.theme() オブジェクトに与えることで使いまわしている. また, credit* はデータの出処を明記するためにgeom_text() を入れてある*6ので, 公表するのでなければ省いて構わない. また, coord_equal() 関数は, 縦横比を固定する関数である.
この方法でできた図が以下である. 基本的に色が薄いほど数値が高い. ただし年収の灰色は欠損値である. 中央区・千代田区・港区の平均年収の高さと投票率の高さが印象的になる.
23区に限定した地図
この図では, 面積の広い西部の市町村が強調され*7, 23区の様子がわかりづらい. そこで, 23区以外を切り落として再度表示させてみる. 地図の切り落としは, データフレーム型でも基本的なやり方は変わらない.
label.23.gg.df <- label.gg.df[grep(pattern = "区",x = label.gg.df$SIKUCHOSON),]
今回は市区町村名の漢字表記データも残しておいたので, grep() で名前に "区" を含むものだけを抜き出すことで, 23区のデータのみ抜き出している.
note2 <- note note2[,c("x","y")] <- c(139.55, 35.52) credit1 <- geom_text(aes(x,y, label=text1, hjust=0, vjust=0),data=note2, size=3) credit2 <- geom_text(aes(x,y, label=text2, hjust=0, vjust=0),data=note2, size=3) credit3 <- geom_text(aes(x,y, label=text3, hjust=0, vjust=0),data=note2, size=3) ggplot(tokyo23.gg.df) + geom_polygon(aes(long,lat, group=group, fill=age), color="white") + geom_text(data =label.23.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + scale_fill_continuous(name="平均年齢") + labs(title="市区町村別平均年齢") + g.theme + credit1 + coord_equal() ggsave(filename = "tokyo23_avg_age.png", width = 10.24, height = 6.4, unit ="in" ) ggplot(tokyo23.gg.df) + geom_polygon(aes(long,lat, group=group, fill=mean.income), color="white") + geom_text(data =label.23.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + scale_fill_continuous(name="平均世帯年収\n(万円)") + labs(title="市区町村別世帯年収分布") + g.theme + credit2 + coord_equal() ggsave(filename = "tokyo23_avg_income.png", width = 10.24, height = 6.4, unit ="in" ) ggplot(tokyo23.gg.df) + geom_polygon(aes(long,lat, group=group, fill=vote.ratio), color="white") + geom_text(data =label.23.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + scale_fill_continuous(name="投票率\n(%)") + labs(title="2014年都知事選\n市区町村別投票率") + g.theme + credit3 + coord_equal() ggsave(filename = "tokyo23_vote.png", width = 10.24, height = 6.4, unit ="in" ) ggplot(tokyo23.gg.df) + geom_polygon(aes(long,lat, group=group, fill=retired.ratio), color="white") + geom_text(data =label.23.gg.df, aes(x,y,label=SIKUCHOSON), size=3) + scale_fill_continuous(name="棄権率\n(%)") + labs(title="2014年度都知事選\n市区町村別棄権率") + g.theme + credit3 + coord_equal() ggsave(filename = "tokyo23_retire.png", width = 10.24, height = 6.4, unit ="in" )
23区に拡大すると, 棄権率と平均収入を交互に見比べると色が反転しているように見える. 収入の高い区ほど投票率が高い (=年収の低い区ほど投票率が低い) 傾向がよく分かり, 興味深い.
補足: 平均世帯年収の計算方法
平均世帯年収の計算は, 基本的には東京の年収地図と同じである. 不明な世帯を除外し, 1500万円以上の階級値を2000万円に設定したところも同じである. それ以外の階級値は中央値としたので, 数値は若干異なるかもしれないが, 大きな影響はでないと思う.
2017/10/21追記: 実際にはこのような方法では想像より誤差がありそうである. 以下の記事でこのような概算で発生する問題について書いた. ill-identified.hatenablog.com
参考
- stack overflow rgdal package installation
- Rjp Wiki ShapeFileライブラリ
- 古谷知之 (2011) 『Rによる空間データの統計分析』, 朝倉書店, ISBN: 978-4-254-12815-4
- 谷村晋 (2010)『地理空間データ分析』, 金明哲編「Rで学ぶデータサイエンス」7巻, 共立出版, ISBN: 978-4-320-01927-0
山本義郎・飯塚誠也・藤野友和 (2013)『統計データの視覚化』, 金明哲編「Rで学ぶデータサイエンス」12巻, 共立出版, ISBN: 978-4-320-11016-8 (10章のみ)
R-bloggers (2010) 『Maps with ggplot2』
- 舞田俊彦 (2012) 『東京の年収地図』, 「データえっせい」, 2012年12月2日
*1:その他にも無償提供されたデータは存在するが, 必ずしもデータ構造が同じとは限らないため, 以降のコードで再現できない可能性がある.
*2:前者には1箇所だけミスがあり, readShape... 関数をsdped パッケージとしている (正確にはmaptools) ので注意.
*3:先にも言ったように素人なの修正するのが面倒臭い
*4:Maps with ggplot2 の著者がこれを使わなかったのは, 著者の64bit mac os でうまく行かなかったかららしい.
*5:実は R-bloggers の方法ではこのへんがうまく行かなかったのだが, 修正が面倒だった
*6:これが面倒なのだが, 脚注を付ける関数がないらしいので苦労している. なにかいい方法はないだろうか?
