ggplot2 で沖縄をずらして日本地図を描きたい
ill-identified
最終更新日: 2021/06/02 作成日: 2021/5/23
Rpubs は検索がしづらいので質問や指摘は以下のほうに書いていただけると反応しやすいです.
https://ill-identified.hatenablog.com/entry/2021/05/30/180927
1 はじめに
全国の天気予報図で降水確率や最高気温を表示する時, 見やすさのため離れた沖縄だけ移動して表示することはよくある. ggplot2 でこういった地図を描画する場合, 地図データのほうを手動で書き換える必要がある. アメリカ合衆国の地図でも本土から離れたアラスカやハワイを移動させて表示させる例はよくあり, 例えば fiftystater というパッケージですでに関数が用意されているが, 日本地図で同様のことをする関数はまだない1ため自作する必要がある.
twitter で指摘を受けたように, kuniezu パッケージにも同様な機能が用意されています (紹介しなかったのは単に見落としていただけです). ただし, 現時点ではこのパッケージはJGD2011の地図データを想定している (日本国内で使用されているものであれば多くはJGD2011だと思います) ので厳密に言うと想定している条件が異なります.
タイトルには ggplot2 と書いたが, 実際には plot() でも同様に (部分的に) 可能である.
なお, 以前書いた『[小ネタ] R でコロプレス図 (色分け地図) をなるべく簡単に描く』の内容が理解できる程度にRに習熟していることを前提としている.
2 大まかな技術的解説
sf, ggplto2 パッケージを使用し, 地図データを編集してグラフにする. sf パッケージは地図の情報を座標で表現されたポリゴンデータで扱うため, 平行移動, 拡大縮小, 回転, せん断, といったアフィン変換でできる範囲ならどんな変換もできる2. よって, 面積が小さい沖縄の位置を動かすだけでなく同時に拡大することもできる.
地図データは geom_sf() だけで表示することができるが, 地図の一部を編集したことを表すためのアウトライン (天気予報なんかでよく見るようなあれ) をつけたいとも思うだろう. 今回は ggplot2 を使った地図の作成なので, ggplot2 が用意する関数で様々な補助線を引くことができる. 今回は annotate() 関数を使用して線分を引く. これは見た目上は geom_segment() を使っているのと同じだが, アウトラインはデータではないため annotate() を使うことでコードを読みやすくしている.
3 実際のコード
まずは最低限必要なパッケージを読み込む (後で jpndistrict も読み込む). ggthemes は地図の表示に適したテーマ関数が用意されているので「最低限必要」なものとした.
require(NipponMap)
require(tidyverse)
require(sf)
require(ggthemes)再現可能性のため, 自分のセッション情報を掲載しておく.
sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] LC_CTYPE=ja_JP.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=ja_JP.UTF-8 LC_COLLATE=ja_JP.UTF-8
## [5] LC_MONETARY=ja_JP.UTF-8 LC_MESSAGES=ja_JP.UTF-8
## [7] LC_PAPER=ja_JP.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=ja_JP.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggthemes_4.2.4 sf_0.9-8 forcats_0.5.1 stringr_1.4.0
## [5] dplyr_1.0.5 purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
## [9] tibble_3.1.2 ggplot2_3.3.3 tidyverse_1.3.1 NipponMap_0.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 lubridate_1.7.10 class_7.3-19 assertthat_0.2.1
## [5] digest_0.6.27 utf8_1.2.1 R6_2.5.0 cellranger_1.1.0
## [9] backports_1.2.1 reprex_2.0.0 evaluate_0.14 e1071_1.7-6
## [13] httr_1.4.2 pillar_1.6.0 rlang_0.4.11 readxl_1.3.1
## [17] rstudioapi_0.13 jquerylib_0.1.3 rmarkdown_2.7 textshaping_0.3.3
## [21] munsell_0.5.0 proxy_0.4-25 broom_0.7.6 compiler_4.0.5
## [25] modelr_0.1.8 xfun_0.22 pkgconfig_2.0.3 systemfonts_1.0.1
## [29] htmltools_0.5.1.1 tidyselect_1.1.0 bookdown_0.22 fansi_0.4.2
## [33] withr_2.4.2 crayon_1.4.1 dbplyr_2.1.1 grid_4.0.5
## [37] jsonlite_1.7.2 gtable_0.3.0 lifecycle_1.0.0 DBI_1.1.1
## [41] magrittr_2.0.1 units_0.7-1 scales_1.1.1 KernSmooth_2.23-20
## [45] cli_2.4.0 stringi_1.5.3 fs_1.5.0 xml2_1.3.2
## [49] bslib_0.2.4 ellipsis_0.3.2 ragg_1.1.2 generics_0.1.0
## [53] vctrs_0.3.8 tools_4.0.5 glue_1.4.2 hms_1.0.0
## [57] yaml_2.2.1 colorspace_2.0-0 rvest_1.0.0 classInt_0.4-3
## [61] knitr_1.33 haven_2.4.0 sass_0.3.1地図に適したテーマ関数として theme_void() を設定する3. 背後のグリッドが気にならないならこの操作は無視してもよい. テーマ関数はグラフ個別に指定することも可能だが, 煩雑になるのでまとめて theme_set() で指定する.
theme_set(theme_void())ただし, Mac や Windows ユーザーで, かつ日本語を表示させたい場合はさらなる下準備4をした上で base_family = の指定が必要になる.
地図データを取得する
map_nipponmap <-
read_sf(system.file("shapes/jpn.shp", package = "NipponMap")[1],
crs = "WGS84")必要なデータと結合する. ここでは jiscode で結合できるデータがあると仮定している.
df <- left_join(
map_nipponmap,
tibble(jiscode = formatC(1:47, width = 2, format = "d", flag = "0"), x = 1:47),
by = "jiscode")沖縄を左上に移動して表示し, また, 沖縄のみ大きさを3倍に拡大してみよう.
この「沖縄だけ移動」と「アウトラインを描画」の処理をそれぞれ関数にする.
(公開後, twitter での意見を参考に R コードの可読性のため微修正)
shift_okinawa <-
function(data,
col_pref = "name",
pref_value = "Okinawa",
geometry = "geometry",
zoom_rate = 3,
pos = c(4.5, 17.5)) {
row_okinawa <- data[[col_pref]] == pref_value
geo <- data[[geometry]][row_okinawa]
cent <- sf::st_centroid(geo)
geo2 <- (geo - cent) * zoom_rate + cent + pos
data[[geometry]][row_okinawa] <- geo2
return(sf::st_as_sf(data))
}
layer_autoline_okinawa <- function(
x = c(129, 132.5, 138),
xend = c(132.5, 138, 138),
y = c(40, 40, 42),
yend = c(40, 42, 46),
size = ggplot2::.pt / 15
){
ggplot2::annotate("segment",
x = x,
xend = xend,
y = y,
yend = yend,
size = .pt / 15
)
}shift_okinawa() はデフォルトでは NipponMap パッケージから取得したデータの入力を前提としているが, 後述する jpndistrict などでは列名が異なったり, ローマ字表記か漢字表記かなどで違いが発生する. そういった場合は pref_name に都道府県を識別する列名を, pref_value に県名を与えれば, 同様に機能する.
さらに, shift_okinawa() を実行すると,
Warning in st_centroid.sfc(geometry): st_centroid does not give correctのような警告が表示されるかもしれないが, コロプレス地図で座標の正確さは要求されないと思うので無視してよい.
3.1 補足: 座標調整の内部処理について
これは, 拡大のためにまず沖縄の座標を一旦ゼロが中心になるよう平行移動してから 3 を掛け, それから座標の相対位置を元に戻し, さらに追加でずらすことで実現している.
\[ ((x, y) - (\bar{x}, \bar{y})) \times 3 + (\bar{x}, \bar{y}) + (a, b) \]
layer_autoline_okinawa() 内部で使われている annotate() は, 名前通りデータの描画ではなく補足情報を与えるための関数である. これは実質的に, アウトラインのデータを作り geom_segment() を呼び出して与えるのと同じである. \((x, y)\) から \((\mathit{xend}, \mathit{yend})\) の2点間を結ぶ線分を描画する関数で, x, xend, y, yend にベクトルを与えれば複数の線分を引くこともできる.
それっぽく見えるような座標を設定したつもりだが, もし気に入らないなら引数でも座標を入力できるようにしてある.
3.2 実行例
これらの関数を使ってみる. 色分けに上記の地図データに最初から用意されている population 列を使用した.
df %>% shift_okinawa() %>%
ggplot() +
layer_autoline_okinawa() +
geom_sf(aes(fill = population), size = .pt / 10) +
# 以下はデザインの調整
scale_fill_continuous_tableau(labels = function(x) x / 1e6) +
labs(fill = "人口 (100万人)") +
theme(legend.position = "right")
上記の関数を使いたくないなら, 例えば以下のように書くことができる. tidyverse でワンライナーで書いたほうが試行錯誤しながらの作業に向いていると思うので, そのようなコードの例にしてみる.
df %>% mutate(
geometry = if_else(
name == "Okinawa",
st_cast((geometry - st_centroid(geometry)) * 3.0 + st_centroid(geometry) + c(4.5, 17.5), "GEOMETRY"),
st_cast(geometry, "GEOMETRY")
)
) %>%
st_as_sf() %>%
ggplot() +
geom_sf(aes(fill = population), size = .pt / 10) +
annotate("segment",
x = c(129, 132.5, 138),
xend = c(132.5, 138, 138),
y = c(40, 40, 42),
yend = c(40, 42, 46),
size = .pt / 15
) +
scale_fill_continuous_tableau(labels = function(x) x / 1e6) +
labs(fill = "人口 (100万人)") +
theme(legend.position = "right")
st_cast(..., "GEOMETRY") というのは, どうも現状 sf のクラスのキャストがうまくいかないのが原因でエラーが発生するようなので, 手動で stc_geometry に変換するために必要である. もしかするともう少しシンプルな方法があるかもしれないし, 今後の更新で改善されるかもしれない.
3.3 補足: plot() を使う場合
沖縄の座標を移動する処理はグラフィックデバイスとは関係ないので同じように使える. なぜかグラフのオーバーレイはことごとく無効化されるのでアウトラインを重ねる方法はよくわからない. どうしても plot() でやりたい場合, NipponMap 本来の方法を使うべきだろう (参考: https://oku.edu.mie-u.ac.jp/~okumura/stat/nippon.html).
plot(map_nipponmap %>% shift_okinawa() %>% .["population"])
4 jpndistrict を使う場合
色分け地図の作成には, ここまで高精度な地図でないとだめという場面はあまりない5と思うが, 一応. なお, 都道府県より小さな範囲で細かい地図が必要, あるいは行政区画だけでなく河川や等高線が必要といったというときは, 国交省の国土数値情報ダウンロードサービスを利用するという手もある.
以降は上記と同じグラフを作るため, グラフ描画部分を関数化する.
plot_japan <- function(data){
ggplot(data) +
geom_sf(aes(fill = x), size = .pt / 20) +
layer_autoline_okinawa()
}jpndistrict を読み込む. 最近 CRAN から除外されてしまったので以下のページの指示に従いインストールする.
https://github.com/uribo/jpndistrict
library(jpndistrict)地図データの読み込み. NipponMap とは列名が異なることに注意.
map_jpndistrict <-map_df(1:47, ~ jpn_pref(pref_code = ., district = F)) %>%
st_simplify(preserveTopology = F, dTolerance = .01)“st_simplify does not correctly simplify …” というメッセージが表示されるが無視してよい. st_simplify(...) は意図的に地図の解像度を劣化させる処理であり, このままでは以降の描画に時間がかかりすぎるため挿入している. もし無劣化の地図データが必要ならばこの処理は省くべきである. また劣化の程度を表す 0.01 という値はおぼろげながら浮かんだ数字にすぎず, 場合によっては実際の出力を確認して調整する必要がある.
ここでは prefecture 列 (漢字表記の都道府県名) で結合できるデータがあると仮定する.
df_jpndistrict <- left_join(
map_jpndistrict,
tibble(prefecture = map_jpndistrict$prefecture, x = 1:47),
by = "prefecture"
)沖縄県を移動して表示する. 今度は拡大なし.
df_jpndistrict %>% shift_okinawa(col_pref = "prefecture", pref_value = "沖縄県", zoom_rate = 1) %>%
plot_japan() +
scale_fill_continuous_tableau()
4.1 小笠原諸島を消去する
jpndistrict のデータでは島しょ部の情報も細かい. しかし都道府県別の色分け地図ではそこまで正確である必要はない.
そこで, まずは島しょ部を除去した map_jpndistrict_crop を用意してグラフを作成する6.
以降では
although coordinates are longitude/latitude, st_intersection assumes that they are planarというメッセージが表示されるが今回は無視してよい.
map_jpndistrict_crop <- map_jpndistrict
row_tokyo <- map_jpndistrict_crop$prefecture == "東京都"
map_jpndistrict_crop$geometry[row_tokyo] <- st_crop(
map_jpndistrict$geometry[row_tokyo],
c(xmin = 0, xmax = 142, ymin = 35, ymax = 90)
)
# 他のデータと結合
df_jpndistrict_crop <- left_join(
map_jpndistrict_crop,
tibble(prefecture = map_jpndistrict$prefecture, x = 1:47),
by = "prefecture"
)
df_jpndistrict_crop %>%
shift_okinawa(col_pref = "prefecture", pref_value = "沖縄県", zoom_rate = 1) %>%
plot_japan() +
scale_fill_continuous_tableau()
上記は prefecture 列の値が "東京都" であるかどうかで東京都の地図情報を判別する前提のコードである.7 現時点での sf の仕様では, mutate(if_else(prefecture == "東京都", ..., ...)) のように書くのは難しいようである. (これも改善案を知っている人がいれば指摘歓迎)
4.2 小笠原諸島の位置をずらす
次は島しょ部を沖縄と同じようにシフトさせ, 地図を見やすくしてみる. 小笠原諸島と, 南鳥島をそれぞれ移動させる. 以下はだいぶ長くなってしまったが, 先ほどの例と同じように prefecture == "東京都" の行を取り出し, さらにポリゴンデータから島しょ部を分離して独立した行として追加したのが map_jpndistrict_shift である.
# 島しょ部をシフトさせるための地図
map_jpndistrict_shift <- map_jpndistrict
# 島しょ部を分離する
map_tokyo_islands <- filter(map_jpndistrict_shift, prefecture == "東京都")
map_tokyo_islands <- mutate(map_tokyo_islands, geometry = st_crop(
filter(map_jpndistrict_shift, prefecture == "東京都")$geometry[1],
c(xmin = 142, xmax = 180, ymin = 0, ymax = 35)
),
islands = "島しょ部"
)
# さらに南鳥島を分ける
map_tokyo_mts <- map_tokyo_islands
map_tokyo_mts <- mutate(
map_tokyo_mts,
geometry = st_crop(
map_tokyo_islands$geometry[1],
c(xmin = 153, xmax = 155, ymin = 23, ymax = 25)),
islands = "南鳥島"
)
map_tokyo_islands$geometry[1] <- st_crop(
map_tokyo_islands$geometry, c(xmin = 0, xmax = 152, ymin = 24, ymax = 90)
)
# 元の地図から島しょ部を除去
map_jpndistrict_shift$geometry[13] <- st_crop(
map_jpndistrict_shift$geometry[13],
c(xmin = 0, xmax = 142, ymin = 35, ymax = 90)
)
map_jpndistrict_shift <- bind_rows(
map_jpndistrict_shift,
map_tokyo_islands,
map_tokyo_mts
)同じように視覚化したいデータフレームを結合させて色分け地図にする. 小笠原諸島と南鳥島をそれぞれ位置調整し, それぞれ3倍と10倍に拡大し, さらにアウトラインを表示するのでこのコードも長くなっている.
df_jpndistrict_shift <- left_join(
map_jpndistrict_shift,
tibble(prefecture = map_jpndistrict$prefecture, x = 1:47),
by = "prefecture"
)
df_jpndistrict_shift %>% shift_okinawa(col_pref = "prefecture", pref_value = "沖縄県", zoom_rate = 1) %>%
# 小笠原諸島移動&拡大
mutate(geometry = if_else(
prefecture == "東京都" & islands == "島しょ部" & !is.na(islands),
(geometry - st_centroid(geometry)) * 3 + c(2.5, 6) + st_centroid(geometry),
geometry
)) %>%
# 南鳥島移動&拡大
mutate(geometry = if_else(
prefecture == "東京都" & islands == "南鳥島" & !is.na(islands),
(geometry - st_centroid(geometry)) * 10.0 + c(-9.5, 5) + st_centroid(geometry),
geometry
)) %>%
plot_japan() +
annotate("segment",
x = c(142.5, 142.5, 143.5),
xend = c(142.5, 147, 147),
y = c(27, 37, 30.7),
yend = c(37, 37, 30.7),
size = .pt / 15
) +
scale_fill_continuous_tableau()
こういうページを見つけることができたが, sf 登場以前で, なおかつ応用しづらいので新たに書いた. https://www.mk-mode.com/blog/2014/12/03/gis-japan-okinawa-division/↩︎
https://r-spatial.github.io/sf/articles/sf3.html こういった用語がわからない場合は https://imagingsolution.blog.fc2.com/blog-entry-284.html の図を参照.↩︎
この関数は ggthemes パッケージによって提供される. なお,
theme_map()という関数も存在するが, ヘルプによるとgeom_map()関数に対応したテーマである. これは sf パッケージが登場する前に使われていたもので, 今回のgeom_sf()を使った描画とは相性が悪い (と個人的には思う) ので今回は使用しないことにした.↩︎この下準備の詳細は https://speakerdeck.com/ktgrstsh/display-cjk-font-in-any-gpraphic-device-and-platform-2020, https://uribo.hatenablog.com/entry/2021/03/29/202756 などを参考にする.↩︎
理由は以前ここに書いた通り. https://ill-identified.hatenablog.com/entry/2020/12/07/134705↩︎
以下のようにこの処理は複雑なので, できれば最初から都道府県より細かい単位で分割された地図データを用意したほうが簡単である. 例えば
jpn_pref()にdistrict = Tを指定すると, より細かい行政区分でデータが出力される.↩︎