cucumber flesh

Rを中心としたデータ分析・統計解析らへんの話題をしていくだけ

hex mapの決定版になりそうなhexmapr

mapping sf ggplot2

以前、このような記事を書きました。

uribo.hatenablog.com

ここに書いた通り、私は簡単に六角形の地図を描画できるパッケージを探しています。

今回紹介するのは、前回のhexamapmakerとは別の方法で六角形を作るhexmaprパッケージです。geospatial polygonsを正方形や六角形に変換してくれます。

hexmaprこのパッケージはCRANには登録されていません。利用の際はGitHubからインストールします。

devtools::install_github("sassalley/hexmapr")

使い方

今回は神奈川県を例に、hexmaprを使った可視化をやってみましょう。hexmaprを使った作図の方法をまとめると次のような手順になります。

  • read_polygons() などの関数を使い、対象のgeom情報をもったSpatialPolygonDataFrameを作成
  • 作図する領域(グリッドの数)などをget_shape_details()で算出する
  • calculate_cell_size()で与えたポリゴンデータとget_shape_details()の情報を元に各ポリゴンのサイズと配置を決定する
library(ggplot2)
library(viridis)
library(dplyr)
library(hexmapr)

まずは神奈川県のshapefileを国土数値情報の行政区域データ (平成28年)より取得することにします。ウェブからファイルをダウンロードしても良いのですが、ここではkokudosuuchiを利用してダウンロード先のURLを取得し、Rからファイルのダウンロード、圧縮ファイルの解凍を済ませます。

file.url <- kokudosuuchi::getKSJURL(identifier = "N03", prefCode = 14) %>% 
  arrange(desc(year)) %>%
  slice(1L) %>% 
  use_series(zipFileUrl)

download.file(file.url,
              destfile = basename(file.url))
unzip(basename(file.url))

original_shapes <- read_polygons("N03-17_14_170101.shp")

original_details <- get_shape_details(original_shapes)
original_details$nhex
# [1] 1292

new_cells_hex <- calculate_cell_size(original_shapes, 
                                     original_details,
                                     learning_rate = 0.03, 
                                     grid_type = "hexagonal", 
                                    seed = 72)

plot(new_cells_hex[[2]])

f:id:u_ribo:20171108213320p:plain

calculate_cell_size()は引数seedをもっています。これはポリゴンの算出に乱数を利用しているためで、乱数の値によりピリゴンの配置が異なるようになっています。そのため開発者は乱数を変えて、望ましい図を得ることを推奨しています。次のように良い感じの出力が得られる乱数を特定しておくと良いでしょう。

par(mfrow = c(2, 3), mar = c(0, 0, 2, 0))
for (i in 1:6){
    new_cells <-  calculate_cell_size(original_shapes, 
                                      original_details,0.03, 
                                      "hexagonal", 
                                      i)
    plot(new_cells[[2]], main = paste("Seed", i))
}

f:id:u_ribo:20171108213423p:plain

上の例では、個人的にseed 5の結果を推したいところです。

sfオブジェクトから作成する場合

先ほど用いたhexmaprのread_polygons()はsfパッケージの読み込み関数をラップしているのですが、出力されるオブジェクトはSpatialPolygonsDataFrameオブジェクトです。sfオブジェクトからhexmaprを利用するには、この形式に変換する必要があります。

以下は既存のsfオブジェクトを利用したいという時のコードです。といっても、asを使って変換するだけです。

sf.pref14 <- st_read("N03-17_14_170101.shp",
             options = c("ENCODING=CP932"))

original_shapes <- sf.pref14 %>% 
  as("Spatial")

残りのcalculate_cell_size()以下は先ほどと同じです。今度はhexではなくregular gridで作図してみましょう。

new_cells <-  calculate_cell_size(original_shapes,
                                  original_details,
                                  0.03, 
                                      "regular", 
                                      5)

plot(new_cells[[2]])

f:id:u_ribo:20171108213545p:plain

一つの市区町村で一つのポリゴンに限定する

これまでの神奈川県を例にしたhexmaprによる出力では、市区町村のポリゴンごとに一つの六角形あるいは四角形を生成する、というものでした。なので、飛び地や島を含んだ市区町村は複数行にまたがって記録されるということになっています。

POLYGONをMULTIPOLYGONに変換することで対応します(これがベストな方法なのかは不安なところ)。

sf.pref14 %<>% 
  mutate(city_name = paste(
    if_else(is.na(N03_003), 
            "",
            paste(N03_003)),
    N03_004
  ))

sf.pref14 %>% 
  # 市区町村ごとにポリゴンを結合する
  # (複数のPOLYGONをMULTIPOLYGONにする)
  group_by(city_name) %>% 
  do(out = st_union(.) %>% st_buffer(dist = 0.0001)) -> df.tmp

original_shapes <- df.tmp$out %>% 
  purrr::reduce(c) %>% 
  st_sf()

original_shapes$city_name <- as.character(df.tmp$city_name)

original_shapes %<>% 
  mutate(area = st_area(.))

original_shapes %<>% as("Spatial")
original_details <- get_shape_details(original_shapes)

new_cells_hex <- calculate_cell_size(original_shapes, 
                                     original_details,
                                     learning_rate = 0.03, 
                                     grid_type = "hexagonal", 
                                    seed = 5)

plot(new_cells_hex[[2]])

f:id:u_ribo:20171108213612p:plain

ggplot2による可視化

# 可視化のために利用
# 別な方法がありそうな
shp_clean <- function(shape){
  shape@data$id = rownames(shape@data)
  shape.points = ggplot2::fortify(shape, region = "id")
  shape.df = plyr::join(shape.points, shape@data, by = "id")
}

result_df_hex <- assign_polygons(original_shapes, new_cells_hex) %>% 
  shp_clean()

# quartzFonts(YuGo = quartzFont(rep("IPAexGothic", 4)))
# theme_set(theme_gray(base_size = 12, base_family = "IPAexGothic"))
ggplot(result_df_hex) +
  # 市区町村の面積に応じて色付けを行う
  geom_polygon(aes(x = long, y = lat, fill = as.numeric(area), group = group)) +
  geom_text(aes(V1, V2, label = city_name), size = 0.8, 
            color = "white", family = "IPAexGothic") +
  coord_equal() +
  scale_fill_viridis() +
  guides(fill = FALSE)

f:id:u_ribo:20171108213626p:plain

なんだかんだで結構しんどいのでした...

Enjoy!