PostgreSQL+PostGISに格納したGTFS形式の公共交通オープンデータを操作する

#はじめに
この記事では、PostgreSQL+PostGIS環境にインポートしたGTFSデータの簡単な操作を解説します。記事の前提として、GTFSファイルを[GTFSDB](https://github.com/OpenTransitTools/gtfsdb)というツールを用いてPostgreSQL+PostGISにインポートした環境を利用します。この方法の詳細は、過去記事である「[GTFSファイルをマージしてPostgreSQL+PostGISに投入する方法](https://qiita.com/niyalist/items/22b743d973e4a6aeef86)」、「[群馬県公共交通オープンデータ（GTFS）を一括してデータベースに格納する方法](https://qiita.com/niyalist/items/d9bf70080068a68a8587)」を参考にしてください。後者の記事からは、データベースのダンプファイルをダウンロード出来ます。

この記事は、GTFSやSQLについての基礎的な知識はあると想定しています。GTFSに関しては、仕様書の他に「[やんばる急行バスと学ぶ標準的なバス情報フォーマット(GTFS-JP)](https://qiita.com/kumatira/items/c1e179c1119c9d1c1236)」といった記事も参考になります。SQLに関しては多くの参考書が出ていますが、私自身は「[10年戦えるデータ分析入門 SQLを武器にデータ活用時代を生き抜く・青木 峰郎](https://www.amazon.co.jp/dp/4797376279/)」という書籍がわかりやすく感じました。

## なぜGTFSDBを用いるのか？
GTFS形式自体が、リレーショナルデータベースをCSVファイルへダンプした形式であるため、特別なツールを使わなくても簡単にデータベースに格納できます。それにも関わらずGTFSDBを用いるのには、以下の理由があります。もちろん、高度なアプリケーション開発を目指す場合は独自にインポート操作を行った方がいいでしょう。

### PostGISのジオメトリ型に対応しているため

PostGISには空間データを表すジオメトリ型というデータ形式があります。GTFSファイルを直接インポートしただけでは、緯度経度が単なる数値としてしか登録されませんが、GTFSDBを用いると、ジオメトリ型に変換され登録されます。`stops.txt`中の駅・バス停の緯度経度はもちろん、`shape.txt`に路線形状が登録されている場合、点の集合から線(`LINESTRING`)が作成され、`patterns`というテーブルに格納されます。これによって、地図表示や距離計算、距離など地理条件による検索が容易に出来るようになります。以下は群馬県公共交通オープンデータのバス停データと路線形状をpgAdmin4で示した例です。

![stops_and_shapes.jpg](https://qiita-image-store.s3.amazonaws.com/0/53416/be2c6a11-7494-0319-18d5-5393721c5ad2.jpeg)


### 公共交通の運行パターンに対応しているため

公共交通の運行パターンは、平日/休日/祝日により異なるのが一般的です。さらに細かく見ていくと、盆暮れ正月、特定のイベント開催日などに独自の運行パターンがあったり、地方によっては「隔週水曜日」など様々なパターンが存在します。GTFSでは`calendar.txt`、`calendar_dates.txt`という2つのファイルでこれらを表現していますが、正直、その解釈はやや面倒です。GTFSDBでは、インポート時にこれを解釈し「何月何日はどのパターンか」という情報に落とし込みます。そのため、解釈に悩むことなく容易にデータを利用出来るようになります。

# データベースに問い合わせる

GTFSデータを扱うシナリオを想定し、データベースへの問い合わせを実行します。

## 特定の地点から半径1km以内の駅・バス停を取得する

GTFSというかPostGISの入門的なお題ですが、必要な場面は多いと思います。ここでは群馬県庁（139.06083 36.39111）付近のバス停を探してみます。バス停の緯度経度、距離も出力します。

```sql
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を作成してここにアクセスしています。これを見ると、しっかり条件に合ったバス停が選択されていることが分かります。
<img width="1468" alt="スクリーンショット 2018-12-27 14.11.23.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/53416/876cb83e-997c-0ec6-c204-4c03dcd8e615.jpeg">

なおバス停の表現には、乗り場（ポール）が複数ある場合に、代表点1つで表現する方法、ポールひとつひとつを表現する方法、ポールひとつひとつを表現するうえに、代表点を置き親子関係を表現する方法があります。今回例に用いている群馬県のデータにおいても混在しているので、実際のアプリケーション開発においては、前処理の段階で形式を揃えた方がいいかもしれません。

## ある日の運行パターンを一覧する

GTFSでは、`service_id`というIDで運行パターンを表現しています。それでは、特定の日に該当する`service_id`は何でしょうか。ここでは、GTFSDBがインポート時に追加した`universal_calendar`というテーブルを参照します。`universal_calendar`は、簡単な構造でservice_idと日付(date)のペアの一覧を格納したテーブルになっているので、日付を検索条件にして結果を求めればいいということが分かります。そのため、例えば12月26日に運行する `service_id` の一覧は以下のように求められます。

```sql
select 
	service_id
from
	universal_calendar
where
	date='2018-12-26'
```

## ある日運行されるtripを一覧する

日本語では「便（びん）」、つまり「海街線15:26発 市民ホール行き」というような、運行の1単位をGTFSではtripと表現します。tripは、始発駅（バス停）から終着駅（バス停）まで特定の時刻で走ることが1単位であり、毎日運行されるダイヤなら、翌日また同じtripが運行されることになります。実際には、平日だけの便、休日だけの便、特定の日だけ走る臨時便などがあるので、ここでは、特定の運行日から運行される便の一覧を検索します。

ここでポイントは、`service_id`というフィールドに対してある検索条件で絞り込むことです。whereの条件中に、先ほど書いた「特定の日の運行パターン一覧を求めるSQL」をそのまま突っ込んでしまいましょう（サブクエリといいます）。

```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）し、その結果を日付ごと集計します。なお、期間を年末年始の期間に限定しています。

```sql
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
```

ちなみに、群馬県オープンデータでこれを実行すると以下のような結果になります。やはり元日がいちばん少なく、思ったより、日々の変動があります。

![number_of_trips.png](https://qiita-image-store.s3.amazonaws.com/0/53416/ba83cc34-82c9-c1bd-e00f-9e5c2dfbeb9d.png)

## 特定のバス停の時刻表を取得する

バス停のIDが分かったとして、そのバス停を通過するバスの時刻表を取得します。ここでは、路線名、行き先も取得したかったので、`routes`、`trips`テーブルもjoinしています。注意する必要があるのは、stop_headsignが設定されている場合、trip_headsignを上書きするため、CASE式を用いてどちらかをheadsignとして出力するようにしています。`service_id`に関しては、先ほどと同じサブクエリの考え方です。

```sql
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`を変えて結果を確認してみてください。
<img width="781" alt="スクリーンショット 2018-12-27 17.46.06.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/53416/c53249fa-4c81-e22a-ab92-70074bed1709.jpeg">



## 近傍のバス停の時刻を束ねて出力する

群馬県のデータには、「県庁前」という名の付いたバス停が10箇所（親バス停を除く）あります。実際には、これをまとめた時刻表が欲しくなります。考え方としては、位置で検索する方法と名称で検索する方法がありますが、実際にバスを利用する状況を考えると、位置で検索した方がわかりやすいのではないでしょうか。そこで、「特定の地点から半径1km以内のバス停を取得する」を思いだし、where節のサブクエリに設定します。今回は、距離を半径250mにしましょう。

```sql
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
```

この検索結果は以下のようになります。「県庁前」バス停の他に「市役所・合庁前」が含まれていますね。

<img width="1300" alt="スクリーンショット 2018-12-27 17.56.24.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/53416/9dd09edf-b190-bbfa-754d-067e65089aa9.jpeg">
これは、地図を見れば明らかです。県庁の位置がちょっとずれてる気もしますが、実際、県庁から「市役所・合庁前」も歩いてすぐなので、これで良しとします。

<img width="1468" alt="スクリーンショット 2018-12-27 18.00.12.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/53416/85ae54c3-bc07-24a7-3925-d3a6f0e1240a.jpeg">

## ある場所からある場所へのバスの一覧を取得する

2地点間を直接結ぶバスの情報を取得します。具体的には前橋駅（139.07326 36.38328）から群馬県庁（139.06083 36.39111）へのバスを例に考えますが、この区間は、実際多くの会社の多くのバスが走っており、とても分かりにくい状況です。この区間の時刻表情報を取得しましょう。

基本的な考え方としては、「前橋駅付近を出発するバス（origin）」と「県庁付近に到着するバス（destination）」をそれぞれ検索し、内部結合（inner join）することで両方を通るバスを抽出します。SQLのwith句を使って可読性を高めています。

ただこれだけでは、250m以内に複数のバス停があった場合に複数の時刻が出てしまいます。そこで、Window関数のrank()を用いて、同一trip_idのデータをdistanceでランク付けして、最も近いバス停だけを選ぶようにしています。厳密にはこの方法では、循環路線があった場合に正しい結果が得られない可能性がある（図を参照）のですが、そこは、SQLで処理するものではない気もして、このままにしています。（引き算で乗車時間を計算してもいいですが、24時以降の時刻表現がありうることに注意してください）

![無題1.jpg](https://qiita-image-store.s3.amazonaws.com/0/53416/654fb491-7002-3bac-e13a-9c335ac44700.jpeg)


```sql

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本もバスを走らせるすごい区間です。かなり本数が多いのと、同じ時刻で発車するバスもあれば、長く空く時間もあるのが気になりますね。

<img width="1200" alt="スクリーンショット 2018-12-27 22.50.15.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/53416/142a88aa-a9d8-e14d-9d48-52b6bd840acb.jpeg">

<img width="271" alt="スクリーンショット 2018-12-27 22.53.33.jpg" src="https://qiita-image-store.s3.amazonaws.com/0/53416/89be5338-1158-e13e-8ff4-c6a9941a869d.jpeg">


さてこの結果、正しいでしょうか？実は[群馬県庁のサイトに前橋駅から群馬県庁へのバスの時刻一覧が掲載](http://www.pref.gunma.jp/04/h2110131.html)されてて、一致することを確認しました。せっかく整備したオープンデータですので、是非まずは、県庁のWeb担当者の方に使って頂きたいです。このSQLの出力結果を元にデータを成形すれば、だいぶ業務が楽になるのではないでしょうか。