2021/08/22追記
以前書いた、{jpndistrict}のjpn_pref()の方法では日本地図を書くためのデータを作ろうとしたときにErrorが出るようになってしまってうまく行かなかった。R 4.1.1で久しぶりに日本地図を書こうと思ったらうまく行かなかったのでおそらくR 4.1.xからかそれぐらいにアップデートされた周辺パッケージによるものかと思うが？残していたR 4.0.5だとうまく行った。
ちゃんと検証していないのでどの環境（多分R 4.1.ｘ以降）で再現するかがわからないが、{jpndistrict}のjpn_pref()で一部都道府県で市町村を省いたデータを作るためのdistrict = FALSEと引数をするとErrorが出てしまう。
library(jpndistrict)
## This package provide map data is based on the Digital Map 25000 (Map
## Image) published by Geospatial Information Authority of Japan (Approval
## No.603FY2017 information usage <https://www.gsi.go.jp>)
# 岩手県（JISコード03）を市町村なしでデータを出そうとするとErrorがでる jpn_pref(3, district = FALSE)
## Error in s2_geography_from_wkb(x, oriented = oriented, check = check): Evaluation error: Found 12 features with invalid spherical geometry.
## [2] Loop 45 is not valid: Edge 0 has duplicate vertex with edge 2
## [3] Loop 23 is not valid: Edge 1 crosses edge 3
## [6] Loop 12 is not valid: Edge 1 is degenerate (duplicate vertex)
## [9] Loop 10 is not valid: Edge 1 is degenerate (duplicate vertex)
## [10] Loop 1 is not valid: Edge 0 is degenerate (duplicate vertex)
## [24] Loop 12 is not valid: Edge 0 has duplicate vertex with edge 2
## [25] Loop 1 is not valid: Edge 1 is degenerate (duplicate vertex)
## [26] Loop 6 is not valid: Edge 1 is degenerate (duplicate vertex)
## [27] Loop 10 is not valid: Edge 1 is degenerate (duplicate vertex)
## [28] Loop 20 is not valid: Edge 0 is degenerate (duplicate vertex)
## ...and 2 more.
R（パッケージ仕様？）のジオデータの取り扱いの仕様変更に伴うもの？よくわからない。
日本地図を作るためのジオデータを国土地理院のページからダウンロードしたものを都道府県ごとのデータに整形するが必要になり、自前でデータを整形するならシェープファイルを読むの手法でやる必要がある、が結構処理に時間がかかるのでこだわりがなければシェープファイルを読むの途中に乗っている整形済みジオデータをつかうのがよいと思う（他人のふんどしで戦う）。
追記終了
Rの{ggplot2}で日本地図と日本の一部地域（東北地方や関東地方など）の地図を書きたいと思った。とりあえず地図を描くだけの簡単な方法をメモしておく。
地理データのことは正直良くわからないがとりあえず白地図を{ggplot2}で描くため調べたことをまとめておく。ベストプラクティスかどうかはわからない。
- 地図データを作る（全国版）
- 日本地図を描く
- 例：東北地方の地図を描く
- 例：関東甲信越地方の地図を描く（離島を除いた版も描く）
- 座標系をUTMに変換する
- そもそも一地方だけならデータを作るところからサクッとやれる
- 蛇足：{jpndistrict}がパッケージをロードしなければ使えないのは……
- 補足：地図データの入手
- 参考
地図データを作る（全国版）
地図を描画するミソは地図データ（座標データ）を持ってくるまたは作るところだと思う。
今回は地図データを作るのに{tidyverse}と{jpndistrict}、あと生データのままだと境界線が精緻すぎるのでデフォルメするために{rmapshaper}のms_simply()を使う1。
# `{tidyverse}`や`{sf}`がインストールされていないなら次のコードを実行してインストールする。 install.packages("tidyvese") install.packages("rmapshaper") # `{jpndistrict}`がインストールされていないなら次のコードを実行してインストールする。 # 現在`{jpndistrict}`はCRANに登録されていないのでGithubからインストールする。 # `{remotes}`がインストールされていないのなら`install.packages("remotes")`で先にインストールする。 remotes::install_github("uribo/jpndistrict")
# パッケージをロードする。 # 特に、`{jpndistrict}`はロードしていないと使えない（蛇足後述）。 library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## √ ggplot2 3.3.3 √ purrr 0.3.4
## √ tibble 3.1.0 √ dplyr 1.0.5
## √ tidyr 1.1.3 √ stringr 1.4.0
## √ readr 1.4.0 √ forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(jpndistrict)
## This package provide map data is based on the Digital Map 25000 (Map
## Image) published by Geospatial Information Authority of Japan (Approval
## No.603FY2017 information usage <https://www.gsi.go.jp>)
library(rmapshaper)
## Registered S3 method overwritten by 'geojsonlint':
## method from
## print.location dplyr
とりあえず、日本地図全体を描く。全国版の地図データを作る2。
# 地図データを`japan_map`として`{jpndistrict}`の`jpn_pref()`で作る。 # 1:47で都道府県JISコードの1北海道から47沖縄県までのデータを作るように指定。 japan_map <- 1:47 %>% # `jpn_pref()`に都道府県JISコードを渡して、`district = FALSE`で市町村データを除くようにする。 # `district = TRUE`のままだと、市町村までのデータになる。 map(~ jpn_pref(pref_code = ., district = FALSE)) %>% # 1北海道～47沖縄までのデータをくっつける。 reduce(rbind) %>% # データが精緻すぎるのでgeometryをデフォルメする。今回は東京都離島を描画したい（後述する方法で除く方法を描くため）ので約1/100にデータをデフォルメする。 # ちなみにもとから1/1000にデフォルメすると結構な離島が消える。 ms_simplify(keep = 0.01, keep_shapes = TRUE)
一応データを見てみる。
head(japan_map)
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 139.333 ymin: 37.73284 xmax: 148.8932 ymax: 45.55566
## Geodetic CRS: WGS 84
## # A tibble: 6 x 3
## pref_code prefecture geometry
## <chr> <chr> <MULTIPOLYGON [arc_degree]>
## 1 1 北海道 (((146.9328 44.60531, 146.9472 44.63215, 146.9541 44.626~
## 2 2 青森県 (((139.8578 40.60251, 139.86 40.61657, 139.8918 40.63429~
## 3 3 岩手県 (((140.7654 39.21918, 140.77 39.22447, 140.7863 39.22765~
## 4 4 宮城県 (((140.4448 38.14382, 140.4525 38.1527, 140.4666 38.1577~
## 5 5 秋田県 (((139.6906 39.99485, 139.6988 40.00803, 139.7094 40.002~
## 6 6 山形県 (((139.5386 38.55663, 139.5566 38.5761, 139.5562 38.5840~
日本地図を描く
日本地図を描くのはこのデータで{ggplot2}のgeom_sf()で地図を描くだけ。テーマを調整しないで日本地図全体を描く。
japan_map %>% ggplot() + geom_sf()
ちなみに、geometoryをデフォルメしてないと精緻すぎて描画するまで時間がかかる。
例：東北地方の地図を描く
あとは描画したい地域だけのデータを抽出すれば、任意の地域の地図だけ書ける。例えば東北地方だけなら次のように書ける。
japan_map %>% # 都道府県JISコードの`pref_code`が文字列なので数値に変換する。 mutate(pref_code = as.numeric(pref_code)) %>% # 都道府県JISコードで2青森県～7福島県を指定する。 filter(pref_code >= 2 & pref_code <= 7) %>% ggplot() + geom_sf()
都道府県JISコードが並んでいれば楽ちんである（filter()はその範囲が指定されていれば良いので色々書き方はある）。
例：関東甲信越地方の地図を描く（離島を除いた版も描く）
例えば関東甲信越地方なら、都道府県JISコードが関東地方と、甲信越（長野県、新潟県、山梨県）でちょっと離れているので、ちょっとめんどくさい。
japan_map %>% mutate(pref_code = as.numeric(pref_code)) %>% # 都道府県JISコードで8茨木県～14神奈川県を指定する。 # 残りの三県もJISコードで指定してもよいが、`prefcture`で直接県名との一致で条件一致させてみる。 filter((pref_code >= 8 & pref_code <= 14) | prefecture %in% c("長野県","新潟県","山梨県")) %>% ggplot() + geom_sf()
東京都の離島まで描出されるのでやや見にくい（前述の通り意図的に離島まで描出した）がちゃんと関東甲信越市方が描画される。
離島も東京都なので安易に描画しないようにすべきではない気がするが、図の大半が海になっているので見慣れた関東甲信越市方の地図にするために比較的小さな離島を除くために{rmapshaer}のms_filter_islands()を使う。
japan_map %>% mutate(pref_code = as.numeric(pref_code)) %>% filter((pref_code >= 8 & pref_code <= 14) | prefecture %in% c("長野県","新潟県","山梨県")) %>% # `min_area`で除外する島の最大の大きさを指定する。東京都の目立つ離島で最も大きい大島が約91km^2なので100km^2（指定はm^2）で指定しておく。 ms_filter_islands(min_area = 100000000) %>% ggplot() + geom_sf()
いい感じに見慣れた関東甲信越地方の地図になった。佐渡ヶ島も除きたければmin_areaを任意でさらに大きくするといい。
座標系をUTMに変換する
ところで、座標系をUTM（Universal Transverse
Merctor）に変換するには途中に、{sf}をつかってst_transform(data, "+proj=utm +zone=54 +datum=WGS84 +units=km")をかませるといい3。
ここまでの記事の流れで、{rmapshaper}をインストールしていると依存関係で{sf}がインストールされているはずだが、インストールしていなかったらinstall.packages("sf")でインストールする。
library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
japan_map %>% st_transform("+proj=utm +zone=54 +datum=WGS84 +units=km") %>% ggplot() + geom_sf()
そもそも一地方だけならデータを作るところからサクッとやれる
さて、そもそも東北地方だけなど一地域だけの地図を描画して今後地図を描画する予定がないなら、地図データを任意の地域だけで作るのもまた一つの方法である。
# 2:7で都道府県JISコードの2青森県から7福島県までのデータを作るように指定 touhoku_map <- 2:7 %>% map(~ jpn_pref(pref_code = ., district = FALSE)) %>% reduce(rbind) %>% ms_simplify(keep = 0.01, keep_shapes = TRUE) touhoku_map %>% ggplot() + geom_sf()
蛇足：{jpndistrict}がパッケージをロードしなければ使えないのは……
{jpndistrict}をつかって地図データを作るときに、jpn_pref()を一回しか使わないのでjpndistrict::jpn_pref()として実行したいがならぬものはならぬのです。
試しにやってみる。
# 明示的に`{jpndistrict}`をアンロードしておく。 detach("package:jpndistrict", unload = TRUE) # 都道府県JISコード1の北海道の地図データを作ろうとしてみる。 hokkaidou_map <- jpndistrict::jpn_pref(1, district = FALSE)
## Error in dplyr::filter(jpnprefs, jis_code == code_validate(code)$code): オブジェクト 'jpnprefs' がありません
エラーでストップする。
これは、多分この関数内で参照しているデータがR/sysdat.rdaに格納されているのではなくdata/jpnprefs.rdaに格納されているためである。
{jpndistrict}の作成者でもあるu_ribo氏のブログにもあるように、
パッケージを利用するユーザーが利用できるように提供するデータはdata/フォルダに置かれる。data()関数で呼び出されるデータはこれに当たる。dataフォルダには.rdataもしくは.rdaファイル（data()関数で呼び出すことを想定する）を置く。これらのファイルに保存されたRオブジェクトはパッケージの呼び出しとともに利用可能になる。
とのことで、data/以下に格納されているデータはパッケージがロードされるまで基本的いは参照できないため、::で関数を呼び出してもデータが参照できずエラーで止まってしまう様子。こんなの気にするのはパッケージ作ったりする人だけかもしれないが。
最初::で関数呼び出したコードがエラーでて、library()でロードしたときはうまくいくのがわからなくてちょっとハマった。
補足：地図データの入手
地図データの入手は今回はjpndistrict::jpn_pref()をつかって作った。
このデータを使う限り、“Map data is based on the Digital Map 25000 (Map Image) published by Geospatial Information Authority of Japan”と地図データの出典を記載できる。しかし、上記の関数を使う手法はやや重たい感が否めない。
もちろん、国土地理院のページから生データをダウンロードして処理しても良いと思うし、その他でも活用可能なデータを公開しているところは多い。
また、すでに公開されている整形済みデータを活用するのも一つの手だと思う。例えば、三重大学奥村晴彦氏のページで公開されているデータで
# japan2が地図データ japan2 = read_sf("https://oku.edu.mie-u.ac.jp/~okumura/stat/data/japan2.geojson")
でもよいと思う。
ただ、もともとのデータがデフォルメされすぎているデータだと全国地図だとそんなに違和感がないが、地方レベルにすると違和感が出てくるデータがあるので注意が必要だとおもう（適所の問題）。
地方レベルの地図を描画するのに向いていないデータの一例を次に示すが、データが悪いのではなく適所ではないということであしからず。
{NipponMap}の都道府県の境界データを使った場合の{ggplot2}での全国版の白地図は次のようになる（{NipponMap}をインストールしていなかったらinstall.packages("NipponMap")でインストール）4。
map <- read_sf(system.file("shapes/jpn.shp", package = "NipponMap")[1], crs = "+proj=longlat +datum=WGS84") map %>% ggplot() + geom_sf()
全国版だとそんなに違和感はないしデータが軽いので描画が早い5。このデータは、地方データがあるので地方でデータを指定すると、東北地方などの地図が簡単に書ける。
map %>% filter(region == "Tohoku") %>% ggplot() + geom_sf()
アップになるとちょっとデフォルメされすぎている感が強い。特に岩手から宮城のあたりのリアス式海岸のギザギザがギザギザハートすぎる印象があるし、何よりも図のサイズのせいもあるが宮城県の牡鹿半島が島になっているように見える（ギリギリつながっているが）（網地島はこんなに大きくない）。
地図の範囲に合わせた適所なデータを使うように心がけたい。
参考
シェープファイルを読む は三重大学奥村晴彦氏のページで、座標系をUTMに変換したり、奥村氏が変換済みのデータのことについてもこのページに記載がある。
GitHub - uribo/jpndistrict: Create Japansese Administration Area{jpndistrict}についてはこちらにある。u_ribo氏は{kuniezu}など地理データに係るRのパッケージを多数作成しており、ブログ上の記事も多い。大変参考になる。なお現在、{jpndistrict}、{kuniezu}ともにCRANではARCHIVE扱いとなっているが、solarisと古いmac版でビルドエラーでARCHIVE扱いになっているので通常に使うにはARCHIVE版でもGithub版でも特に問題なく使えている。
R でコロプレス図 (色分け地図) をなるべく簡単に描く -ill-identified diary簡易の地図を描くのに大変参考にした。その他のリンクも参考になった。
-
{sf}のst_simplify()が使われる記事を見かけることが多いが、km単位にしないと警告文がでるし、{rmapshaper}のms_simplify()のほうがよいアルゴリズムらしい。↩ -
このコードは3.2.Data Visualization [2]を参考（ほぼ借用）にした。
{tidyverse}を使ってmap()やreduce()で書いているが、lapply()やdo.call()をつかってもいいと思う。多分後者関数のほうが早いと思う。追記：lapply()やdo.call()を使うと次のようになる、japan_map <- rmapshaper::ms_simplify(do.call(rbind, lapply(1:47, jpn_pref, district = FALSE)), keep = 0.01, keep_shapes = TRUE)こちらのほうが圧倒的に早いと思う。：追記終了↩ -
らしい。ヘルプを読んでもよくわからない点だがこれでそれっぽくなる。↩
-
{ggplot2}を使わないで白地図を描きたかったら、NipponMap::JapanPrefMap(inset = FALSE)だけで描ける。↩ -
これはR でコロプレス図 (色分け地図) をなるべく簡単に描く -ill-identified diaryに記載されていたコードの一部から借用↩
