はじめに

NTTドコモクロステック開発部の中村と申します。本記事では、バスのオープンデータを使って、遅延が発生しやすい時間帯や場所をアニメーションで可視化してみようと思います。
※普段の業務では交通や人流データの分析やシミュレーションを行っていまして、過去の記事もあるのでぜひ覗いてみてください！

1. バスのリアルタイムロケーションデータを取得

オープンデータの公開元

公共交通オープンデータ協議会様が公共交通オープンデータセンターを運用されていて、バスをはじめ数多くの公共交通のデータが公開されています。
※データを使用する際の基本ライセンスやガイドラインはこちらから確認できます。

今回は、東京都交通局様のバスに関するデータを使わせていただこうと思うので、リンク先ページの左カラムのフィルタで、組織＝東京都交通局、タグ＝バスを指定します。すると、下記のようにバス停情報、バス路線情報...とバスに関するデータが複数件ヒットします。

バスロケーション情報が時々刻々と変わるバスの位置情報のデータ、それ以外は時刻表やバス停の緯度経度情報といった静的なマスターデータのようです。

本記事の執筆時点で、バスロケーション情報はJSON形式とProtocol Buffers形式の２種類で提供されています。恥ずかしながら、著者はProtocol Buffersなるものを執筆時に初めて知ったので、勉強を兼ねてProtocol Buffers形式で情報を取得していきます。
画面上で「バスロケーション情報」を選択し進めて行くと下記の画面にたどり着きます。

ライセンス：Creative Commons Attribution 4.0 International (CC BY 4.0)として公開頂いているので、本情報を取得、加工させて頂こうと思います。

赤枠の通り、APIにてデータが提供されているので、画面に沿ってユーザ登録しアクセストークンを取得してAPIを叩くことで定期的にデータを取得できます。本記事もこのAPIで取得することを前提にこの先進んでいきます。

(ちょっと脱線)Protocol Buffersとは？

Protocol Buffersがどんなものかについては、こちらの記事やこちらの記事にて勉強させて頂きました。
JSONとの明確な違いとしては、JSONはデータをそのまま保持するが、Protocol Buffersはバイナリ形式で記録され、データ構造を別途保持する、ということと理解しました。データのぱっと見のわかりやすさはJSONに軍配があがりますが、こちらの記事によると、Protocol Buffersの方がパース速度が圧倒的に速いらしいので、一長一短ありそうです。

バスロケーション情報をpythonで取得

前置きが長くなりましたが、オープンデータが取得方法もわかったので、いよいよバスロケーションデータをpythonで取得していきます。
バスに関するデータ仕様は「動的バス情報フォーマット(GTFSリアルタイム)」として標準化されていて、データ構造を表す「gtfs-realtime.proto」は、上記リンクからダウンロードできるのでこれで元のバスロケーション情報を取得できそうです。
※余談ですがGTFSのGは開発元であるGoogleの頭文字だったそうですが、現在は無用な誤解や忌避を避けるため、Genral Transit Feed Specificationの略ということになっているそうです（Wikipedia情報）

この「gtfs-realtime.proto」を使って情報を取得しても良いのですが、GTFSリアルタイムデータをパースできるpythonのライブラリを見つけたので、こちらを使っていきます。ライセンスは、Apache License 2.0です。

下記のようにpipでインストールできます。

pip install gtfs-realtime-bindings

まずはAPIでデータを取得してdataframeに変換しようと思います。ソースコードはこちらの記事を大いに参考にさせて頂きました。

import pandas as pd
import urllib.request, urllib.error
from google.transit import gtfs_realtime_pb2
import datetime

token = "***" #ユーザ登録で取得したアクセストークン
toei_bus_location_URL = "https://api.odpt.org/api/v4/gtfs/realtime/ToeiBus?acl:consumerKey="
feed = gtfs_realtime_pb2.FeedMessage()
API_Endpoint = toei_bus_location_URL+token

column = ["id","trip_id","route_id","direction_id","lat","lon","current_stop_sequence","timestamp","stop_id"]

def get_gtfs_rt():
    result = []
    now = datetime.datetime.utcnow().replace(microsecond=0) + datetime.timedelta(hours=9)#現在時刻を取得
    now_str = now.strftime('%Y%m%dT%H%M%S')#現在時刻を文字型に変換

    with urllib.request.urlopen(API_Endpoint) as res:
        feed.ParseFromString(res.read())
        for entity in feed.entity:
                record = [
                entity.id,                            #車両ID
                entity.vehicle.trip.trip_id,          #一意に求まるルート番号
                entity.vehicle.trip.route_id,         #路線番号（≒系統）
                entity.vehicle.trip.direction_id,     #方向（上り下り）
                entity.vehicle.position.latitude,     #車両経度
                entity.vehicle.position.longitude,    #車両緯度
                entity.vehicle.current_stop_sequence, #直近で通過した停留所の発着順序
                entity.vehicle.timestamp,             #タイムスタンプ
                entity.vehicle.stop_id,               #直近で通過した停留所
                ]
                result.append(record)

    df = pd.DataFrame(result,columns=column)
    df["timestamp"] = pd.to_datetime(df.timestamp, unit='s', utc=True).dt.tz_convert('Asia/Tokyo')  #タイムスタンプ情報をUNIX時間から日本時間に変換
    df["timestamp"] = df["timestamp"].dt.tz_localize(None)  #Timezone情報を削除

APIで取得できるタイムスタンプはUNIX時間（エポック秒）になっていて、そのままだと扱いにくいので日本時間に変換しています。

これでpythonでリアルタイムなデータを取得できるようになりました。後は、scheduleライブラリ等を使えば、一定間隔でリアルタイムデータを取得できます。

1日分のロケーション情報をdataframeにまとめたのがこちらです。以降の記事ではこのdataframeを変数dfとして扱っていきます。

日によって異なりますが、この日の運行分（2023/11/14運行分）は画像の通り300,000レコード程ありました。

trip_idとcurrent_stop_sequenceで重複するレコードを抽出したところ8,000件程に減ったので、頻繁に（数秒～数十秒おきに）位置情報が取得されているというよりは、停留所を通過した、など一定のトリガーで取得されているようです。

次に、ロケーション情報が更新されるtimestampがどの時点のものか調べるため、下記のコードで「停留所」と「バスのロケーション」との距離を算出します。そのためには、バス停の緯度経度が必要なので、先程の公共交通オープンデータセンターから静的データ（GTFSデータ）を取得します。ToeiBus-GTFS.zipというファイルをダウンロードして解凍し、その中のstop.txtに各バス停の緯度経度が記録されています。

df= df.drop_duplicates(["trip_id","current_stop_sequence"],keep='first')#同じ区間に２つ以上のレコードがあるデータを削除

df_stops =  pd.read_csv("stops.txt",usecols=["stop_id","stop_lat","stop_lon"])#停留所情報の静的データをインポート
df = df.merge(df_stops,on="stop_id",how='left')#stop_idでマージ

df["distance"] = df.apply(lambda x:geodesic((x["lat"], x["lon"]), (x["stop_lat"], x["stop_lon"])).m, axis=1)#timestampが記録された時の緯度経度とstop_idの緯度経度の距離を算出

df['distance'].hist()#距離のヒストグラムを表示

算出した距離のヒストグラムはこのようになりました。

停留所の位置とバスのロケーションはほぼ一致しているようです。（もしくは「バスのロケーション＝通過した停留所の位置」として返却している可能性もありますね。）このことからバス停を通過したタイミングでtimestampが記録され、ロケーション情報が連携されている、と考えて良さそうです。

2. 遅延が発生しやすい場所を抽出（Python）

dataframeに時刻表からの遅延を追加する

次に、バス停の時刻表上の出発時刻とバスロケーション情報の当該バス停を通過したときのtimestampを照らし合わせて、計画通りか遅延しているかを判別していきます。

先程ダウンロードしたToeiBus-GTFS.zipの中のstop_times.txtが、いわゆる時刻表に相当するのでこのファイルを使っていきます。 ※なお、stop_times.txtで記録されているdeparture_timeはHH:MM:SS形式ですが、日を跨いだレコードは下図の通り24:00:00といった値です。

このままではdatetime型として扱えないため、str2datetimeといった関数を作って一工夫しています。(こちらのQAを参考にさせて頂きました。)

def str2datetime(s,date):#24時以降の時刻を正しくdatetime型に変換する関数
    if s == s:
        H,M,S = map(int, (s[:2], s[3:5], s[6:8]))
        y=date.year
        m=date.month
        d=date.day
        if H>=24:
            H=H-24
            d=d+1
        return datetime.datetime(y,m,d) + datetime.timedelta(hours=H, minutes=M, seconds=S)
    else:
        return

df_dia =  pd.read_csv("stop_times.txt").loc[:,"trip_id":"stop_headsign"].drop(['stop_id'],axis=1)#時刻表データの取り込み
df = df.merge(df_dia, how='left',left_on=["trip_id","current_stop_sequence"],right_on=["trip_id","stop_sequence"])#trip_idとstop_sequenceで結合

date = df["timestamp"][0] #運行日を取得（11/14 25時の運行であれば、運行日は11/15ではなく11/14）
df["departure_time"] = df["departure_time"].apply(str2datetime, args=(date,))

df["late_time"] = df["timestamp"]-df["departure_time"]#timestampと出発時刻の差を計算
df["late_time"] = df["late_time"].apply(lambda x: x.seconds if (x.days>=0) else 0)#差が正の場合は秒換算した値を代入、負の場合は0を代入　#後述の可視化第一弾で使用

df["time_period"]= df["departure_time"].dt.round("H")#出発時刻を1時間単位で格納
df_g = df.groupby(['route_id','direction_id','time_period'])["trip_id"].nunique().reset_index()#各ルートが1時間当たりに何便運航しているかを集計
df_g.rename({'trip_id':'trip_count'},axis=1,inplace=True)

df = df.merge(df_g,on=['route_id','direction_id','time_period'])
df['late'] = df['late_time'].apply(lambda x:x>=300 if 1 else 0)#5分以上遅延している場合は、1（遅延）それ未満であれば0とする。
df['late_ratio']=df['late']/df['trip_count']#1時間当たりの運行本数に対する遅延便数　#後述の可視化第二弾で使用

これで、遅延時刻が記録されたdataframeを作ることができました！このdataframeをcsvとして出力し、QGISにインポートしていきます。

3. 地図上に渋滞区間を可視化

QGISにデータをインポート

QGISを起動し、データソースマネージャーを開きます。出力したcsvを選択し、図のようにX属性、Y属性、CRSが設定されていることを確認します。（もしうまく読み込まれていなければ手動で選択します。）

追加ボタンを押し、CSVデータをインポートします。またOpenStreetMapも表示させると、図のように各バスの位置情報が地図上に重畳されました。紫色のポチポチが一日分のバスのロケーション情報を表しています。東京都心部と小平市以西に運行されていることがわかりますね。

QGISで遅延箇所を可視化（第一弾）

目的：遅延時間と場所の可視化

時刻表よりも出発が遅くなったバスを可視化していきます。 csvデータのレイヤプロパティを開き、シンポロジの設定で「連続値による定義」を選択します。さらに、同じくレイヤプロパティで、時系列を選択し、図のように設定します。画面上に時系列コントローラが表示されていない場合は、上の時計のアイコンをクリックすることで表示させることができます。

遅延時間を円の大きさで表現してみました。円の中心はバスのロケーション情報です。時系列を設定することで、図のように15時台のデータのみを表示、といったことが可能となります。

15時台は錦糸町や御茶ノ水付近の円が大きいので、その辺りでは遅延が発生していたということがわかります。時系列コントローラの保存ボタンを押すと、設定した時間間隔のキャプチャ画像が生成されます。一連の画像をgifとしてつなげると下図のようになります。

時間に応じて円の大きさ（遅延時間）が変わっていって面白いですね。場所によって、朝と夕の時間帯のみ円が大きい（遅延が発生している）区間もあれば、全くと言っていいほど遅延が発生しない区間、一日を通して円が大きい区間もあるなど、場所によって様々な特性があることが分かります。地図を使って可視化すると全体像を把握するのにとても便利ですね。