盆暗の学習記録

データサイエンスを中心として,日々学んだことの備忘録としていく予定です。初心者であり独学なので内容には誤りが含まれる可能性が大いにあります。

Rで小地域データの塗り分け地図を描く

小地域(町字レベルの細かい地域区分のデータ)のシェープファイル(.shp)を使用してコロプレス図(塗り分け地図)を描く方法のメモ。

f:id:nigimitama:20180925013333p:plain

1. データの取得

まず,シェープファイルを取得します。

e-statに行き,「統計GIS」から

f:id:nigimitama:20180924215229p:plain

「境界データダウンロード」をクリックしていった先からダウンロードできます。

f:id:nigimitama:20180924215324p:plain

今回は「小地域」の「国勢調査」のデータを使ってみます。

f:id:nigimitama:20180924220136p:plain

f:id:nigimitama:20180924220155p:plain

データの説明書である「定義書」のpdfをダウンロードし,「小地域(町丁・字等別)」をクリックします。

f:id:nigimitama:20180924220324p:plain

シェープファイルが欲しいので「Shape形式」のものを選ぶのですが,地図の作成方式の違いから「世界測地系緯度経度・Shape形式」と「世界測地系平面直角座標系・Shape形式」があります。

世界測地系緯度経度・Shape形式」を選びます。

f:id:nigimitama:20180924220505p:plain

港区のデータを選んでみます。

f:id:nigimitama:20180925002532p:plain

世界測地系緯度経度・Shape形式」をクリックすれば,シェープファイルが入ったzipファイルがダウンロードされます。

2. データの読み込みとデータの確認

{rgdal}パッケージを使用して読み込んでみます。

library(rgdal)
shape <- readOGR("h27ka13103.shp",
                 stringsAsFactors = FALSE, encoding = "UTF-8")

RStudioのEnvironment部で,読み込んだシェープファイルの中身を確認してみます。

f:id:nigimitama:20180924223947p:plain

f:id:nigimitama:20180925011036p:plain

このpolygonsには地図の図形が入っています。

そしてdataというデータフレームには,自治体コードや自治体名,町字,面積といった地域のデータと,人口,世帯数といった国勢調査のデータが入っています。

> head(shape@data)
     KEY_CODE PREF CITY S_AREA PREF_NAME CITY_NAME     S_NAME KIGO_E HCODE      AREA
0 13103001001   13  103 001001    東京都      港区   芝1丁目   <NA>  8101  68807.76
1 13103001002   13  103 001002    東京都      港区   芝2丁目   <NA>  8101 165462.02
2 13103001003   13  103 001003    東京都      港区   芝3丁目   <NA>  8101 153048.48
3 13103001004   13  103 001004    東京都      港区   芝4丁目   <NA>  8101  94327.10
4 13103001005   13  103 001005    東京都      港区   芝5丁目   <NA>  8101 215053.82
5 13103002001   13  103 002001    東京都      港区 海岸1丁目   <NA>  8101 391922.17
(以下省略)

塗分ける前の地図データを確認してみます。

{leaflet}で描画します。

# シェープファイルの単純表示
library(leaflet)
shape %>% 
leaflet() %>% 
  addTiles() %>% 
  setView(lat = 35.65, lng = 139.75, zoom = 12) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addPolygons(fillOpacity = 0.5,
              weight = 1,
              fillColor = "lightblue")

こんな感じです。

f:id:nigimitama:20180925005129p:plain

3. 塗り分け地図の作図

地図作成の準備

まず,塗り分け用に人口密度データを作ります。(人口で塗り分けると面積の差が出て自治体同士を正しく比較できないため)

定義書によれば,shape@data$JINKOが人口データで,shape@data$AREAが面積(㎡)のようなので,これを使います。

# 人口密度データの作成
PopDensity <- as.numeric(shape@data$JINKO) / shape@data$AREA * 1000000 # 1平方キロあたりにするために*1000000

つぎに,塗り分け地図の装飾用のセッティングを行います。

plot用のカラーパレットを定義して…

# 塗る色(連続値のカラーパレット)をセット
pal <- colorNumeric("Blues", domain = PopDensity, reverse=F)

マウスオーバー時のラベルを定義します。

# マウスオーバー時の表示内容を設定
labels <- sprintf("<strong>%s</strong><br/>%5.1f",
                  paste0(shape@data$MOJI),
                  PopDensity) %>% lapply(htmltools::HTML)

ここのラベルの定義は

> head(labels)
[[1]]
<strong>芝1丁目</strong><br/>22366.7

こういうhtmlコードにしたいからやりました。

そして,leafletでプロットします

地図のプロット

# 地図にプロット
shape %>% 
  leaflet() %>% 
  # setView() : 地図をズームした状態で表示する
  setView(lat = 35.65, lng = 139.75, zoom = 12) %>% 
  # addProviderTiles() : 背景のタイルを指定
  addProviderTiles(providers$CartoDB.Positron) %>% 
  # addPolygons() : 塗り分け地図の描画
  addPolygons(fillOpacity = 0.5,
              weight=1,
              color = "#666",
              fillColor = ~pal(PopDensity), # ここで人口密度データを使う
              label = labels,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "15px",
                direction = "auto"),
              highlight = highlightOptions(
                weight = 5,
                color = "#666",
                dashArray = "",
                fillOpacity = 0.7,
                bringToFront = TRUE)) %>% 
  # addLegend() : 凡例の設定
  addLegend("bottomright", pal = pal, values = ~PopDensity, # ここで人口密度データを使う
            title = "人口密度(人口/km2)")

f:id:nigimitama:20180925010851p:plain

結果はこんな感じ。マウスオーバーで詳細を表示してくれます。

Windowsでそのまま作るとMSゴシックのギザギザとしたフォントがなんだか醜いですね(フォント厨並みの感想)

オプションでfontなどを指定していきましょう

地図の装飾

# 地図にプロット
map <- shape %>% leaflet() %>% 
  # setView() : 地図をズームした状態で表示する
  setView(lat = 35.65, lng = 139.75, zoom = 13) %>% 
  # addProviderTiles() : 背景のタイルを指定
  addProviderTiles(providers$CartoDB.Positron) %>% 
  # addPolygons() : 塗り分け地図の描画
  addPolygons(fillOpacity = 0.5,
              weight=1,
              color = "#666",
              fillColor = ~pal(PopDensity), # ここで人口密度データを使う
              label = labels,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal",
                             "font-family" = "Yu Gothic",
                             "box-shadow" = "3px 3px rgba(0,0,0,0.25)",
                             "padding" = "3px 8px"),
                textsize = "15px",
                direction = "auto"),
              highlight = highlightOptions(
                weight = 5,
                color = "#666",
                dashArray = "",
                fillOpacity = 0.7,
                bringToFront = TRUE)) %>% 
  # addLegend() : 凡例の設定
  addLegend("bottomright", pal = pal, values = ~PopDensity, # ここで人口密度データを使う
            title = "人口密度(人口/km2)")

ラベルの文字に関する装飾はlabelOptions = labelOptions(style = list())CSSを指定することができます。

凡例(legend)は関数としてOptionを受け入れていないので,leafletの地図オブジェクトを作った後{htmltools}CSSを追加します。

# legendも游ゴシックにする
library(htmltools)
browsable(
  tagList(list(
    tags$head(
      tags$style(
        ".leaflet .legend {
            font-family : Yu Gothic;
         }
        "
      )
    ),
    map
  ))
)

f:id:nigimitama:20180925013333p:plain

見た目がよくなりました。