Posted at

【まとめ用】ggplot2を使って地図を書く

More than 3 years have passed since last update.

Rを使って地図を表示することがしばしばある。その度にやり方を忘れていて、調べ直すことが多い。Google先生に尋ねると、当然のごとく結果が変わっている。時間が経つと、より評価され・より良い内容の記事が検索のトップに並ぶからだ。なので毎回方法が変わってしまう。自分なりの手法を確立しておきたい。


Rを用いた地理情報データとデータマッピングの流れ

まずは先行して書かれている既存のブログ記事などを参照に、Rを使って地図を描画する手順を整理してみる。大まかな流れとして、


  1. 必要なファイル(Shapefile)を読み込む

  2. 描画用のパッケージ関数を使って地図を描く

  3. 目的に応じてデータを地図に追加する(マッピング)

となっている。順を追って手を動かしてみよう。なお個人的に、地図をプロットする際にはこうしたい、という願望があるので、



  1. Shapefileは使いたくない。GeoJSONかTopoJSONファイルを使用したい


    • Shapefileの仕様が謎。複数のファイルが必要だったりして複雑

    • JSONベースのファイルは軽量かつ読みやすい



  2. プロットには{ggplot2}を使用したい


  3. 複雑ではなく、限りなくシンプルなコードにしたい


これらをできるだけ達成できるよう頑張る。


地図描画のためのパッケージあれこれ

適宜インストール。

library(magrittr)

library(sp)
library(ggplot2)
library(rgdal)
library(broom)
library(ggmap)
library(viridis)


1 必要なファイルを用意し、Rで扱えるようにする

幾つかの方法があるが、ここでは{rgdal}readOGR関数を利用してファイルを読み込む。readOGRでは、ベクター形式のファイルを読み込むことができる。readOGRの引数をみると、

args(readOGR) %>% as.list() %>%  names()

##  [1] "dsn"                       "layer"                    

## [3] "verbose" "p4s"
## [5] "stringsAsFactors" "drop_unsupported_fields"
## [7] "input_field_name_encoding" "pointDropZ"
## [9] "dropNULLGeometries" "useC"
## [11] "disambiguateFIDs" "addCommentsToPolygons"
## [13] "encoding" "use_iconv"
## [15] "swapAxisOrder" "require_geomType"
## [17] "integer64" ""

が用意されている。このうち重要になってくるのは、描画するソースを定義するdsnlayerである。最低限、これだけを指定すれば読み込んでくれる。今回は{rgdal}に内臓されているデータを使うことにする。次のコマンドを実行して出力されるディレクトリにデモ用のShapefileがある。

path <- vector()

system.file("vectors", package = "rgdal") %>% {
path <<- .
print(.)
}

## [1] "/Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/vectors"

次のようなディレクトリ構造になっている(一部のみ出力)。

list.files(path) %>% grep(pattern = "^trin_inca_pl03", x = ., value = TRUE)

## [1] "trin_inca_pl03.dbf" "trin_inca_pl03.shp" "trin_inca_pl03.shx"

trin_inca_pl03.shpに付属して、幾つかのファイルが用意されている。Shapefileを出力させる場合、これらのファイルも必要になる。では、このtrin_inca_pl03ファイルをRに取り込んでみる。

map <- data.frame()

readOGR(dsn = path, layer = "trin_inca_pl03") %>% {
map <<- .
print(class(.))
str(., max.level = 2)
}

## OGR data source with driver: ESRI Shapefile 

## Source: "/Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/vectors", layer: "trin_inca_pl03"
## with 3 features
## It has 12 fields
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 2 obs. of 12 variables:
## ..@ polygons :List of 2
## ..@ plotOrder : int [1:2] 1 2
## ..@ bbox : num [1:2, 1:2] -87.1 5.5 -79 10.1
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

readOGRで読み込んだオブジェクトは{sp}パッケージで定義されるSpatialPointsDataFrameというものになる。S4メソッドで定義されているので、各オブジェクトへのアクセスは@スロットを利用する。strで示したように、SpatialPointsDataFrameは5つのスロットを持つ。@dataは属性値を持つデータフレームであり、@coordsは位置座標の行列オブジェクトである。

dplyr::tbl_df(map@data)

## Source: local data frame [2 x 12]

##
## ENGL_NAME SCI_NAME AUTHORITY FAMILY PRESENCE
## (fctr) (fctr) (fctr) (fctr) (int)
## 1 Wandering Tattler Tringa incana (Gmelin, 1789) Scolopacidae 1
## 2 Wandering Tattler Tringa incana (Gmelin, 1789) Scolopacidae 1
## Variables not shown: ORIGIN (int), COMPILER (fctr), SCALE (fctr), TAX_COM
## (fctr), DIST_COM (fctr), REFERENCES (fctr), REVIEWERS (fctr)

地図を描画する場合、元になるファイルが必要になるが、多くの場合Shapefileを使うことになるだろう。しかし、readOGRではShapefile以外のJSON(GeoJSON, TopoJSON)形式のファイルでも読み込むことが可能である。いずれの方法もShapefileの時と同じである。

# topojson

map <- data.frame()
# 国土交通省 国土地理院 20万分1地勢図図を元にしたファイルを用意する
file.path <- "https://gist.github.com/uribo/b09d642351c03dd975aa/raw/6cd2f7e70922b24f89b442c96c001dc3ab794317/japan_one_twenty_map.topojson"
readOGR(dsn = file.path,
layer = "OGRGeoJSON",
stringsAsFactors = FALSE) %>% {
map <<- .
print(class(.))
str(., max.level = 2)
}

## OGR data source with driver: GeoJSON 

## Source: "https://gist.github.com/uribo/b09d642351c03dd975aa/raw/6cd2f7e70922b24f89b442c96c001dc3ab794317/japan_one_twenty_map.topojson", layer: "OGRGeoJSON"
## with 130 features
## It has 2 fields
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 130 obs. of 2 variables:
## ..@ polygons :List of 130
## .. .. [list output truncated]
## ..@ plotOrder : int [1:130] 75 78 59 41 79 93 8 11 61 31 ...
## ..@ bbox : num [1:2, 1:2] 123 24 149 46
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

GeoJSONあるいはTopoJSONの場合、readOGRlayer引数に渡す値がOGRGeoJSONとなっているが、これは

ogrListLayers(dsn = file.path)[1]

## [1] "OGRGeoJSON"

で出力されるものである。出力したいlayers名が不明な場合、こちらで確認すると良いだろう。

対象のファイルがGeoJSONあるいはTopoJSON形式である場合、{rgdal}の代わりに{geojsonio}を利用することも可能である。

geojsonio::topojson_read(file.path) %>% {

print(class(.))
str(., max.level = 2)
}


2 地図をプロットする

readOGRにより、SpatialPolygonsDataFrameオブジェクトを作成したら、ひとまず地図がプロットできるようになる。plot(map)のように、plot関数にオブジェクト名を指定するだけで地図が描画される。ggplotによるプロットは次のようにすれば実行されるが、これはポリゴンとして描画しているため解像度が荒い。

## # ggplot2::ggplot による描画: よくない例。荒い。data.frame, SpatialPolygonsDataFrameどちらでもおk

ggplot() + geom_polygon(data = map,
aes(x = long, y = lat, group = id))

unnamed-chunk-11-1.png

{ggplot2}を使って描画させる場合は、SpatialPolygonsDataFrameをdata.frameに変換する必要がある。その際、{broom}tidyを利用すると便利である(これまではggplot2::fortifyが使われてきたが、fortify は将来廃止予定とのことであり、 代わりに{broom}パッケージを使うことをHadleyは推奨している)。

broom::tidy(map) %>% {

map <<- .
print(class(.))
dplyr::glimpse(.)
}

## [1] "data.frame"

## Observations: 650
## Variables: 7
## $ long (dbl) 148, 149, 149, 148, 148, 148, 149, 149, 148, 148, 147, 1...
## $ lat (dbl) 46.00000, 46.00000, 45.33333, 45.33333, 46.00000, 45.333...
## $ order (int) 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4,...
## $ hole (lgl) FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, ...
## $ piece (fctr) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ group (fctr) 1.1, 1.1, 1.1, 1.1, 1.1, 2.1, 2.1, 2.1, 2.1, 2.1, 3.1, ...
## $ id (chr) "1", "1", "1", "1", "1", "2", "2", "2", "2", "2", "3", "...

tidy処理を施したオブジェクトは共通の列名を持つデータフレームとなる。このオブジェクトを用いてggplotの上に地図をプロットしていく。

ggplot() + 

geom_map(data = map,
map = map,
aes(x = long, y = lat, map_id = id),
fill = "white", color = "black")

ggplotdata引数及びgeom_mapmapに与えるデータフレームはtidy処理を施したオブジェクトである。審美的属性に与える要素も共通である。これが基本型となるが、これだと地図が塗りつぶしされていたり、プロットが正方形になっていたりと、地図として格好が悪いので、次のように幾つかの属性を付加しておくと良い。

ggplot() +

geom_map(data = map,
map = map,
aes(x = long, y = lat, map_id = id),
color = "black", fill = "white", # 境界線を黒く、領域を白で塗りつぶす
size = 0.5) + 
coord_map(projection = "mercator") -> plot_map # 投影法を指定(メルカトル図法)
# ggthemes::theme_map() も地図描画用のggplotテーマなので、指定しておくと良いかもしれない。
# 背景が白くなり、軸などが表示されなくなる

plot_map

gg_map_basic-1.png

plot_mapとして、再帰的に地図を描画できるようにした。ここから必要なデータを重ねたりしていく。


3 マッピング

いよいよ本番。自分のデータを地図上に追加していく。まず、データを用意するところから始め、マッピングを進めていく。地図上にデータを重ねていくマッピングの手法として、ラベルの出力、ポイントの表示、値や水準による塗り分けがあると考えられるので、それぞれについて日本の地図を例を元にして説明する。


データの用意

マッピングを行う際に必要になるのは、どの座標(緯度と経度)にどのような値(文字や種類)を出力するか、塗り分けの基準となる値や水準といったものなので、これらを含んだデータフレームを用意しておく。

df_jp <- data.frame(address          = c("Tokyo, Japan", 

"Okayama, Japan",
"Sapporo, Japan",
"Okinawa, Japan"),
id = c(56, 90, 26, 125),
population = c(8945695, 709584, 1913545, 315954),
# https://ja.wikipedia.org/wiki/都道府県庁所在地と政令指定都市の人口順位 から人口を入力
stringsAsFactors = FALSE) %>%
mutate_geocode(address) %>% # addressを元に、緯度経度を求める
dplyr::mutate(address = gsub(pattern = "[[:punct:]] Japan$", replacement = "", address))

dplyr::tbl_df(df_jp)

## Source: local data frame [4 x 5]

##
## address id population lon lat
## (chr) (dbl) (dbl) (dbl) (dbl)
## 1 Tokyo 56 8945695 139.6917 35.68950
## 2 Okayama 90 709584 133.9350 34.66167
## 3 Sapporo 26 1913545 141.3469 43.06417
## 4 Okinawa 125 315954 127.8014 26.33583


ラベル及びポイントの表示

緯度・経度、出力する文字列はそれぞれ、df_jp$lat, df_jp$lon, df_jp(addressである。また、ポイントの大きさを人口df_jp)population`を元にして表示してみる。

plot_map + 

# ポイントの出力
geom_point(data = df_jp,
aes(x = lon, y = lat),
fill = "skyblue",
shape = 21, alpha = 0.8, size = 2 * log(df_jp$population) / 10) +
geom_text(data = df_jp,
aes(label = address, x = lon + 0.2, y = lat + -0.4),
size = 3,
color = "tomato",
angle = 20)

unnamed-chunk-15-1.png


塗り分け地図(コロプレス図)

塗り分けは少し複雑になる。ベースとしている地図での塗り分けの基準と、データフレームでの基準を一致させていないと、塗り分けされる箇所が変わってしまう。

plot_map + 

geom_map(data = df_jp,
map = map,
aes(fill = population, map_id = id, stat = "identity")) +
scale_fill_viridis()

unnamed-chunk-16-1.png

こうした塗り分け地図を描画する専門のパッケージもある。

https://github.com/trulia/choroplethr


おまけ: 拡大・縮小できる地図

実際に動作しているのはこちらを参考にしてほしい。

RPubs - Choropleth Map of Japan with svgPanZoom, viridis and lineworkmaps packges.

library(svgPanZoom)

library(SVGAnnotation)

svgPanZoom((plot_map), 

controlIconsEnabled = TRUE)

次回は応用編とさまざまな地図描画手法の紹介です


参考