[会津若松市オープンデータ]車両のプローブデータの可視化

はじめに

車両からとれるGPSや速度、加速度の時系列データ（プローブデータ）の可視化方法についてご紹介いたします。

プローブデータを集計、解析することで、
車両がどこでどんな挙動をしているかを把握することができ、
ヒヤリハットや速度違反の発生個所の確認することができます。

今回は、特定の運行のプローブデータの可視化をJupyterNotebookを使用して実行します。

データのサンプルは、
会津若松市で公開されているパトロールカーやバスのプローブデータを使用します。

目次

1. 会津若松市オープンデータ
2. データの前処理
3. 時系列プロット
4. 地図プロット
5. おわりに

会津若松市オープンデータ

会津若松市オープンデータ活用実証事業において収集した、道路パトロール車および市内循環バス「エコろん号」の走行情報履歴データを利用します。

データとしては、時系列で以下のレコードが記録されています。
測定車両名（car_name）
測定ID（measurement_data_id）
測定日時（measurement_datetime）
UNIX時間（measurement_ms）
緯度（latitude）
経度（longitude）
測位誤差（gps_error_meter）
左右加速度（accel_x_transverse）
前後加速度（accel_y_longitudinal）
上下加速度（accel_z_vertical)

データの前処理

必要なライブラリをインポートします

import pandas as pd
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.express as px
import folium
from folium.plugins import TimestampedGeoJson

上記サイトからファイルをダウンロードするとtxt形式のTSVになっているので、
TSVに合わせてデータフレームに変換します

# TSVの読み込み
df_trip = pd.read_csv('./data/20210330_151839_aizu_BL-01_20.txt', sep='\t')
df_trip.head()

UNIX時間をdatetimeに変換

# UNIX時間を変換
df_trip['datetime'] = pd.to_datetime(df_trip['measurement_ms'], unit='ms', utc=True).dt.tz_convert('Asia/Tokyo')

方位のカラムがないので、一つ前のレコードの緯度経度から向きを計算して方位を作成します。

def calculate_bearing(lat1, lon1, lat2, lon2):
    lat1 = np.radians(lat1)
    lon1 = np.radians(lon1)
    lat2 = np.radians(lat2)
    lon2 = np.radians(lon2)
    
    delta_lon = lon2 - lon1
    x = np.sin(delta_lon) * np.cos(lat2)
    y = np.cos(lat1) * np.sin(lat2) - (np.sin(lat1) * np.cos(lat2) * np.cos(delta_lon))
    
    initial_bearing = np.arctan2(x, y)
    
    # Convert from radians to degrees
    initial_bearing = np.degrees(initial_bearing)
    
    # Normalize to 0-360
    bearing = (initial_bearing + 360) % 360
    
    return bearing

df_trip['bearing'] = calculate_bearing(
    df_trip['latitude'].shift(), 
    df_trip['longitude'].shift(), 
    df_trip['latitude'], 
    df_trip['longitude']
)

時系列プロット

Plotlyを使用して、データをプロットします。
共通のX軸にDatetimeを指定して、それぞれのカラムをY軸に設定します。

def make_timeseries_graph(df, x_col, y_cols):
    """
    共通のx軸で複数の折れ線をプロット
    """
    # サブプロットを作成し、x軸を共通にする
    fig = make_subplots(rows=len(y_cols), cols=1, shared_xaxes=True, 
                        subplot_titles=(y_cols),
                        vertical_spacing=0.1)

    for i, col in enumerate(y_cols):
        # プロット（サブプロット）
        fig.add_trace(go.Scatter(x=df[x_col], y=df[col], mode='lines+markers', name=col), row=i+1, col=1)

    # レイアウト設定
    fig.update_layout(
        height=200*len(y_cols),  # グラフの高さを調整
        title=f"Time Series of {', '.join(y_cols)}",
        showlegend=False,  # レジェンドをオフに
        plot_bgcolor='white',  # プロットの背景色を白に
    )

    # 各サブプロットのグリッド線の色をグレーに設定
    fig.update_xaxes(showgrid=True, gridcolor='lightgray')
    fig.update_yaxes(showgrid=True, gridcolor='lightgray')
    
    # グラフを表示
    fig.show()
    
make_timeseries_graph(df_trip, 'datetime',['speed', 'accel_x_transverse', 'accel_y_longitudinal', 'accel_z_vertical', 'gps_error_meter'])

Plotlyでグラフを作成すると、拡大縮小やプロットの値などをインタラクティブに確認することが可能です。

地図プロット

Foliumを使用して、緯度経度を地図にプロットします。
地図のタイルには国土地理院の衛星データを使用します。

10秒ごとに車両の進行方向をプロットしたいため、
カスタムのマーカーをHTMLで定義して、CSSで方向を指定します。

# カスタム三角形マーカーをHTMLとCSSで作成
def create_triangle_marker(points, rotation, m):
    triangle_html = f"""
    <div style="
        width: 0;
        height: 0;
        border-left: 5px solid transparent;
        border-right: 5px solid transparent;
        border-bottom: 10px solid orange;
        transform: rotate({rotation}deg);
        ">
    </div>
    """
    icon = folium.DivIcon(html=triangle_html)
    folium.Marker(location=points, icon=icon).add_to(m)

df_trip['points'] = df_trip.apply(lambda x: [x['latitude'], x['longitude']], axis=1)
m = folium.Map(tiles='https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg', attr='国土地理院', max_native_zoom=22)
folium.PolyLine(locations=df_trip['points'].to_list(), color="orange").add_to(m)
# 10レコードごとに向きを追加
for i,row in df_trip[::10].iterrows():
    create_triangle_marker(row['points'], row['bearing'], m)
m.fit_bounds(df_trip['points'].to_list())
m

スライダーを追加することで車両の軌跡をアニメーションで確認することも可能です。

m = folium.Map(tiles='https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg', attr='国土地理院')
m.fit_bounds(df_trip['points'].to_list())

data = []
for _, row in df_trip.iterrows():
    feature = {
        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [row["longitude"], row["latitude"]]
        },
        "properties": {
            "time": row["measurement_ms"],
            "popup": f"Speed: {row['speed']} m/s",
            "icon": "circle",
            "iconstyle": {
                "fillColor": "blue",
                "fillOpacity": 0.6,
                "stroke": "true",
                "radius": 5
            }
        }
    }
    data.append(feature)

TimestampedGeoJson(
    {
        "type": "FeatureCollection",
        "features": data,
    },
    transition_time=1,
    add_last_point=False,
    period="PT1S",
    auto_play=False,
    loop=False,
).add_to(m)
m

おわりに

今回は、1運行をFoliumとPlotlyを使用して可視化を行いました。

次回は、複数の運行からの統計分析やマップマッチングあたりの記事を作成したいと考えております。