kuniezu: 日本の国土地理を扱いやすくするRパッケージをCRANに登録しました
はじめに
kuniezuパッケージ (v0.1.0) をCRANにリリースしました。
このパッケージは、私が業務や趣味で日本国内の地理空間データを扱う時に作っていた関数を一つのパッケージに整理したものです。 空間的に世界規模のデータを扱うのではなく、日本国内に限った話であれば、日本に即した仕様や座標参照系を利用した方が良いことがあります。 そうした日本の地理空間データを処理する際に利用することがある機能や、あると便利なデータセットを提供できるように努めています。
ゆるゆると開発するつもりでいましたが、Twitterで開発中であることをつぶやいたところ想定以上に反響が大きかったのでCRAN登録を急ぐことにしました。
仕事でも使うし、誰かの役に立つかと思い、日本の地理空間データのための処理を関数化し、パッケージとして開発始めました。🗾ブログに書いた、度分秒での経緯度の変換や平面直角座標系の特定とか。機能は少ないですが使ってみてください。（フィードバック、機能要望もぜひ） https://t.co/VUaPyC6MaT pic.twitter.com/bwaFL0vJqw
— Uryu Shinya (@u_ribo) 2020年4月27日
この記事では、そんなkuniezuパッケージの機能紹介を行います。
まずはみなさんパッケージがインストールされていないと思うので、利用されたい方は以下のコマンドを実行してインストールを済ませてください。
install.packages("kuniezu")
library(kuniezu) # 使い方を説明するために以下のパッケージも読み込んでおきます library(sf) library(ggplot2) library(leaflet)
使い方
最初のリリースでは以下の機能を実装しました。
- DMS表記を十進数表記に変換
- 日本測地系2011における平面直角座標系の特定
- 南西諸島・小笠原諸島を移動した日本地図の描画
- 地理院タイルをleafletで簡単に利用できるように
- 国土地理に関するいくつかのデータセット
順に解説します。
parse_*_dohunbyo(): DMS表記を十進数表記に変換
地球上の任意の座標を示す際に使われる緯度と経度の表記には、十進数で表される場合と「度(Degree)・分(Minute)・秒
(Second)」を用いて表記される場合があります。
後者はDMS表記と呼ばれ、十進数で139.7414、35.6581と表される座標に対してE139°44′28.8869″、N35°39′29.1572″のように表現します。
更に日本では度分秒の記号に対して漢字使って東経139度44分28.8869秒、北緯35度39分29.1572秒等の表記が用いられることがしばしばあります。
この日本独自のDMS表記を十進数に変換する関数がparse_*_dohunbyo()です（アスタリスクの部分にはlon、latがそれぞれ指定できます）。
lonは東経、西経を示す(180~-180）経度、latが北緯、南緯（90~-90）の緯度を扱うために使い分けます。
# 日本経緯度原点（東京都港区麻布台二丁目）の座標 parse_lon_dohunbyo("東経139度44分28.8869秒")
## [1] 139.7414
parse_lat_dohunbyo("北緯35度39分29.1572秒")
## [1] 35.6581
これは過去のブログ記事に書いた関数をパッケージに含めたものです。
加えて、十進数に変換する前段階として漢字を記号に処理するreplace_dohunbyo_kanji()も用意しています。
replace_dohunbyo_kanji("東経139度44分28.8869秒")
## [1] "E139°44’28.8869."
日本測地系2011における平面直角座標系の特定
日本の現行の測地系である日本測地系2011には、19の平面直角座標系の区分がなされています。
平面直角座標系は狭い範囲の測量や距離の算出に用いるのに適しています。
任意の経緯度を平面直角座標に変換するには、まず座標が平面直角座標系のどの系番号に含まれるのかを知っておく必要が生じます。
これを求めるのがst_detect_jgd2011()です。
sfパッケージでおなじみのst_*()が関数名につくように、関数にはsfcオブジェクトを与えて実行します。
# 日本経緯度原点の座標 st_detect_jgd2011(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326))
## [1] 6677
move_jpn_rs(): 南西諸島の一部・小笠原諸島を移動した日本地図の描画
日本だけを地図上に描画する際、表示する画面の都合からか、南西諸島の一部や小笠原諸島の位置を海上に移動することがあります（例えばNHKの気象情報のページや白地図）。 これを実現するには、予め移動対象の地物に対する操作が必要となります。
move_jpn_rs()はこれらの島々の位置を変更するための関数です。
地図描画のために奄美群島と沖縄県、小笠原諸島の位置を変更する処理が実行されます。
また、この地図化もパッケージの関数で行えるようにggplot2ベースのgeom_jpsegment()を用意しました。
使い方は次のように、まず都道府県や市区町村レベルでポリゴンが分かれた地理空間データに対してmove_jpn_rs()を適用、その後マッピングを行います。
move_jpn_rs(jgd2011_bbox) %>% ggplot() + geom_sf() + geom_jpsegment()
デフォルトの設定で移動した地物がわかるように境界線が引かれます。
この関数の元も、かつて記事に書いたものです。
少し前に @shoei05 さんがブラッシュアップしてくれたのをきっかけにパッケージへ含めました。
@u_ribo いつも勉強させていただいています。以前の記事（https://t.co/iJxfDxFhJ3 ）に関連して、コロプレスマップなど作図する場合は、sf_ja_omit47 とするのではなく、やや強引ですが同一セット内の $geometry に座標を足してあげると良さそうだったので、少しスクリプトをいじってみました。 pic.twitter.com/EI9NvHMpnD
— Shoei/将英 (@shoei05) March 26, 2020
一方で、まだこの関数はさらなる磨き上げが必要で、解決策がわかる人はぜひ教えてください。
地理院タイルをleafletで簡単に利用できるように
インタラクティブなウェブ地図を利用したことがある方ならわかるかと思いますが、leafletでは任意の地図タイルを背景画像として扱うことができます。 国土地理院が発行する地理院タイルもまた、この背景画像として利用可能で、 一覧ページに示されるように多種多様な地図タイルが整備されています。
kuniezuパッケージではこの地理院タイルをleafletで簡単に導入できるようになっています。 これらはgsi_tilesオブジェクトを通して呼び出します。
現在のバージョンでは「近年の災害」をはじめとして、すべての地理院タイルに対応しているわけではありません。 kuniezuで利用可能な地理院タイルの一覧は以下のコードで確認できます。
names(gsi_tiles)
## [1] "standard" "pale" "english"
## [4] "lcm25k_2012" "lcm25k" "ccm1"
## [7] "ccm2" "vbm" "vbmd_bm"
## [10] "vbmd_colorrel" "vbmd_pm" "vlcd"
## [13] "lum200k" "lake1" "lakedata"
## [16] "blank" "seamlessphoto" "ortho"
## [19] "airphoto" "gazo4" "gazo3"
## [22] "gazo2" "gazo1" "ort_old10"
## [25] "ort_USA10" "ort_riku10" "relief"
## [28] "anaglyphmap_color" "anaglyphmap_gray" "hillshademap"
## [31] "earthhillshade" "slopemap" "slopezone1map"
## [34] "afm" "lcmfc2" "lcmfc1"
## [37] "swale" "cp" "jikizu2015_chijiki_d"
## [40] "jikizu2015_chijiki_i" "jikizu2015_chijiki_f" "jikizu2015_chijiki_h"
## [43] "jikizu2015_chijiki_z" "jikizu_chijikid" "jikizu_chijikii"
## [46] "jikizu_chijikif" "jikizu_chijikih" "jikizu_chijikiz"
使い方は次のようになります。leafletおよびmapviewでの例を紹介します。
# 出力結果は省略します gsi_tiles$standard %>% addCircles( data = sf::st_transform(extreme_points %>% purrr::reduce(c), crs = 4326)) mapview::mapview(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326), map = gsi_tiles$pale)
なお、これらの地理院タイル利用する際はタイルごとに設定されている利用規約に従ってください。
国土地理に関するデータセット
これまでの例で出てきたものがありますが、日本の国土地理を扱う上であると便利なデータ等をパッケージで提供するようにしました。 現在、次のようなデータがあります。
- extreme_points… 東西南北端点の経度緯度データ
- jgd2011_bbox… 日本測地系2011の平面直角座標系の区域を示すポリゴンデータ
- jp47prefectural_offices… 47都道府県の庁舎の位置データ
こんなデータや機能が欲しいという方はぜひ issue に書き込んでください。 バグ報告もありがたいです。 GitHub以外ではTwitterでコメントいただいても対応します。
Enjoy!
ある座標からの指定半径に含まれるメッシュコードを知る
新型コロナウイルスのデータを扱う際に、メッシュコード（標準地域メッシュ）が利用されることがあります。 特にNTTドコモ「モバイル空間統計」分析レポートのデータは、内閣官房の新型コロナウイルス感染症対策のページにも掲載されているように 全国各地の人口変動を分析するのに欠かせないデータとなっています。
人流データを使った分析は、先日発表された「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日） (PDF)の中でも 行われており、本文中に以下の記述があります。
渋谷駅周辺と難波駅周辺から半径 1 ㎞圏内においては、10 歳台および 20 歳台の若者を中心として昼夜問わず接触頻度が 80％以上、減少したことがうかがえる。
また、GitHubにアップロードされている
この結果の補足資料を見るとモバイル空間統計のメッシュデータを使った分析と可視化の事例が確認できます（3. 各エリアの接触頻度と変化率）。
補足資料の図2など、対象の駅周辺でのメッシュごとの接触の変化率を表示する際、駅周辺の半径1km圏内とそれに含まれる500mメッシュの枠が示されています。
緊急事態宣言中の人流データを整理した資料をクラスター対策班がまとめたので共有します。
— 新型コロナクラスター対策専門家 (@ClusterJapan) 2020年5月1日
難しい内容になっていますので、読むのに時間がかかるかもしれませんが、質問いただければ私たちとしてできるだけ回答します。がんばります。
▼緊急事態宣言中の人流データーまとめhttps://t.co/NKiVE8oFYC
この図を見た直後、これはRでできるぞー！と思いました（もともとRでやられている、より効率的な処理を実行されているのかもしれませんが...）。 ブログの主題にあるように「ある座標からの指定半径に含まれるメッシュコードを知る」ことができれば簡単に実行できます。
真似事になりますが、メッシュコードを使った分析例として優れており、応用範囲の広いトピックスですのでRでやる方法を整理しておきます。
Rコード
library(dplyr) library(mapview) library(sf) library(jpmesh)
まず対象の座標からPOINTのsfcオブジェクトを作成します。
sfパッケージではst_point()に任意の座標をベクトルで与えることでPOINTの作成が行われます。
これを測地系WGS84の座標参照系、EPSG:4326として扱えるようにst_sfc()で指定します。
st_point()で定義する座標は、500mメッシュコードに変換した時に 533935961 となる渋谷駅の経緯度です。
p_shibuya_st <- # 目測で座標を決めているので厚生労働省のポイントとはズレがあります st_point(c(139.70172, 35.65846)) %>% st_sfc(crs = 4326)
続いてこの座標点を起点としたバッファ領域（緩衝帯）を生成します。これにはst_buffer()を使いますが、EPSG:4326の地球楕円体モデルは先述の通りWGS84です。
これは水平位置を表す経緯度と垂直位置を表す高度との組み合わせからなる、3次元の地理座標系です。
2地点間の距離やバッファを求める際には地物を平面に投影し、XY座標で表す投影座標系である平面直角座標系を利用すると良いので st_transform()により座標参照系の変更を行います。
ここでは日本測地系2011の平面座標系の一つであるEPSG:6677（福島県、栃木県、茨城県、埼玉県、千葉県、群馬県、神奈川県、島嶼部を除いた東京都が含まれる）を指定しました。
b1km_shibuya_st <- p_shibuya_st %>% st_transform(crs = 6677) %>% st_buffer(dist = units::set_units(1, km)) %>% st_transform(crs = 4326)
上記の例ではバッファの範囲を半径1kmとしました。処理のあと、st_transform()により再び元の座標参照系に戻しています。
一度確認してみましょう。白地図では位置関係が掴みにくいので地図タイルを背景に、データを重ます。
m <- mapview(p_shibuya_st) + mapview(b1km_shibuya_st) m
座標とバッファ領域が表示されました。
続いて、この範囲に含まれる500mメッシュを特定します。 ここではメッシュコードが未知のものとして座標を与えて探し出すところから行う例を示します。 座標の位置する80kmメッシュから、それに含まれる10kmメッシュコードの生成、ポリゴン化を以下のコードで行います。 なお1kmメッシュから探索を開始しないのは、候補となるメッシュの数を減らし、実行時間を短縮するためです。
mesh_candidate <- p_shibuya_st %>% coords_to_mesh(geometry = ., mesh_size = 80) %>% mesh_convert(to_mesh_size = 10) %>% export_meshes() nrow(mesh_candidate)
10kmメッシュのポリゴンが用意できたら、バッファのポリゴンに対する空間関係を調べます。
st_join()は2つの地物の空間関係をもとにデータを結合するために使われます。
今回は2つの地物（10kmメッシュ、バッファポリゴン）に共有部分がある場合にデータを結合するように join = st_intersectsを引数で与えます（st_join()の既定値）。
さらに結合後のデータを共有部分だけのものにしたいのでleft = FALSEを指定しました。
mesh_candidate <- mesh_candidate %>% st_join(st_sf(b1km_shibuya_st), join = st_intersects, left = FALSE) mesh_candidate
実行結果を見ると、10kmメッシュのデータが2件に絞られています。 これはバッファポリゴンと交差するポリゴンを持っている10kmメッシュコードが2つあることを示します。
残りの処理は対象のメッシュに対してスケールダウン、バッファとの空間関係を調べて抽出するを繰り返すだけです。 ただし、10kmから1km、500mのスケールダウンを何度も行うのは手間ですので、ここで10kmメッシュに含まれる500mメッシュを用意して処理を簡略化します。 1つの10kmメッシュコードに含まれる500mメッシュの数は400、そのため800個の500mメッシュから最終的なメッシュコードを抽出することになります。
mesh_target <- mesh_candidate %>% pull(meshcode) %>% purrr::map( ~ mesh_convert(.x, to_mesh_size = 0.5) ) %>% purrr::reduce(c) %>% unique() %>% export_meshes() %>% st_join(st_sf(b1km_shibuya_st), join = st_intersects, left = FALSE) nrow(mesh_target)
最終的な渋谷駅周辺の半径1kmに含まれる500mメッシュの数は22個となりました。 確認のためにもう一度データを表示してみます。
m + mapview(mesh_target)
バッファ領域とそれに含まれるメッシュの抽出ができました。 いくらか手間がかかりますがRでもできます、という話でした。
Enjoy!
おまけ: 出力のための調整
このマップをさらに報告書などのPDFに添付することがある場合を想定し、その際の調整方法を書いておきます。 具体的にはマップに対して
- 凡例の変更... 元のマップには凡例もないので載せる
- ファイルへの出力... 厚生労働省の補足資料では切り抜きが雑なので統一する
を行います。これらもRで実行します。
まずはマッピングするデータを用意します。手元にデータがないので図2 左の値をデータ化します。
df_value <- tibble::tibble( meshcode = as.character(c(533945052, 533945061, 533935953, 533935954, 533935963, 533935964, 533935951, 533935952, 533935961, 533935962, 533935971, 533935853, 533935854, 533935863, 533935864, 533935873, 533935851, 533935852, 533935861, 533935862, 533935754, 533935763)), value = c(-51, -60, -22, -49, -53, -54, -35, -60, -62, -57, -33, -35, -55, -60, -52, -31, -15, -35, -38, -38, -27, -35))
続いて凡例で使われるカラーパレットを定義しておきます。変化率は正負いずれの値も取り得るとし、 原点(0)を中心に赤系（増加）と青系（減少）を示す10刻みのカラーパレットを作成します。
class_int <- classInt::classIntervals(seq.int(-100, 100, by = 10), n = 20, style = "fixed", fixedBreaks = seq.int(-100, 100, by = 10)) # カラーパレットの確認 pals::pal.bands(pals::ocean.balance(20))
df_map <- mesh_target %>% left_join(df_value, by = "meshcode") %>% mutate(value_class = cut(value, class_int$brks, include.lowest = TRUE)) res_map <- df_map %>% mapview(zcol = "value_class", layer.name = "変化率", col.regions = pals::ocean.balance(20), homebutton = FALSE) %>% leafem::addStaticLabels( data = df_map, textsize = "18px", label = paste0(df_map$value, "%"))
mapshot(res_map, file = "out.png", remove_controls = c("zoomControl", "layersControl", "homeButton"))
以上です。
R Markdownでのハイライトを行うflairパッケージを使ってみた
R MarkdownによるHTML出力を行う時、コードのシンタックスハイライトを取り入れることができます。一方でコードの特定の箇所を強調したり装飾することは困難でした(xaringanパッケージでのプレゼンテーション用の出力では可能)。
ですがそれも過去の話。今日紹介する flair パッケージを使うと、そうした要望に答えられます。
公式のドキュメントが整備されているのでそっちを読んだ方が手っ取り早い、正確な気もしますが備忘録として日本語で整理しておきます。
flairパッケージでできること
- R MarkdownによるRコードのHTML出力の部分ハイライト
- ハイライト対照のコードはコードチャンク、文字列で指定
decorate()でチャンクオプションの制御flair()- 特定の文字列、関数、引数をハイライトできる派生関数も豊富
- ハイライトの文字の色、背景色、斜体や太字などの文字装飾を調整可能
なお、ハイライトは便利でわかりやすさを促進する可能性がありますが、ハイライトによって見えにくくなったり混乱を招く恐れがあることも注意が必要です。これについて、公式のドキュメント内で
However, please remember to be judicious in your color choices, and to keep in mind how your colors appear to colorblind individuals.
の一文があるのがとても好印象でした。用法要領を守って適切に使っていきたいですね。
導入
flairは先日CRANに登録されました。Rの標準パッケージ追加関数であるinstall.packages()を使ってインストールが可能です。開発版を入れたい人は
remotes::install_github("kbodwin/flair)を実行してください。
library(flair) library(dplyr)
flairではdecorate()とflair()が主要な関数となります。あとで触れますがflair()にはいくつかの派生系が用意されています。基本的な使い方を見ていきましょう。
まずはいつも通りチャンクを使ってRコードを書きます。このコードがハイライトの対象となります。この時、チャンクラベル（チャンク名）をつけておくのが大事です。
```{r how_to_pipe, include = FALSE}
# how_to_pipeがチャンク名になります
iris %>%
group_by(Species) %>%
summarize(sepal_length_avg = mean(Sepal.Length))
```
次の処理がハイライトのための操作となります。まずdecorate()で先ほどの、ハイライトさせたいコードが書かれてたチャンクラベルを文字列で指定、そのオブジェクトをflair()に渡します。flair()の第一引数にdecorate()の結果が入ることになるのでパイプ演算子 (%>%)を使うと楽です。複数のハイライトを行う場合も、パイプ演算子でつなげて書けるので便利です。
そしてflair()関数ではハイライトさせる文字列を指定します。
decorate("how_to_pipe", eval = FALSE) %>% flair("%>%")
iris %>%
group_by(Species) %>%
summarize(sepal_length_avg = mean(Sepal.Length))上記のハイライトされている部分が出力結果です。KnitするとpreとcodeタグからなるHTMLが出力される仕組みです。実際の使い所としてはdecorate()、flair()のコードは見せなくても良いので、そちらのチャンクコードをecho=FALSEにしておくと良いと思います。
decorate()では、引数内でチャンクオプションの値を制御できます。上記のコードでは、irisデータの集計結果を出力しないためにeval = FALSEを与えました。
また、decorate()ではチャンクラベルではなく、Rコードを直接記述してハイライトの対象としても良いです。今度はチャンクコードのRコードも実行する様にeval = FALSEを与えずに実行してみます。
decorate("mean(1:10)") %>% flair("mean")
mean(1:10)## [1] 5.5
R MarkdownファイルをKnitする前に、どのような出力になるかを確認することも可能です。コンソールで上記のコードを実行してみましょう(Knitする前のチャンク名をどうして把握できているのか、よくわかっていない)。RStudioを使っている場合、Viewerパネルにハイライトされた結果が表示されます。
flairの派生系
先ほどはflair()で特定の文字列をハイライトする例でしたが、関数や引数をハイライトするための関数も用意されています。それぞれflair_funs()とflair_args()です。このほかにも指定した行をハイライトするflair_lines()、関数に与えた引数と値のためのflair_input_vals()、正規表現によるパターン指定が可能なflair_rx()などがあります。
decorate("how_to_pipe", eval = FALSE) %>% flair_funs()
iris %>%
group_by(Species) %>%
summarize(sepal_length_avg = mean(Sepal.Length))decorate("how_to_pipe", eval = FALSE) %>% flair_args()
iris %>%
group_by(Species) %>%
summarize(sepal_length_avg = mean(Sepal.Length))文字装飾
flairパッケージではハイライトしたときの見た目を変更する機能が備わっています。これはflair()内で指定しますが、論理値で指定するものと値を指定するものがあります。前者は太字
(bold)、下線 (underline)、斜体 (italic)とするか、後者はハイライトの色
(color) および背景色 (bg_color)などです。
decorate("how_to_pipe", eval = FALSE) %>% flair_lines(2) %>% flair_funs(color = "LightCoral") %>% flair("%>%", color = "Azure", background = "Navy")
iris %>%
group_by(Species) %>%
summarize(sepal_length_avg = mean(Sepal.Length))こんなてんこ盛りな装飾もできます。
Enjoy!
