Rを使って地図を表示することがしばしばある。その度にやり方を忘れていて、調べ直すことが多い。Google先生に尋ねると、当然のごとく結果が変わっている。時間が経つと、より評価され・より良い内容の記事が検索のトップに並ぶからだ。なので毎回方法が変わってしまう。自分なりの手法を確立しておきたい。
Rを用いた地理情報データとデータマッピングの流れ
まずは先行して書かれている既存のブログ記事などを参照に、Rを使って地図を描画する手順を整理してみる。大まかな流れとして、
- 必要なファイル（Shapefile）を読み込む
- 描画用のパッケージ関数を使って地図を描く
- 目的に応じてデータを地図に追加する（マッピング）
となっている。順を追って手を動かしてみよう。なお個人的に、地図をプロットする際にはこうしたい、という願望があるので、
-
Shapefileは使いたくない。GeoJSONかTopoJSONファイルを使用したい
- Shapefileの仕様が謎。複数のファイルが必要だったりして複雑
- JSONベースのファイルは軽量かつ読みやすい
プロットには
{ggplot2}を使用したい複雑ではなく、限りなくシンプルなコードにしたい
これらをできるだけ達成できるよう頑張る。
地図描画のためのパッケージあれこれ
適宜インストール。
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" ""
が用意されている。このうち重要になってくるのは、描画するソースを定義するdsnとlayerである。最低限、これだけを指定すれば読み込んでくれる。今回は{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は５つのスロットを持つ。@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()
# 国土交通省 国土地理院 ２０万分１地勢図図を元にしたファイルを用意する
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の場合、readOGRのlayer引数に渡す値が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))
{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")
ggplotのdata引数及びgeom_mapのmapに与えるデータフレームは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
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)
塗り分け地図（コロプレス図）
塗り分けは少し複雑になる。ベースとしている地図での塗り分けの基準と、データフレームでの基準を一致させていないと、塗り分けされる箇所が変わってしまう。
plot_map +
geom_map(data = df_jp,
map = map,
aes(fill = population, map_id = id, stat = "identity")) +
scale_fill_viridis()
こうした塗り分け地図を描画する専門のパッケージもある。
おまけ: 拡大・縮小できる地図
実際に動作しているのはこちらを参考にしてほしい。
RPubs - Choropleth Map of Japan with svgPanZoom, viridis and lineworkmaps packges.
library(svgPanZoom)
library(SVGAnnotation)
svgPanZoom((plot_map),
controlIconsEnabled = TRUE)
次回は応用編とさまざまな地図描画手法の紹介です
参考
- Rの基本グラフィックス機能またはggplot2を使って地図を描くには - verum ipsum factum
- [R] Rで学ぶ都知事選のデータ可視化【地理データ編】 - ill-identified diary
- choroplethrで大阪市のコロプレス図を描く - Technically, technophobic.
- plotting polygon shapefiles · hadley/ggplot2 Wiki
- Overcoming D3 Cartographic Envy With R + ggplot | rud.is
- Charting/Mapping the Scottish Vote with R (an rvest/dplyr/tidyr/TopoJSON/ggplot tutorial) | rud.is
- Mapping in R using the ggplot2 package | Technical Tidbits From Spatial Analysis & Data Science
- Mapping
- RPubs - TopoJSON map of Japan