ある座標からの指定半径に含まれるメッシュコードを知る
新型コロナウイルスのデータを扱う際に、メッシュコード（標準地域メッシュ）が利用されることがあります。 特に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"))
以上です。
