複数の点から凸包？を生成するには{concaveman}が便利
Isochrone mapやサービスエリアマップなどの作成で点の集合の外側を結んでポリゴンを作成したい。（凸包というのかな？）
sf::st_convex_full()でも生成できるがmapbox社製concavemanのRインタフェースがパッケージ化されていたので試してみた。
- GitHub - mapbox/concaveman: A very fast 2D concave hull algorithm in JavaScript
- GitHub - joelgombin/concaveman: A very fast 2D concave hull algorithm
サンプルデータの取得
サンプルデータとして東京公共交通オープンデータチャレンジで公開されているドコモ・バイクシェアのポートデータを使用する。
library(dplyr) consumerKey <- "******" # シェアサイクルポートリスト parklist <- glue::glue("https://api-tokyochallenge.odpt.org/api/v4/docomo-cycle/parklist?acl:consumerKey={consumerKey}") %>% jsonlite::fromJSON() %>% # 不正データを除く dplyr::filter(lat > 0.0) %>% sf::st_as_sf(coords = c("lon", "lat")) %>% sf::st_set_crs(4326) parklist %>% head() %>% knitr::kable()
| basho | limit | ent_id | area_id | park_id | end_time | cycle_num | image_url | start_time | updated_at | setting_num | cycle_num_lv | park_name_en | park_name_jp | setting_num_lv | geometry |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 秋葉原・神田/Akihabara・Kanda | 0 | TYO | 1 | 00010909 | 00:00 | 4 | https://tcc.docomo-cycle.jp/cycle/park_img/00010909.jpg | 00:00 | 1595570524 | 33 | 1 | A3-42.Uchikanda 3Chome | A3-42.内神田3丁目 | 4 | POINT (139.7701 35.6926) |
| 高田馬場・早稲田/Takadanobaba・Waseda | 0 | TYO | 5 | 00010854 | 00:00 | 1 | https://tcc.docomo-cycle.jp/cycle/park_img/00010854.jpg | 00:00 | 1595570524 | 4 | 1 | D3-12.Bellus Vita | D3-12.Bellus Vita | 4 | POINT (139.7089 35.71087) |
| 大島/Ojima | 0 | TYO | 4 | 00010915 | 00:00 | 32 | https://tcc.docomo-cycle.jp/cycle/park_img/00010915.jpg | 00:00 | 1595570524 | 0 | 4 | H1-133.kitasuna2chome park | H1-133.北砂2丁目公園 | 0 | POINT (139.8258 35.679) |
| 飯田橋・九段下/Iidabashi・Kudanshita | 0 | TYO | 1 | 00010001 | 00:00 | 4 | https://tcc.docomo-cycle.jp/cycle/park_img/00010001.jpg | 00:00 | 1595570524 | 6 | 2 | A1-01.Chiyoda City Office | A1-01.千代田区役所 | 3 | POINT (139.7533 35.6938) |
| 飯田橋・九段下/Iidabashi・Kudanshita | 0 | TYO | 1 | 00010003 | 00:00 | 7 | https://tcc.docomo-cycle.jp/cycle/park_img/00010003.jpg | 00:00 | 1595570524 | 0 | 4 | A1-03.Tokyo Kusei Kaikan building(Metro A5 Exit) | A1-03.東京区政会館（メトロA5出口） | 0 | POINT (139.7475 35.70107) |
| 飯田橋・九段下/Iidabashi・Kudanshita | 0 | TYO | 1 | 00010004 | 00:00 | 10 | https://tcc.docomo-cycle.jp/cycle/park_img/00010004.jpg | 00:00 | 1595570524 | 0 | 4 | A1-04.Tokyo Kusei Kaikan building (Western parking lot) | A1-04.東京区政会館（西側駐輪場） | 0 | POINT (139.747 35.70135) |
地図上に可視化してみると次のようになる。
pal <- leaflet::colorFactor(c(RColorBrewer::brewer.pal(8, "Set1"), RColorBrewer::brewer.pal(5, "Pastel1")), parklist$area_id) # 可視化 map <- leaflet::leaflet(options = leaflet::leafletOptions(zoomControl = FALSE)) %>% leaflet::addTiles("http://cyberjapandata.gsi.go.jp/xyz/blank/{z}/{x}/{y}.png", attribution = "<a href='http://maps.gsi.go.jp/development/ichiran.html' target='_blank'>地理院タイル</a>") map %>% leaflet::addCircleMarkers(data = parklist, radius = 3, color = ~pal(area_id), weight = 1, fillColor = ~pal(area_id), fillOpacity = 0.5, popup = ~park_name_jp) %>% leaflet::addLegend("bottomright", pal = pal, values = ~area_id, data = parklist)
sf::st_convex_hull()の場合
まずは千代田区のポートを例に{sf}パッケージのst_convex_hull()でポリゴンを生成してみる。
# ポリゴン作成 convex <- parklist %>% dplyr::filter(area_id == target_area_id) %>% sf::st_union() %>% sf::st_convex_hull() # 可視化 map %>% leaflet::addPolygons(data = convex, color = pal(target_area_id), weight = 1.0, fillColor = "Transparent") %>% leaflet::addCircleMarkers(data = dplyr::filter(parklist, area_id == target_area_id), radius = 3, color = ~pal(area_id), weight = 1, fillColor = ~pal(area_id), fillOpacity = 0.5, popup = ~park_name_jp)
悪くはないけど少し膨らんで見える。
concaveman::concaveman()の場合
次に{concaveman}パッケージのconcaveman()で ポリゴンを生成してみる。
# ポリゴン作成 concave <- parklist %>% dplyr::filter(area_id == target_area_id) %>% concaveman::concaveman() # 可視化 map %>% leaflet::addPolygons(data = concave, color = pal(target_area_id), weight = 1.0, fillColor = "Transparent") %>% leaflet::addCircleMarkers(data = dplyr::filter(parklist, area_id == target_area_id), radius = 3, color = ~pal(area_id), weight = 1, fillColor = ~pal(area_id), fillOpacity = 0.5, popup = ~park_name_jp)
求めていたのが出来た！
concaveman::concaveman()には引数にconcavityとlength_thresholdの2つがある。
- concavity
- 相対的な凹みの尺度（デフォルトは
2）1は比較的詳細な形状になる無限大にすると凸包になる- 1未満の値にすると奇妙な形状になる
- 相対的な凹みの尺度（デフォルトは
- length_threshold
concavityを0.5, 1, 10で実行した結果は以下の通り。
0.5は一筆書きみたいになっているので巡回セールスマン問題に使えるかもしれない。
シェアサイクルの最適なポート密度は1平方キロメートルあたり10〜16とされている。*1 密度の計算に総面積や可住地面積を使用すると過小な値となってしまうため適切にサービスエリア面積を求める必要がある。上記で生成したポリゴンを元にサービスエリア面積を計算し、エリア毎にポート密度を求めてみる。
parklist %>% dplyr::group_by(area_id) %>% tidyr::nest() %>% dplyr::mutate( n = purrr::map_int(data, NROW), service_area = purrr::map_dbl(data, function(x){ x %>% concaveman::concaveman() %>% sf::st_area() %>% units::set_units(km^2) }), density = n/service_area ) %>% dplyr::ungroup() %>% dplyr::select(-data) %>% dplyr::arrange(area_id) %>% knitr::kable()
| area_id | n | service_area | density |
|---|---|---|---|
| 1 | 94 | 8.741680 | 10.753082 |
| 2 | 58 | 5.975986 | 9.705511 |
| 3 | 117 | 14.755185 | 7.929416 |
| 4 | 130 | 30.739930 | 4.229027 |
| 5 | 83 | 15.018784 | 5.526413 |
| 6 | 66 | 8.091223 | 8.156987 |
| 7 | 79 | 32.209518 | 2.452691 |
| 8 | 48 | 11.296027 | 4.249282 |
| 9 | 47 | 14.334039 | 3.278908 |
| 10 | 81 | 12.612764 | 6.422066 |
| 11 | 19 | 7.938291 | 2.393462 |
| 12 | 30 | 8.240241 | 3.640670 |
| 13 | 8 | 2.644939 | 3.024644 |
昼間人口の多い千代田区・中央区・港区・文京区はポート密度が高く、大田区・練馬区の密度が低い。人口やトリップ数も考慮する必要があるため一概には言えないがコロナ禍で密を避ける交通手段としてシェアサイクルが注目されているのでポートが充実することに期待したい。
