Rで日本地図を描く方法
変更履歴・概要
2023/6/13 初版
2023年6月に行った作業のメモをもとに整理したものです。それぞれの環境によって違いがあるかもしれませんが，ご容赦ください。
できるようになること
Rで日本地図を描けるようになります。WEB地図画像を（背景として）利用する方法と，ベクタ境界データをダウンロードして利用する方法を紹介します。また，Rのパッケージに組み込まれた境界データも紹介します。
WEB地図の利用（leaflet）
leafletはWEB地図を作成するためのJavaScriptのオープンソースライブラリです。leafletパッケージはleaflet.jsをRでも使えるようにしたもので，JavaScriptを使わなくてもRだけで利用できます。このパッケージを利用して，日本地図の画像をRで表示してみます。
leafletパッケージの詳細は，以下のページを参考にしてください。
https://rstudio.github.io/leaflet/
https://kazutan.github.io/kazutanR/leaflet_d.html
HTMLウィジェットとして出力されるので，rmarkdownでの出力形式をHTMLとすると，動的な地図を作成できます。が，Markdown形式への出力では何も作成されません。
動的な地図はあきらめて，静的な地図画像を出力してみます。実はknitr（v1.13以降）ではHTML以外への出力時にHTMLウィジェットを出力させようとすると，webshotパッケージにより静的なスクリーンショットを出力します（https://bookdown.org/yihui/bookdown/html-widgets.html）。
webshotパッケージとPhantomJSは以下によりインストールできます。
install.packages("webshot")
webshot::install_phantomjs() # PhantomJSをインストール
webshot::is_phantomjs_installed() # インストールされているか確認
デフォルトではOpenStreatMapが表示されます。
library(leaflet)
m <- leaflet() %>%
addTiles() %>%
setView(135, 35, zoom=5)
m
本記事のテーマは日本地図なので，地理院タイルの白地図を利用してみます。
attr <- "<a href='http://maps.gsi.go.jp/development/ichiran.html' target='_blank'>地理院タイル</a>"
m <- leaflet() %>%
addTiles("https://cyberjapandata.gsi.go.jp/xyz/blank/{z}/{x}/{y}.png", attribution = attr) %>%
setView(135, 35, zoom=5)
m
なお，leafletパッケージの出力として静的な地図画像を出力したい場合は，mapviewパッケージのmapshot()により出力できます。
mapview::mapshot(m, file="leaflet.png")
ベクタ境界データの利用（tmap，ggplot2）
都道府県のベクタ境界データとして公的機関から入手できるものとしては，国土数値情報の行政区域データがありますが，全国で427MBと重いです。また，総務省統計局が整備している政府統計の統合窓口から国勢調査等の小地域データをダウンロードできますが，こちらは都道府県単位でしかダウンロードできません。
ここでは，国土地理院が作成した地球地図日本の行政界データを利用してみます。QGISによる作業説明を参考にさせていただきました。
必要なパッケージ
ベクタ境界データはsfパッケージで読み込みます。Geocomputation with Rを一読されることをおすすめします。必要なパッケージは以下の通りです。
library(sf)
library(tidyverse)
library(tmap)
地図データの読み込み
地球地図日本ダウンロードから第2.1版ベクタ（2015年公開）行政界のShapeファイルをダウンロードして展開します。展開したフォルダにあるpolbnda_jpn.shpをread_sf()で読みこみます（gm_jpn.pathにファイルのパスを指定してください）。
gm_jpn <- read_sf(gm_jpn.path)
class(gm_jpn)
## [1] "sf" "tbl_df" "tbl" "data.frame"
head(gm_jpn)
## Simple feature collection with 6 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 140.6929 ymin: 41.71093 xmax: 144.4865 ymax: 43.957
## Geodetic CRS: ITRF94
## # A tibble: 6 × 10
## f_code coc nam laa pop ypc adm_code salb soc
## <chr> <chr> <chr> <chr> <int> <int> <chr> <chr> <chr>
## 1 FA001 JPN Hokkai Do Sapporo Shi 1.93e6 2014 01100 UNK JPN
## 2 FA001 JPN Hokkai Do Hakodate Shi 2.74e5 2014 01202 UNK JPN
## 3 FA001 JPN Hokkai Do Otaru Shi 1.27e5 2014 01203 UNK JPN
## 4 FA001 JPN Hokkai Do Asahikawa Shi 3.49e5 2014 01204 UNK JPN
## 5 FA001 JPN Hokkai Do Muroran Shi 9.13e4 2014 01205 UNK JPN
## 6 FA001 JPN Hokkai Do Kushiro Shi 1.80e5 2014 01206 UNK JPN
## # ℹ 1 more variable: geometry <POLYGON [°]>
市区町村レベルの境界データなので，都道府県レベルに集約します。都道府県名はnamにあるので，ディソルブ（dissolve）します。
gm_jpn <- read_sf(gm_jpn.path) %>%
group_by(nam) %>% summarize()
head(gm_jpn)
## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 132.0119 ymin: 32.89807 xmax: 141.6831 ymax: 41.5562
## Geodetic CRS: ITRF94
## # A tibble: 6 × 2
## nam geometry
## <chr> <MULTIPOLYGON [°]>
## 1 Aichi Ken (((137.1258 34.6534, 137.1252 34.65187, 137.1254 34.64987,…
## 2 Akita Ken (((140.6984 39.46827, 140.7065 39.4742, 140.72 39.47933, 1…
## 3 Aomori Ken (((140.8723 41.01786, 140.8718 41.01773, 140.8718 41.0162,…
## 4 Chiba Ken (((139.9716 35.671, 139.9787 35.66987, 139.9792 35.66987, …
## 5 Ehime Ken (((133.1697 34.29093, 133.1697 34.29033, 133.1722 34.2868,…
## 6 Fukui Ken (((135.8062 35.6322, 135.806 35.6356, 135.8059 35.63927, 1…
地図表示
これで都道府県レベルに集約できたので，地図として表示してみます。
base Rのplot()
plot(gm_jpn)
tmapパッケージ
tm_shape(gm_jpn, projection = "EPSG:3857") + tm_polygons()
ggplot2パッケージ
ggplot() + geom_sf(data=gm_jpn) +
coord_sf(crs=sf::st_crs("EPSG:3857"),default_crs = sf::st_crs("EPSG:4326")) +
theme_bw()
なお，tmapとggplot2では多くの世界地図で用いられるメルカトル投影の座標系に再投影しています。他の座標系でも表示できます。日本でよく利用されるものはhttps://lemulus.me/column/epsg-list-gisにわかりやすく整理されています。
（参考1）一部の都道府県だけ着色する
pref_choice <- c('Hyogo Ken', 'Iwate Ken', 'Okinawa Ken', 'Tokyo To', "Hokkai Do", "Kyoto Fu")
gm_jpn %>%
mutate(pref=if_else(condition = nam %in% pref_choice, true = T, false = F)) %>%
ggplot() + geom_sf(aes(fill = pref), color="black", show.legend = FALSE) +
coord_sf(crs=sf::st_crs("EPSG:3857"), default_crs = sf::st_crs("EPSG:4326")) +
scale_fill_manual(values = c("white", "red")) +
theme_bw()
（参考2）地球地図日本データにある人口データ（2014年）を表示させる
read_sf(gm_jpn.path) %>%
# -99999999（不明）と-89999999（メインポリゴン付与）のポリゴンは0とする
mutate(Pop=if_else(condition = pop < 0, true = 0, false = pop)) |>
group_by(nam) %>%
summarize(Pop=sum(Pop)) %>%
tm_shape(projection = "EPSG:3857") +
tm_fill(col="Pop") + tm_borders(col="black")
沖縄県だけをずらして描画
次に見やすさのため，沖縄県だけをずらして描画してみます。
最初に（しなくてもよいのですが）地図表記を簡略にするため10 sfパッケージ（v.1.0.0以降）では Google のS2球面幾何エンジンを利用できるので，地理座標系（経緯度）データのままで面積計算が可能です。デフォルトではオンですが，S2球面幾何エンジンがオンになっていることを確認しておくほうがよいと思います。
# S2球面幾何エンジンがオンかを確認
sf_use_s2()
## [1] TRUE
# ディゾルブ後のデータをpolygonに分割，属性の割り当てで警告がでる
gm_jpn_OnlyLgIs <- gm_jpn %>% st_cast("MULTIPOLYGON") %>%
st_cast("POLYGON")
## Warning in st_cast.sf(., "POLYGON"): repeating attributes for all
## sub-geometries for which they may not be constant
# ポリゴン面積を計算し，10km2より大きなもののみを選択
gm_jpn_OnlyLgIs <- mutate(gm_jpn_OnlyLgIs, area=st_area(gm_jpn_OnlyLgIs)) %>%
filter(area > units::set_units(1*10^7, m^2))
まずはtmapパッケージを利用してみます。Geocomputation with Rの第9章を参考にしています。
# 沖縄県とそれ以外の都道府県のデータに分割
map2 <- gm_jpn_OnlyLgIs %>% filter(nam == "Okinawa Ken") %>%
tm_shape(projection = "EPSG:3857") +
tm_polygons()+
tm_layout(frame=F)
map1 <- gm_jpn_OnlyLgIs %>% filter(nam != "Okinawa Ken") %>%
tm_shape(projection = "EPSG:3857") +
tm_polygons()+
tm_layout(frame=F)
# 2つの地図を結合して表示。大きさや位置は適当なので，必要に応じて見直してください
map1
print(map2, vp = grid::viewport(0.25, 0.8, width = 0.45, height = 0.45))
grid::grid.lines(x = c(0.01, 0.47), y = c(0.65, 0.65), gp = grid::gpar(lty = 2))
grid::grid.lines(x = c(0.47, 0.47), y = c(0.99, 0.65), gp = grid::gpar(lty = 2))
次はggplot2パッケージでの作成例です。shift_okinawa()とlayer_autoline_okinawa()はggplot2 で沖縄をずらして日本地図を描きたいのコードを利用させていただきました。
gm_jpn_OnlyLgIs %>% shift_okinawa(col_pref="nam",
pref_value="Okinawa Ken", zoom_rate = 1) %>%
ggplot() + geom_sf() +
coord_sf(crs=sf::st_crs("EPSG:3857"), default_crs = sf::st_crs("EPSG:4326")) +
layer_autoline_okinawa() +
theme_bw() +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())
なお，kuniezuパッケージを利用すると簡単にできそうです（が，CRANからはダウンロードできません）。
その他パッケージの利用
パッケージのなかには境界データを含むものがあります。利用することで（比較的簡単に）地図を作成できます。
jpndistrictパッケージ
国土数値情報の行政区域データを元にした都道府県・市町村の境界データ（と思われます）。利用条件は国土数値情報の利用約款によるようです。現在，CRANからはダウンロードできません。詳しくは以下のページで確認してください。
https://github.com/uribo/jpndistrict
NipponMapパッケージ
都道府県単位で，コロプレス図（階級区分図）を描くためのパッケージ。都道府県境界は作者が修正するなどしており，正確ではないので面積や距離を測るには不適。
https://cran.r-project.org/web/packages/NipponMap/index.html
http://www.eeso.ges.kyoto-u.ac.jp/emm/materials/rmemo/japanprefmap
https://okumuralab.org/~okumura/stat/nippon.html
mapdataパッケージ
mapsパッケージの補足で，都道府県境界を含むデータ（legacy packageだそうです）。元データ等は不明です。
https://cran.r-project.org/web/packages/mapdata/index.html
spDataパッケージ
Natural Earthのworldデータが用意されています。Geocomputation with Rではこのパッケージのデータを使って解説されています。
https://cran.r-project.org/web/packages/spData/index.html
Session info
sessionInfo()
## R version 4.3.0 (2023-04-21 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Japanese_Japan.utf8 LC_CTYPE=Japanese_Japan.utf8
## [3] LC_MONETARY=Japanese_Japan.utf8 LC_NUMERIC=C
## [5] LC_TIME=Japanese_Japan.utf8
##
## time zone: Etc/GMT-9
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods
## [7] base
##
## other attached packages:
## [1] tmap_3.3-3 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
## [5] dplyr_1.1.2 purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
## [9] tibble_3.2.1 ggplot2_3.4.2 tidyverse_2.0.0 sf_1.0-12
## [13] leaflet_2.1.2
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 viridisLite_0.4.2
## [3] farver_2.1.1 blob_1.2.4
## [5] fastmap_1.1.1 mapview_2.11.0
## [7] XML_3.99-0.14 digest_0.6.31
## [9] timechange_0.2.0 lifecycle_1.0.3
## [11] ellipsis_0.3.2 processx_3.8.1
## [13] terra_1.7-29 magrittr_2.0.3
## [15] compiler_4.3.0 rlang_1.1.1
## [17] RODBC_1.3-20 tools_4.3.0
## [19] utf8_1.2.3 yaml_2.3.7
## [21] knitr_1.43 htmlwidgets_1.6.2
## [23] bit_4.0.5 sp_1.6-0
## [25] classInt_0.4-9 here_1.0.1
## [27] RColorBrewer_1.1-3 abind_1.4-5
## [29] KernSmooth_2.23-20 withr_2.5.0
## [31] leafsync_0.1.0 odbc_1.3.4
## [33] grid_4.3.0 stats4_4.3.0
## [35] fansi_1.0.4 e1071_1.7-13
## [37] leafem_0.2.0 colorspace_2.1-0
## [39] scales_1.2.1 dichromat_2.0-0.1
## [41] cli_3.6.1 crayon_1.5.2
## [43] rmarkdown_2.21 generics_0.1.3
## [45] rstudioapi_0.14 tzdb_0.4.0
## [47] tmaptools_3.1-1 DBI_1.1.3
## [49] proxy_0.4-27 stars_0.6-1
## [51] parallel_4.3.0 s2_1.1.4
## [53] base64enc_0.1-3 vctrs_0.6.2
## [55] webshot_0.5.4 jsonlite_1.8.4
## [57] callr_3.7.3 hms_1.1.3
## [59] bit64_4.0.5 crosstalk_1.2.0
## [61] units_0.8-2 lwgeom_0.2-13
## [63] glue_1.6.2 leaflet.providers_1.9.0
## [65] codetools_0.2-19 ps_1.7.5
## [67] stringi_1.7.12 gtable_0.3.3
## [69] raster_3.6-20 munsell_0.5.0
## [71] pillar_1.9.0 htmltools_0.5.5
## [73] satellite_1.0.4 R6_2.5.1
## [75] wk_0.7.3 microbenchmark_1.4.10
## [77] rprojroot_2.0.3 evaluate_0.21
## [79] lattice_0.21-8 highr_0.10
## [81] png_0.1-8 class_7.3-21
## [83] Rcpp_1.0.10 xfun_0.39
## [85] fs_1.6.2 pkgconfig_2.0.3
Discussion