はじめに
この記事では、PostgreSQL+PostGIS環境にインポートしたGTFSデータの簡単な操作を解説します。記事の前提として、GTFSファイルをGTFSDBというツールを用いてPostgreSQL+PostGISにインポートした環境を利用します。この方法の詳細は、過去記事である「GTFSファイルをマージしてPostgreSQL+PostGISに投入する方法」、「群馬県公共交通オープンデータ（GTFS）を一括してデータベースに格納する方法」を参考にしてください。後者の記事からは、データベースのダンプファイルをダウンロード出来ます。
この記事は、GTFSやSQLについての基礎的な知識はあると想定しています。GTFSに関しては、仕様書の他に「やんばる急行バスと学ぶ標準的なバス情報フォーマット(GTFS-JP)」といった記事も参考になります。SQLに関しては多くの参考書が出ていますが、私自身は「10年戦えるデータ分析入門 SQLを武器にデータ活用時代を生き抜く・青木 峰郎」という書籍がわかりやすく感じました。
なぜGTFSDBを用いるのか？
GTFS形式自体が、リレーショナルデータベースをCSVファイルへダンプした形式であるため、特別なツールを使わなくても簡単にデータベースに格納できます。それにも関わらずGTFSDBを用いるのには、以下の理由があります。もちろん、高度なアプリケーション開発を目指す場合は独自にインポート操作を行った方がいいでしょう。
PostGISのジオメトリ型に対応しているため
PostGISには空間データを表すジオメトリ型というデータ形式があります。GTFSファイルを直接インポートしただけでは、緯度経度が単なる数値としてしか登録されませんが、GTFSDBを用いると、ジオメトリ型に変換され登録されます。stops.txt中の駅・バス停の緯度経度はもちろん、shape.txtに路線形状が登録されている場合、点の集合から線(LINESTRING)が作成され、patternsというテーブルに格納されます。これによって、地図表示や距離計算、距離など地理条件による検索が容易に出来るようになります。以下は群馬県公共交通オープンデータのバス停データと路線形状をpgAdmin4で示した例です。
公共交通の運行パターンに対応しているため
公共交通の運行パターンは、平日/休日/祝日により異なるのが一般的です。さらに細かく見ていくと、盆暮れ正月、特定のイベント開催日などに独自の運行パターンがあったり、地方によっては「隔週水曜日」など様々なパターンが存在します。GTFSではcalendar.txt、calendar_dates.txtという2つのファイルでこれらを表現していますが、正直、その解釈はやや面倒です。GTFSDBでは、インポート時にこれを解釈し「何月何日はどのパターンか」という情報に落とし込みます。そのため、解釈に悩むことなく容易にデータを利用出来るようになります。
データベースに問い合わせる
GTFSデータを扱うシナリオを想定し、データベースへの問い合わせを実行します。
特定の地点から半径1km以内の駅・バス停を取得する
GTFSというかPostGISの入門的なお題ですが、必要な場面は多いと思います。ここでは群馬県庁（139.06083 36.39111）付近のバス停を探してみます。バス停の緯度経度、距離も出力します。
select
stop_id,
stop_name,
ST_AsText(geom), --geometry型を文字列に変換して出力。プログラムから処理する場合は単に geom でいい
ST_Distance(ST_GeographyFromText('SRID=4326;POINT(139.06083 36.39111)'), geom) as distance,
stop_desc,
zone_id,
location_type,
parent_station
from
stops
where
ST_DWithin(geom, ST_GeographyFromText('SRID=4326;POINT(139.06083 36.39111)'), 1000.0)
order by
distance
この結果を簡単に確認するためには、QGISを用います。群馬県庁の座標を読み込み、「バッファリング」で半径1kmの円を描画し、そこに全てのバス停データ、色を変えて選択されたバス停データを重ねます。QGISはPostgreSQLに直接アクセスし、データを表示出来ます。前述の結果は、VIEWを作成してここにアクセスしています。これを見ると、しっかり条件に合ったバス停が選択されていることが分かります。
なおバス停の表現には、乗り場（ポール）が複数ある場合に、代表点1つで表現する方法、ポールひとつひとつを表現する方法、ポールひとつひとつを表現するうえに、代表点を置き親子関係を表現する方法があります。今回例に用いている群馬県のデータにおいても混在しているので、実際のアプリケーション開発においては、前処理の段階で形式を揃えた方がいいかもしれません。
ある日の運行パターンを一覧する
GTFSでは、service_idというIDで運行パターンを表現しています。それでは、特定の日に該当するservice_idは何でしょうか。ここでは、GTFSDBがインポート時に追加したuniversal_calendarというテーブルを参照します。universal_calendarは、簡単な構造でservice_idと日付(date)のペアの一覧を格納したテーブルになっているので、日付を検索条件にして結果を求めればいいということが分かります。そのため、例えば12月26日に運行する service_id の一覧は以下のように求められます。
select
service_id
from
universal_calendar
where
date='2018-12-26'
ある日運行されるtripを一覧する
日本語では「便（びん）」、つまり「海街線15:26発 市民ホール行き」というような、運行の1単位をGTFSではtripと表現します。tripは、始発駅（バス停）から終着駅（バス停）まで特定の時刻で走ることが1単位であり、毎日運行されるダイヤなら、翌日また同じtripが運行されることになります。実際には、平日だけの便、休日だけの便、特定の日だけ走る臨時便などがあるので、ここでは、特定の運行日から運行される便の一覧を検索します。
ここでポイントは、service_idというフィールドに対してある検索条件で絞り込むことです。whereの条件中に、先ほど書いた「特定の日の運行パターン一覧を求めるSQL」をそのまま突っ込んでしまいましょう（サブクエリといいます）。
select
*
from
trips
where
service_id in
(
select
service_id
from
universal_calendar
where
date='2018-12-26'
)
1日の便数の変化を調べる
1日に走る便が分かると、この数が日々どう変化するか知りたくなりますね。特に年末年始を挟んだ時期の変化は気になるところです。これはSQLによるデータ分析の入門のような課題ですね。ここでは、universal_calendarテーブルとtripsテーブルを結合（inner join）し、その結果を日付ごと集計します。なお、期間を年末年始の期間に限定しています。
select
uc.date,
count(*)
from
universal_calendar as uc
inner join trips as trips
on uc.service_id = trips.service_id
where
uc.date >= '2018-12-17' and
uc.date <= '2019-01-20'
group by
uc.date
order by
uc.date
ちなみに、群馬県オープンデータでこれを実行すると以下のような結果になります。やはり元日がいちばん少なく、思ったより、日々の変動があります。
特定のバス停の時刻表を取得する
バス停のIDが分かったとして、そのバス停を通過するバスの時刻表を取得します。ここでは、路線名、行き先も取得したかったので、routes、tripsテーブルもjoinしています。注意する必要があるのは、stop_headsignが設定されている場合、trip_headsignを上書きするため、CASE式を用いてどちらかをheadsignとして出力するようにしています。service_idに関しては、先ほどと同じサブクエリの考え方です。
select
trips.trip_id,
routes.route_id,
routes.route_short_name,
routes.route_long_name,
departure_time,
case
when stop_headsign is not null then stop_headsign
else trip_headsign
end as headsign
from
stop_times
inner join trips
on stop_times.trip_id = trips.trip_id
inner join routes
on trips.route_id = routes.route_id
inner join stops
on stop_times.stop_id = stops.stop_id
where
stops.stop_id = '135_02' and
service_id in
(
select
service_id
from
universal_calendar
where
date='2018-12-26'
)
order by
routes.route_id,
departure_time
結果は、以下のように得られます。出力結果はroute_idでまずソートし、次にdeparture_timesでソートしていますが、GTFS-JPの特徴として、route（路線）を細かく分割してデータを作成しているので、ほぼ同一の路線も別扱いとなり、やや見にくい結果になったかも知れません。stop_idを変えて結果を確認してみてください。
近傍のバス停の時刻を束ねて出力する
群馬県のデータには、「県庁前」という名の付いたバス停が10箇所（親バス停を除く）あります。実際には、これをまとめた時刻表が欲しくなります。考え方としては、位置で検索する方法と名称で検索する方法がありますが、実際にバスを利用する状況を考えると、位置で検索した方がわかりやすいのではないでしょうか。そこで、「特定の地点から半径1km以内のバス停を取得する」を思いだし、where節のサブクエリに設定します。今回は、距離を半径250mにしましょう。
select
stops.stop_id,
trips.trip_id,
routes.route_id,
stops.stop_name,
routes.route_short_name,
routes.route_long_name,
departure_time,
case
when stop_headsign is not null then stop_headsign
else trip_headsign
end as headsign
from
stop_times
inner join trips
on stop_times.trip_id = trips.trip_id
inner join routes
on trips.route_id = routes.route_id
inner join stops
on stop_times.stop_id = stops.stop_id
where
stops.stop_id in
(
select
stop_id
from
stops
where
ST_DWithin(geom, ST_GeographyFromText('SRID=4326;POINT(139.06083 36.39111)'), 250.0)
)
and
service_id in
(
select
service_id
from
universal_calendar
where
date='2018-12-26'
)
order by
stops.stop_id,
routes.route_id,
departure_time
この検索結果は以下のようになります。「県庁前」バス停の他に「市役所・合庁前」が含まれていますね。
これは、地図を見れば明らかです。県庁の位置がちょっとずれてる気もしますが、実際、県庁から「市役所・合庁前」も歩いてすぐなので、これで良しとします。
ある場所からある場所へのバスの一覧を取得する
2地点間を直接結ぶバスの情報を取得します。具体的には前橋駅（139.07326 36.38328）から群馬県庁（139.06083 36.39111）へのバスを例に考えますが、この区間は、実際多くの会社の多くのバスが走っており、とても分かりにくい状況です。この区間の時刻表情報を取得しましょう。
基本的な考え方としては、「前橋駅付近を出発するバス（origin）」と「県庁付近に到着するバス（destination）」をそれぞれ検索し、内部結合（inner join）することで両方を通るバスを抽出します。SQLのwith句を使って可読性を高めています。
ただこれだけでは、250m以内に複数のバス停があった場合に複数の時刻が出てしまいます。そこで、Window関数のrank()を用いて、同一trip_idのデータをdistanceでランク付けして、最も近いバス停だけを選ぶようにしています。厳密にはこの方法では、循環路線があった場合に正しい結果が得られない可能性がある（図を参照）のですが、そこは、SQLで処理するものではない気もして、このままにしています。（引き算で乗車時間を計算してもいいですが、24時以降の時刻表現がありうることに注意してください）
with origin as
(
select
*,
ST_Distance(ST_GeographyFromText('SRID=4326;POINT(139.07326 36.38328)'), stops.geom) as distance,
trips.trip_id as origin_trip_id
from
stop_times
inner join stops on stop_times.stop_id = stops.stop_id
inner join trips on trips.trip_id = stop_times.trip_id
inner join routes on routes.route_id = trips.route_id
inner join agency on routes.agency_id = agency.agency_id
where
stops.stop_id in
(
select
stop_id
from
stops
where
ST_DWithin(geom, ST_GeographyFromText('SRID=4326;POINT(139.07326 36.38328)'), 250.0) --前橋駅
)
),
destination as
(
select
stop_times.*,
stops.*,
ST_Distance(ST_GeographyFromText('SRID=4326;POINT(139.06083 36.39111)'), stops.geom) as distance
from
stop_times
inner join stops on stop_times.stop_id = stops.stop_id
where
stops.stop_id in
(
select
stop_id
from
stops
where
ST_DWithin(geom, ST_GeographyFromText('SRID=4326;POINT(139.06083 36.39111)'), 250.0) --群馬県庁
)
)
select
*
from
(
select
agency_name,
route_long_name,
case
when origin.stop_headsign is not null then origin.stop_headsign
else origin.trip_headsign
end as headsign,
origin.stop_name as origin_stop_name,
origin.departure_time as origin_departure_time,
destination.stop_name as destination_stop_name,
destination.departure_time as destination_departure_time,
rank() over (partition by origin_trip_id order by origin_trip_id, origin.distance, destination.distance) as distance_rank
from
origin
inner join destination on origin.origin_trip_id = destination.trip_id
where
origin.departure_time < destination.departure_time
and
service_id in
(
select
service_id
from
universal_calendar
where
date='2018-12-26'
)
) as timetable
where
timetable.distance_rank = 1
order by
origin_departure_time
出力結果は、こんな感じになります。7社で計150本もバスを走らせるすごい区間です。かなり本数が多いのと、同じ時刻で発車するバスもあれば、長く空く時間もあるのが気になりますね。
さてこの結果、正しいでしょうか？実は群馬県庁のサイトに前橋駅から群馬県庁へのバスの時刻一覧が掲載されてて、一致することを確認しました。せっかく整備したオープンデータですので、是非まずは、県庁のWeb担当者の方に使って頂きたいです。このSQLの出力結果を元にデータを成形すれば、だいぶ業務が楽になるのではないでしょうか。