標準的な公共交通機関の情報形式 GTFS をRで処理する: gtfsrouter編
この記事はRアドベントカレンダー2019の5日目の投稿です。昨日は id:ando_Roid さんの投稿でした。
はじめに
今回で4回目となる「東京公共交通オープンデータチャレンジ」の応募が始まりました。これは公共交通オープンデータやその他のオープンデータを活用したアプリケーションおよびサービスの提案に関するコンテストです。ここでは、鉄道、バス、航空などの交通機関のデータが、GTFS (General Transit Feed Specification) 形式で公開されています。
GTFSは公共交通機関の時刻表と地理的情報に関するオープンフォーマットとして定義されます。リアルタイムの運行情報を反映した動的データ (GTFS Realtime)と、ダイヤ改正等の事情がない限り安定した静的データがあり、GTFSといえば厳密には静的データの方を指します。複数の事業者・交通機関が共通の形式に従ってデータを整備することで、経路探索や時刻表のアプリケーション開発が容易になるのが利点です。そのため、GTFSを提案したGoogle以外でも広く使われています。日本の交通機関に関しても、GTFSをベースに日本独自の交通事情を加味しながら拡張されたフォーマット (GTFS-JP) が整備されています。
さて、今回はそんなGTFSをRで扱ってみようという記事になります。日本語の内容として id:nonki1974 さんがgtfsrパッケージの利用方法を スライドにまとめています。なのでここでは、もう一つの例として gtfsrouterパッケージを紹介します。
目次です。
gtfsrouterパッケージ
gtfsrouterパッケージの特徴として次の点が挙げられます。
- 2019年12月5日現在、CRANに登録されている。gtfsrは未登録
- GTFSデータをdata.tableオブジェクトとして処理するため、処理速度の高速化が期待できる
- 指定時間内に移動可能な領域 (Isochrone map) の描画ができる
まずはパッケージを利用可能な状態にしておきます。CRANに登録されているバージョンではなくGitHub上の開発版をインストールします。
install.packages("remotes") remotes::install_github("ATFutures/gtfs-router") # gtfs_isochrone()を実行する際に必要です # install.packages(c("alphahull", "geodist"))
library(magrittr) library(gtfsrouter) library(data.table)
ここではデモデータとして用意されているベルリン、ブランデンブルクを運行する交通機関と、日本の岡山電気軌道株式会社が運営する岡電バス http://www.okayama-kido.co.jp/bus/ のGTFSデータを処理する例を見ていきます。なお岡電バスのデータは冒頭の「東京公共交通オープンデータチャレンジ」で提供されているものではなく、元々Creative Commonsライセンス4.0のもとで公開されているデータとなります。
データ公開元: https://loc.bus-vision.jp/ryobi/view/opendata.html
デモデータの処理
デモデータと岡電バスのデータを分けて解説するのは、現行のCRANバージョンでは、transfers.txtが存在しないGTFSデータを読み込むことができないためです。またtransfers.txtがない場合に実行できないいくつかの関数がありますので留意してください (issue #14 参照)。
berlin_gtfs_to_zip() tempfiles <- list.files(tempdir (), full.names = TRUE) filename <- tempfiles [grep ("vbb.zip", tempfiles)]
ここまでがデータを用意するところです。次のextract_gtfs()でGTFSの圧縮ファイル（運行情報に関する各種のcsvファイルが含まれる）を読み込みます。
gtfs <- extract_gtfs(filename)
データはリストに格納されています。どのような情報が格納されているかは、names()を使って調べます。また、各データごとの件数が知りたい時はsummary()を使うと良いでしょう。
names(gtfs)
## [1] "calendar" "routes" "trips" "stop_times" "stops"
## [6] "transfers"
summary(gtfs)
## A gtfs object with the following tables and respective numbers of entries in each:
## calendar routes trips stop_times stops transfers
## 2395 52 28590 22666 957 9833
読み込んだデータの一例として便情報 trips
を表示します。この出力で、各情報はdata.tableオブジェクトであることが確認できます。
head(gtfs$trips)
## route_id service_id trip_id trip_headsign trip_short_name direction_id block_id shape_id wheelchair_accessible bikes_allowed
## 1: 10141_109 139 107928371 S Blankenburg NA 0 NA 1100 1 1
## 2: 10141_109 139 107928370 S Blankenburg NA 0 NA 1100 1 1
## 3: 10141_109 139 107928369 S Blankenburg NA 0 NA 1100 1 1
## 4: 10141_109 139 107928368 S Blankenburg NA 0 NA 1100 1 1
## 5: 10141_109 139 107928367 S Blankenburg NA 0 NA 1100 1 1
## 6: 10141_109 139 107928366 S Blankenburg NA 0 NA 1100 1 1
class(gtfs$trips)
## [1] "data.table" "data.frame"
# 起点と終点、時刻を選んでデータを抽出 gtfs_route(gtfs, from = "Schonlein", to = "Berlin Hauptbahnhof", # 12:02 in seconds start_time = 12 * 3600 + 120, day = "Sunday")
## route_name trip_name stop_name arrival_time departure_time
## 1 U8 S+U Wittenau U Schonleinstr. (Berlin) 12:09:00 12:09:00
## 2 U8 S+U Wittenau U Kottbusser Tor (Berlin) 12:11:00 12:11:00
## 3 U8 S+U Wittenau U Moritzplatz (Berlin) 12:13:00 12:13:00
## 4 U8 S+U Wittenau U Heinrich-Heine-Str. (Berlin) 12:14:30 12:14:30
## 5 U8 S+U Wittenau S+U Jannowitzbrucke (Berlin) 12:15:30 12:15:30
## 6 S3 S Spandau Bhf S+U Jannowitzbrucke (Berlin) 12:20:54 12:21:24
## 7 S3 S Spandau Bhf S+U Alexanderplatz Bhf (Berlin) 12:22:54 12:23:42
## 8 S3 S Spandau Bhf S Hackescher Markt (Berlin) 12:24:54 12:25:24
## 9 S3 S Spandau Bhf S+U Friedrichstr. Bhf (Berlin) 12:26:54 12:27:42
## 10 S3 S Spandau Bhf S+U Berlin Hauptbahnhof 12:29:36 12:30:12
あらかじめ特定の日付や曜日、ルートのパターンを限定しておくことでgtfs_route()の処理を高速化が期待できます。
gtfs_timetable()でデータの抽出を行います。
# 出力結果は省略 gtfs_sunday <- gtfs_timetable(gtfs, day = "Sunday") gtfs_sunday %>% gtfs_route(from = "Schonlein", to = "Berlin Hauptbahnhof", start_time = 12 * 3600 + 120)
data.tableによるGTFSデータの操作
次に岡電バスのGTFSデータを見てみます。ここでは先ほどのgtfs_route()やgtfs_timetable()を使わずに任意のGTFSデータを紐付ける処理を解説します。先述の通り、extract_gtfs()の読み込み結果はdata.tableオブジェクトを格納したリストですので、このデータを処理するためにはdata.tableパッケージの関数についての理解が求められます(dplyrパッケージに慣れている場合はas.data.frame()としても良いです）。
zipファイルをダウンロードすることになりますが、このファイルは展開せずにそのまま利用します。
download.file(url = "http://loc.bus-vision.jp/gtfs/okaden/gtfsFeed", destfile = "data/okaden_GTFS-JP.zip") # ここで警告が出ますが問題はありません gtfs <- extract_gtfs("data/okaden_GTFS-JP.zip")
Warning message:
This feed contains no transfers.txt
まずは便情報 trips
から対象のデータを絞っていきます。tripsには目的地の停留所名があるので、これを元にします。またroute_id、trip_idの情報が他のデータと紐付けるために必要となります。
# 目的地が「岡山駅」の便を検索 # route_id ... ルートを識別するためのID # trip_id... 旅程を識別するID # trip_headsign... 目的地 dt_trips <- gtfs$trips[trip_headsign == "岡山駅"][, list(route_id, trip_id, trip_headsign, shape_id)]
「岡山駅」が目的地となるルートは複数あるので、1つのルートのデータを使うことにします。
dt_trips <- dt_trips[route_id == "10001_1008_1"] %>% unique(by = "route_id") dt_trips
## route_id trip_id trip_headsign shape_id
## 1: 10001_1008_1 347_1114411_20191201 岡山駅 2_1008_1
次にルート routes から路線名を調べます。
gtfs$routes[route_id == dt_trips$route_id][, list(route_id, route_long_name)]
## route_id route_long_name
## 1: 10001_1008_1 大学病院線
続いて、停車地の情報を得れるようにします。stop_times が
trip_idの情報と紐づきます。
dt_stopstimes <- gtfs$stop_times[trip_id == dt_trips$trip_id][, list(trip_id, stop_id)] dt_stopstimes
## trip_id stop_id
## 1: 347_1114411_20191201 18_5
## 2: 347_1114411_20191201 37_1
## 3: 347_1114411_20191201 38_1
## 4: 347_1114411_20191201 39_3
## 5: 347_1114411_20191201 40_1
## 6: 347_1114411_20191201 138_1
## 7: 347_1114411_20191201 63_4
## 8: 347_1114411_20191201 64_1
しかしこれだけではバス停の位置や名前の情報がありません。さらに 停車地
stops のデータを参照する必要が生じます。
dt_stops <- gtfs$stops[stop_id %in% dt_stopstimes$stop_id][, list(stop_id, stop_name, stop_lat, stop_lon)] dt_stops
## stop_id stop_name stop_lat stop_lon
## 1: 18_5 岡山駅 34.665280 133.91846
## 2: 37_1 イオンモール岡山前・源吉兆庵本社前 34.661926 133.91895
## 3: 38_1 山陽新聞社前 34.657890 133.91911
## 4: 39_3 岡山市役所前 34.655470 133.91866
## 5: 40_1 水道局前 34.653824 133.91756
## 6: 63_4 清輝橋 34.651443 133.92488
## 7: 64_1 岡南営業所 34.644432 133.92485
## 8: 138_1 大学病院入口 34.652370 133.92053
これで「大学病院線」の停車地の情報が用意できました。
今度は車両が通過する経路を作りましょう。経路の情報は shapes
にあり、便情報とはshape_idで繋がります。
dt_shapes <- gtfs$shapes[shape_id == dt_trips$shape_id][, list(shape_pt_lon, shape_pt_lat)] dt_shapes
## shape_pt_lon shape_pt_lat
## 1: 133.924891667 34.6442527778
## 2: 133.925005556 34.6464027778
## 3: 133.925002778 34.6468666667
## 4: 133.925036111 34.6476444444
## 5: 133.925097222 34.6482138889
## 6: 133.925211111 34.6504444444
## 7: 133.925475000 34.6512555556
## 8: 133.925400000 34.6513361111
## 9: 133.924158333 34.6516277778
## 10: 133.922086111 34.6520777778
## 11: 133.917530556 34.6531000000
## 12: 133.917425000 34.6532277778
## 13: 133.917550000 34.6536944444
## 14: 133.917827778 34.6540611111
## 15: 133.919227778 34.6562055556
## 16: 133.919277778 34.6564500000
## 17: 133.918972222 34.6630277778
## 18: 133.918952778 34.6640222222
## 19: 133.918997222 34.6646805556
## 20: 133.919166667 34.6654277778
## 21: 133.919100000 34.6655166667
## 22: 133.918891667 34.6654416667
## 23: 133.918411111 34.6647361111
## 24: 133.918325000 34.6647111111
## 25: 133.918236111 34.6647638889
## 26: 133.918277778 34.6648972222
## 27: 133.918613889 34.6653805556
## shape_pt_lon shape_pt_lat
もう一つ、曜日をベースとして停留地を探す処理の例を示します。
trip_ids <- gtfs$trips[.(gtfs$calendar[.(1), on = c("monday")]), on = "service_id"][, trip_id] %>% unique() dt_okaden_monday <- merge(gtfs$stop_times[.(trip_ids), on = .(trip_id)] %>% .[, -c("stop_headsign", "pickup_type", "drop_off_type")], gtfs$stops[, c("stop_id", "stop_name", "stop_lat", "stop_lon")], by = "stop_id") # 並び替え setorder(dt_okaden_monday, trip_id, stop_sequence) dt_okaden_monday
## stop_id trip_id arrival_time departure_time stop_sequence shape_dist_traveled stop_name stop_lat stop_lon
## 1: 64_1 347_1100111_20191201 54540 54540 1 NA 岡南営業所 34.644432 133.92485
## 2: 63_4 347_1100111_20191201 54660 54660 2 NA 清輝橋 34.651443 133.92488
## 3: 138_1 347_1100111_20191201 54720 54720 3 NA 大学病院入口 34.652370 133.92053
## 4: 40_1 347_1100111_20191201 54840 54840 4 NA 水道局前 34.653824 133.91756
## 5: 39_3 347_1100111_20191201 54900 54900 5 NA 岡山市役所前 34.655470 133.91866
## ---
## 32714: 37_1 355_1403536_20191201 67980 67980 29 NA イオンモール岡山前・源吉兆庵本社前 34.661926 133.91895
## 32715: 18_6 355_1403536_20191201 68040 68040 30 NA 岡山駅 34.665115 133.91835
## 32716: 19_5 355_1403536_20191201 68100 68100 31 NA 岡山駅前 34.665688 133.92114
## 32717: 22_1 355_1403536_20191201 68280 68280 32 NA ＮＴＴ岡山前 34.663220 133.92645
## 32718: 27_8 355_1403536_20191201 68640 68640 33 NA 天満屋 34.661713 133.92828
ややこしさもありますが、理解すれば使いやすいデータだと思います。
RでGTFSデータを扱う際のTips
ちょっとしたtipsです。
連続値になっている時刻の変換
通過時刻 stop_times などに記録される arrival_time
は実数になっています。これだと時刻がわかりにくいので、HH:MM:SS形式に変換しておくと確認が楽になります。これにはhms::as_hms()が利用できます。
min(dt_okaden_monday$arrival_time)
## [1] 20640
# 時刻表データに記録される最初と最後の時刻を表示 hms::as_hms(range(dt_okaden_monday$arrival_time))
## 05:44:00
## 22:58:00
sfオブジェクトへの変換
GTFSには、バスの停留所やルートを構成するポイントの情報が記録されています。gtfsrouterやgtfsrパッケージで読み込んだデータも、こうした緯度経度の情報はデータフレームの列に格納されています。そのため、地図上にマッピングしたり地理的な解析を行う際は地理空間オブジェクトへの変換が必要となります。ここでは、先ほどdata.tableを操作して得たデータを元にRの地理空間オブジェクトとして普及しているsfへの変換を試みます。POINTデータとして停留所、LINESTRINGの例として路線をそれぞれ扱います。
library(sf)
sf_okaden_stops <- dt_stops %>% st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326) sf_okaden_stops
## Simple feature collection with 8 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 133.91756 ymin: 34.644432 xmax: 133.92488 ymax: 34.66528
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## stop_id stop_name geometry
## 1 18_5 岡山駅 POINT (133.91846 34.66528)
## 2 37_1 イオンモール岡山前・源吉兆庵本社前 POINT (133.91895 34.661926)
## 3 38_1 山陽新聞社前 POINT (133.91911 34.65789)
## 4 39_3 岡山市役所前 POINT (133.91866 34.65547)
## 5 40_1 水道局前 POINT (133.91756 34.653824)
## 6 63_4 清輝橋 POINT (133.92488 34.651443)
## 7 64_1 岡南営業所 POINT (133.92485 34.644432)
## 8 138_1 大学病院入口 POINT (133.92053 34.65237)
sf_okaden_routes <- dt_shapes %>% as.matrix() %>% st_linestring() %>% st_sfc(crs = 4326)
geometryだけの情報になりますが、sfcオブジェクトとして扱いたい場合があれば
sfheadersパッケージの関数を使うのも一つの手段となります。
dt_stops %>% sfheaders::sfc_point(x = "stop_lon", y = "stop_lat")
## Geometry set for 8 features
## geometry type: POINT
## dimension: XY
## bbox: xmin: 133.91756 ymin: 34.644432 xmax: 133.92488 ymax: 34.66528
## epsg (SRID): NA
## proj4string: NA
## First 5 geometries:
## POINT (133.91846 34.66528)
## POINT (133.91895 34.661926)
## POINT (133.91911 34.65789)
## POINT (133.91866 34.65547)
## POINT (133.91756 34.653824)
最後にleaflet上でマッピングしてみましょう。
mapview::mapview(sf_okaden_stops) + mapview::mapview(sf_okaden_routes)
明日は id:upura さんの予定です。そちらも楽しみですね。 Enjoy!
