2020-05-04

ある座標からの指定半径に含まれるメッシュコードを知る

jpmesh sf mapview

新型コロナウイルスのデータを扱う際に、メッシュコード（標準地域メッシュ）が利用されることがあります。特にNTTドコモ「モバイル空間統計」分析レポートのデータは、内閣官房の新型コロナウイルス感染症対策のページにも掲載されているように全国各地の人口変動を分析するのに欠かせないデータとなっています。

人流データを使った分析は、先日発表された「新型コロナウイルス感染症対策の状況分析・提言」（2020年5月1日） (PDF)の中でも行われており、本文中に以下の記述があります。

渋谷駅周辺と難波駅周辺から半径 1 ㎞圏内においては、10 歳台および 20 歳台の若者を中心として昼夜問わず接触頻度が 80％以上、減少したことがうかがえる。

また、GitHubにアップロードされているこの結果の補足資料を見るとモバイル空間統計のメッシュデータを使った分析と可視化の事例が確認できます（3. 各エリアの接触頻度と変化率）。

補足資料の図2など、対象の駅周辺でのメッシュごとの接触の変化率を表示する際、駅周辺の半径1km圏内とそれに含まれる500mメッシュの枠が示されています。

緊急事態宣言中の人流データを整理した資料をクラスター対策班がまとめたので共有します。

難しい内容になっていますので、読むのに時間がかかるかもしれませんが、質問いただければ私たちとしてできるだけ回答します。がんばります。

▼緊急事態宣言中の人流データーまとめhttps://t.co/NKiVE8oFYC
— 新型コロナクラスター対策専門家 (@ClusterJapan) 2020年5月1日

github.com

この図を見た直後、これは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

f:id:u_ribo:20200504135737p:plain

座標とバッファ領域が表示されました。

続いて、この範囲に含まれる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)

f:id:u_ribo:20200504140722p:plain

バッファ領域とそれに含まれるメッシュの抽出ができました。いくらか手間がかかりますが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))

f:id:u_ribo:20200504140332p:plain

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"))

f:id:u_ribo:20200504140229p:plain

以上です。

2020-04-25

R Markdownでのハイライトを行うflairパッケージを使ってみた

package flair

f:id:u_ribo:20200425120344p:plain

R MarkdownによるHTML出力を行う時、コードのシンタックスハイライトを取り入れることができます。一方でコードの特定の箇所を強調したり装飾することは困難でした(xaringanパッケージでのプレゼンテーション用の出力では可能)。

ですがそれも過去の話。今日紹介する flair パッケージを使うと、そうした要望に答えられます。

github.com

公式のドキュメントが整備されているのでそっちを読んだ方が手っ取り早い、正確な気もしますが備忘録として日本語で整理しておきます。

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!

2020-04-22

tabularmaps: カラム地図で行政区の大まかな配置を可視化する

tabularmaps development mapping

カラム地図と呼ばれるものがあります。これは日本の47都道府県をはじめとした行政区の配置を表上に圧縮表示することで、それぞれの位置関係をわかりやすく伝えるためのプロジェクトです。

f:id:u_ribo:20200422175004p:plain

最近では、カラム地図の開発者の一人、福野泰介さん (@taisukef) による新型コロナウイルス対策ダッシュボードのページでの「現在患者数 / 対策病床数」の可視化で「カラム地図」を見る機会が増えた人もいるのではないでしょうか。

COVID-19 Japan 新型コロナ対策ダッシュボード
厚労省データと、神奈川県が公表した対策病床3,400反映https://t.co/VNM1ngPG95
病床等確保数1,000床＋宿泊療養施設2,400室https://t.co/kVGJbYMNTx \#StopCOVID19JP
対策病床数のオープンデータ化、お願いします！
— taisukef (@taisukef) April 21, 2020

さて、この「カラム地図」をR言語の中で扱うべく、tabularmapsパッケージの開発を進めています。

github.com

以下、この記事ではtabularmapsの使い方を紹介していきます。

パッケージはCRANには登録されていないため、GitHub経由でインストールを行います。以下のコマンドを実行することでパッケージが利用可能になります。

install.packages("remotes")
remotes::install_github("uribo/tabularmaps")

library(tabularmaps)
library(ggplot2)

可視化の際にggplot2パッケージを利用しているので、そちらも読み込んでおく必要があります。

使い方

現在tabularmapsでは

都道府県 (jpn77)
東京23区 (tky23)
ISO-3166による国名 (iso3166)

の3種類に対応しています。括弧内の文字列がパッケージに含まれるデータフレームオブジェクトで、県や区などの対象行政名の位置関係を記録しています。

冒頭の都道府県のカラム地図は次のコマンドにより描画されます。

tabularmap(jpn77, 
           fill = region_kanji, 
           label = prefecture_kanji, 
           size = 3,
           family = "IPAexGothic") +
  theme_tabularmap(base_family = "IPAexGothic") +
  scale_fill_jpregion(lang = "jp",
                      name = "八地方区分")

tabularmap()がカラム地図を表示させる関数として機能します。この関数では第一引数に対象のデータフレームを与えて実行します。そのほかの引数としてggplot2オブジェクトの審美的要素を決定するfill (塗り分けの色基準)、label (表示するラベル文字列)の変数名を与えます。これに加え、ggplot2::geom_text()に渡す値を指定可能です。これはラベルの大きさを調整する size やフォントのための family を使うことを想定しています。

上記のコードでは、tabularmap()に加えてtheme_tabularmap()、scale_fill_jpregion()を実行しています。これらの関数はオプションとして利用可能なもので、よりカラム地図っぽい見た目に調整するためのものです。都道府県を表示させるカラム地図では地方別に塗り分ける色が決まっているため、それに応じた塗り分けを行うようになります。

続いて東京23区の表示をしてみます。先のコードからtabularmap()に与えるデータとラベル用の変数名を変更するだけです。

tabularmap(tky23,
           label = ward_kanji,
           family = "IPAexGothic") +
  theme_tabularmap(base_family = "IPAexGothic") +
  guides(fill = FALSE)

現実の配置や面積がぼやけてしまう欠点もありますが、行政区域のデータを用意する必要がないのはカラム地図の大きな利点です。（カラム地図のレイアウト自体はオープンライセンスで開発されています。）

また本来想定されていたダッシュボードでの利用を考える場合、表示する際のスペースが有効活用できるのも嬉しいです。日本の場合、細長くなるのでどうしても隙間ができてしまいますが、カラム地図だと基本的にグリッドになるので無駄が少ないです。

今後の機能

fukuno.jig.jp

本家のカラム地図では、市区町村を対象にしたレイアウトが全国47都道府県分定まったとのことです。これらの情報を追いながら、tabularmapsパッケージでも市区町村版を追加する作業を進めていくつもりです。

全国のデータを追加するのは時間がかかりそうなので、皆さんからの協力(Pull request)をお待ちしています。

github.com

余談

SNSを見ていると、時々ヘンテコな地図が回ってきます。

「白地図上に47都道府県の名前を書き込んでください」

中学生の地理で習うような問題です。回答を見ていると、有名な東京都や神奈川県、形のわかりやすい北海道などは簡単な様子。一方で北関東の群馬、栃木、茨城の配置、島根と鳥取の位置関係などを間違えている例をしばしば見かけます（私も苦手でした）。またそれ以上にトンデモな解答も…。

行政境界の複雑な幾何学模様を認識し記憶することに対する人間の限界を感じます。それが特に自分とは関わりのない地域であれば尚更です。

また、行政区は人間が勝手に定めたもの。時間の経過とともに境界も変化していきます。

情報を伝えるとき、これらの行政区をどこまで厳密に再現する必要があるでしょうか。

地図を表示するにはポリゴンないしラインデータを用意しなくてはいけない。そう思っていた時期が私にもありました。

ですがそうではなかったのです。

必要なのは、情報を正しく、わかりやすく伝えること。

カラム地図では、現実に対して幾らかの嘘をつくことになりますが、それは許容されうる範囲なのではと思います。

cucumber flesh

Rを中心としたデータ分析・統計解析らへんの話題をしていくだけ

ある座標からの指定半径に含まれるメッシュコードを知る

Rコード

おまけ: 出力のための調整

R Markdownでのハイライトを行うflairパッケージを使ってみた

導入

flairの派生系

文字装飾

tabularmaps: カラム地図で行政区の大まかな配置を可視化する

使い方

今後の機能

余談

参考