所要時間のヒートマップを作成する

Transcript

1. 所要時間のヒートマップを 作成する google directions API & googleway package Tomokazu FUJINO @ R研究集会2018 LT December 28, 2018
2. はじめに → 不動産の賃料推定などで，主要都市などへの所要時間を考慮に 入れたいことがある → 多くの地点から，特定の地点までの所要時間データの作成を自 動化したい → 今回の例：福岡県の 3 次メッシュにおける代表点から，天神へ の所要時間を可視化する 2
3. 道具 → Google Maps Platform の Directions API → R から呼び出すパッケージ：googleway package → 土地利用 3 次メッシュ＠国土数値情報 → kokudosuuchi , sf package → rmapshaper package 3
4. Directions API → Google Map Platform 経路探索 API → 月 40000 回呼び出しま で無料で使える → 日本では電車の経路は使 えない（徒歩，バス， 車，自転車のみ） 4
5. ライブラリ読み込み library(tidyverse) library(sf) library(googleway) 5
6. 利用データの識別子の確認 → 国土数値情報の「土地利用 3 次メッシュ」と「行政区域」の データを API 経由でダウンロードするため，識別子を調べる。 kokudosuuchi::getKSJSummary() %>% filter(title == "土地利用 3 次メッシュ") ## # A tibble: 1 x 5 ## identifier title field1 field2 ## <chr> <chr> <chr> <chr> ## 1 L03-a 土地利用~ 国土（水・~ 土地利用~ ## # ... with 1 more variable: ## # areaType <chr> kokudosuuchi::getKSJSummary() %>% filter(title=="行政区域") ## # A tibble: 1 x 5 ## identifier title field1 field2 ## <chr> <chr> <chr> <chr> 6
7. 福岡県を覆う3次メッシュを取ってくる → st_transform() 関数で平面直角座標系に変換 kokudosuuchi::getKSJURL("L03-a", meshCode = c(5030, 5031, 4930)) %>% filter(year == 2014, datum == 2) %>% pull(zipFileUrl) %>% map(kokudosuuchi::getKSJData) %>% flatten %>% data.table::rbindlist() %>% st_as_sf() %>% st_transform(crs = 2444) -> kyushu.mesh 7
8. 福岡県の行政界を取ってくる kokudosuuchi::getKSJURL("N03", prefCode = 40) %>% filter(year == 2014) %>% pull(zipFileUrl) %>% map(kokudosuuchi::getKSJData) %>% flatten %>% data.table::rbindlist() %>% st_as_sf() %>% st_transform(crs = 2444) -> fukuoka.boundary 8
9. メッシュの切り出しと中心点の計算 福岡県の行政界でメッシュを切り出す st_join(kyushu.mesh, fukuoka.boundary, join = st_intersects) %>% filter(!is.na(N03_007)) -> fukuoka.mesh メッシュの中心点を計算して，出発点として使う fukuoka.mesh %>% mutate(center = st_centroid(geometry)) %>% st_set_geometry("center") %>% st_transform(crs = 4326) %>% mutate(x = st_coordinates(center)[,2], y = st_coordinates(center)[,1]) -> fukuoka.centroid 9
10. 明らかに不要なメッシュを除去 & APIキーの 設定 fukuoka.centroid %>% filter(海水域 < 500000) %>% select(x, y) %>% as.data.frame() %>% select(-center) -> fukuoka.origin set_key("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX") 10
11. API呼び出し → google_directions() 関数を使う → mode = "transit"で徒歩＋バスの経路検索（日本の場合） fukuoka.origin %>% pmap(function(x, y){ tryCatch({ google_directions(origin = data.frame(x, y), destination = c(33.591679, 130.398177), mode = "transit", departure_time = as.POSIXct("2018-12-10 07:00:00")) }, error = function(e){ return(NA) }) }) -> tenjin.duration 11
12. 結果のリストを処理する関数を定義 get_duration <- function(dlist){ if("status" %in% names(dlist)){ if(dlist$status == "OK"){ direction_legs(dlist) %>% select(duration, start_location, end_location) %>% rlist::list.flatten() } } } Unzip <- function(...) rbind(data.frame(), ...) 12
13. データフレームにしてsfにする メッシュと結合するため平面直角座標系にしておく tenjin.duration %>% map(get_duration) %>% do.call(Unzip, .) %>% st_as_sf(coords = c("start_location.lng", "start_location.lat"), crs = 4326) %>% st_transform(crs = 2444) -> tenjin.duration.df メッシュと結合 st_join(fukuoka.mesh, tenjin.duration.df, join = st_intersects) %>% filter(!is.na(duration.value)) -> tenjin.duration.mesh 13
14. 地図の描画 # 小さい離島を除いておく library(rmapshaper) fukuoka.boundary <- ms_filter_islands( fukuoka.boundary, min_area = 2000000 ) ggplot() + geom_sf(data = fukuoka.boundary, aes(fill = NULL)) + geom_sf(data = tenjin.duration.mesh, aes(fill = cut(duration.value, quantile(duration.value)))) + coord_sf(datum = NA) + theme_void() 14
15. 地図の描画 fill (445,4.66e+03] (4.66e+03,6.66e+03] (6.66e+03,8.36e+03] (8.36e+03,3.73e+04] NA 15
16. 車の場合: mode="driving" fill (281,3e+03] (3e+03,3.67e+03] (3.67e+03,4.69e+03] (4.69e+03,1.42e+04] NA 16