読者です 読者をやめる 読者になる 読者になる

choroplethrで大阪市のコロプレス図を描く

R choroplethr

大阪都構想住民投票、盛り上がってましたね。

関連して調べてたら、大阪市のshpファイルは以下で手に入るみたいでした。ライセンスはCC-BYです。

大阪市市政 マップナビおおさかオープンデータ一覧

せっかくなので、これをchoroplethrで使う方法を探してみました。

chroplethrとは

こんな感じの図(コロプレス図)を描くためのパッケージです。

f:id:yutannihilation:20150523173608p:plain:w300

Japan.R 2014の発表スライドとかが参考になると思います。

www.slideshare.net

データのダウンロード

別にここはRでやる必要はないですが、Rでやるとこんな感じです。unzip()は、日本語ファイル名を含むZipアーカイブに対してはうまく動かないので、仕方なくsystem()を使っています。(参考:昨日のメモRpubs

tmp_dir <- tempdir()
zip_file <- file.path(tmp_dir, "24kuikishape.zip")
zip_dir  <- file.path(tmp_dir, "24kuikishape")

download.file("http://www.city.osaka.lg.jp/contents/wdu090/opendata/mapinfo/24kuikishape.zip",
              destfile = zip_file)

system(sprintf('unzip -Ocp932 %s -d %s', zip_file, zip_dir))

list.files(zip_dir)
#> [1] "24区画像.dbf"     "24区画像.prj"     "24区画像.sbn"     "24区画像.sbx"     "24区画像.shp"    
#> [6] "24区画像.shp.xml" "24区画像.shx"

shpファイルの読み込みと整形

ここは、何をやってるのかよく分かりませんが、、https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefilesに手順があるのでそれをそのままなぞります。よく分からないので、おかしい部分あれば教えてください。

library(rgdal)
library(maptools)
library(dplyr)
library(ggplot2)

# レイヤーを調べる
ogrListLayers(file.path(zip_dir, "24区画像.shp"))
#> [1] "24区画像"
#> attr(,"driver")
#> [1] "ESRI Shapefile"
#> attr(,"nlayers")
#> [1] 1

# 「24区画像」しかないのでそれを読み込む
osaka <- readOGR(zip_dir, layer = "24区画像")

osakaは、SpatialPolygonsDataFrameという形式で読み込まれます。これを、ggplot2でプロットできる形式(=choroplethrで使える形式)に変形します。

regionという名前のカラム(ここでは区のコード)は、各区をあらわすキーになります。統計データ側にも同じくregionカラムをつくって、この地図データとマッチさせます。regionintegerでも構わないんですが、下で統計データを処理するときに使ってるplyr::mapvaluescharacterを返すので、それに合わせてcharacterにしています。(データ型も同じでないとマッチしてくれない)

# rownameは、あとでfortifyするのに使う
osaka@data    <- osaka@data %>%
  add_rownames(var = "id") %>%
  mutate(region = as.character(ATTR1),
         region_name = as.character(ATTR2))

head(osaka@data)
#>   id ATTR1  ATTR2 住所コード 町丁目名称 区コード region region_name
#> 1  0 27104 此花区          0       <NA>        0  27104      此花区
#> 2  1 27118 城東区          0       <NA>        0  27118      城東区
#> 3  2 27116 生野区          0       <NA>        0  27116      生野区
#> 4  3 27106   西区          0       <NA>        0  27106        西区
#> 5  4 27122 西成区          0       <NA>        0  27122      西成区
#> 6  5 27124 鶴見区          0       <NA>        0  27124      鶴見区

# idでfortifyする
osaka.points  <- fortify(osaka, region="id")

head(osaka.points)
#>        long       lat order  hole piece group id
#> 1 -48720.97 -145682.5     1 FALSE     1   0.1  0
#> 2 -48720.82 -145682.6     2 FALSE     1   0.1  0
#> 3 -48719.33 -145681.4     3 FALSE     1   0.1  0
#> 4 -48712.64 -145689.8     4 FALSE     1   0.1  0
#> 5 -48712.57 -145689.9     5 FALSE     1   0.1  0
#> 6 -48712.52 -145689.9     6 FALSE     1   0.1  0

# osaka@dataとosaka.pointを合わせる
osaka.df      <- osaka.points %>%
  inner_join(osaka@data, by = "id")

これで地図のデータができました。

Choroplethを継承したR6クラスをつくる

上のデータをつかってコロプレス図を描くR6のクラスをつくります。とりあえずinitialize()だけあればいいみたいですが、各区の名前を表示したい、みたいな場合はrender()を上書きする必要があるみたいです。(例: choroplethr/state.R at f5cc2971dc36c6abcaf51c41ec0ff604b9cba355 · trulia/choroplethr · GitHub)

library(R6)

OsakaChoropleth <- R6Class(
  "OsakaChoropleth",
  inherit = choroplethr:::Choropleth,
  public = list(
    initialize = function(user.df)
    {
      super$initialize(osaka.df, user.df)
    }
  )
)

コロプレス図で可視化したいデータを用意する

大阪市各区の人口構成のデータを加工してCSVにしたやつがあるので、それを使います。作り方とかはREADMEのリンク先をみてください。

github.com

csv_file <- file.path(tmp_dir, "osaka_age_composition.csv")

download.file("https://raw.githubusercontent.com/yutannihilation/osaka_age_composition/master/osaka_age_composition.csv",
              destfile = csv_file, method = "curl")

age_comp.df.raw <- readr::read_csv(csv_file)
head(age_comp.df.raw)
#> Source: local data frame [6 x 5]
#> 
#>   age total male female district
#> 1   0  1113  591    522       北
#> 2   1   969  501    468       北
#> 3   2   819  395    424       北
#> 4   3   870  488    382       北
#> 5   4   719  345    374       北
#> 6   5   657  320    337       北

これを、さっきつくったOsakaChoroplethクラスから使えるような形に変形します。regionはさっきと合わせないといけません。とりあえず、平均年齢を出してみることにします。

age_comp.df <- age_comp.df.raw %>%
  group_by(district) %>% 
  summarise(value = sum(age * total)/sum(total)) %>%
  mutate(district = paste0(district, '区')) %>%
  mutate(region = plyr::mapvalues(district,
                                  from = osaka@data$region_name,
                                  to   = osaka@data$region))

head(age_comp.df)
#> Source: local data frame [6 x 3]
#> 
#>   district    value region
#> 1   中央区 42.05251  27128
#> 2 住之江区 46.55447  27125
#> 3   住吉区 45.90159  27120
#> 4     北区 43.32589  27127
#> 5   城東区 44.68260  27118
#> 6   大正区 47.32154  27108

これをつかってさっきつくったクラスをnew()してrender()すればおわりです。

c  <- OsakaChoropleth$new(age_comp.df)
c$render() + coord_equal()

f:id:yutannihilation:20150523200846p:plain:w450

感想

choroplethrはvignetteがしっかりしてて親切なんですが、自分の地図をつくるときの説明はちょっと分量少ない気がします。さくっと使えるようになるにはもうちょい修業がいりそう。

参考にしたサイト