ICカードデータを用いた鉄道利用ODの時系列クラスタリング

はじめに

e-Statで公開[link]されている第13回大都市交通センサスは，三大都市圏においてICカードを利用した鉄道での移動が記録されています．
入出場駅の情報，入場時間帯，所要時間等の情報が記録されており，多様な分析に活用することができるものです．
本稿では，下記の投稿に着想を得て，鉄道利用ODの時系列クラスタリングを行おうと思います．

データの準備

今回用いるのは，第13回大都市交通センサスの1次ODデータです．
他の公開データとの違いは下の記事をご覧ください．

今回は大阪駅着ODのクラスタリングを行おうと思います．
そこで，公開されている1次ODデータのcsvを結合し，大阪駅着ODを抽出します．
1日目と2日目でそれぞれ6個ずつのcsvが用意されているのですが，各日2個目以降のcsvは列名がセットされていないというゴミ仕様なので，その点気を付けて結合します．

data_Osaka_destination_preparation.py

# ライブラリのインポート
import pandas as pd

# 1日目のデータの結合
df_1 = []

for i in range(1, 7):
    filename = f"1ji_1_{i}.csv"
    if i == 1:
        df = pd.read_csv(filename, header=0)  # 最初の行を列名として読み込む
        columns = df.columns  # 列名を保存
    else:
        df = pd.read_csv(filename, header=None)
        df.columns = columns  # 列名を手動で設定
    df_1.append(df)

combined_df_1 = pd.concat(df_1, ignore_index=True)

# 2日目のデータの結合
df_2 = []

for i in range(1, 7):
    filename = f"1ji_2_{i}.csv"
    if i == 1:
        df = pd.read_csv(filename, header=0)  # 最初の行を列名として読み込む
        columns = df.columns  # 列名を保存
    else:
        df = pd.read_csv(filename, header=None)
        df.columns = columns  # 列名を手動で設定
    df_2.append(df)

combined_df_2 = pd.concat(df_2, ignore_index=True)

# 1日目と2日目のデータの結合
df = pd.concat([combined_df_1, combined_df_2], ignore_index=True)

# 大阪駅着ODの抽出
df_Osaka_origin = df[df['【出場】駅名'] == "大阪"].reset_index(drop=True)

クラスタリング

以下のサイトを参考にクラスタリングを行いました．

今回は，6時から22時までの各入場時間帯の利用者推移でクラスタリングを行います．

clustering.py

# ライブラリのインポート
import pandas as pd
import numpy as np
from tslearn.metrics import dtw 
from tslearn.clustering import TimeSeriesKMeans
import matplotlib.pyplot as plt
import japanize_matplotlib  # グラフに日本語を表示しないなら不要
import matplotlib.cm as cm
import math

df = pd.read_csv("data_Osaka_destination.csv", encoding='shift_jis')

df_OD = df[['【入場】駅名', '【出場】駅名', '【入場】時間帯', '人数']]
df_OD = df_OD.groupby(['【入場】駅名', '【出場】駅名', '【入場】時間帯'], as_index=False)['人数'].sum()
df_OD['OD_pair'] = df_OD['【入場】駅名'].str.cat(df_OD['【出場】駅名'], sep='->')
df_OD = df_OD.drop(['【入場】駅名', '【出場】駅名'], axis=1)
df_OD = df_OD[(df_OD['【入場】時間帯'] >= 6)&(df_OD['【入場】時間帯'] <= 22)]

pivot_df_OD = df_OD.pivot_table(index='OD_pair', columns='【入場】時間帯', values='人数', fill_value=np.nan)
pivot_df_OD = pivot_df_OD.dropna()

クラスタリングの準備として，標準化を行います．

clustering.py

# 行方向に標準化
pivot_df_OD_standardized = pivot_df_OD.sub(pivot_df_OD.mean(axis=1), axis=0)  # 行ごとに平均値を出して，差を求める
pivot_df_OD_standardized = pivot_df_OD_standardized.div(pivot_df_OD.std(axis=1), axis=0)  # 行ごとに標準偏差を出して，商を求める

# 標準化後のグラフ描画（オプション）
plt.figure(figsize=(12, 6))

for od_pair in pivot_df_OD_standardized.index:
    plt.plot(pivot_df_OD_standardized.columns, pivot_df_OD_standardized.loc[od_pair], label=od_pair)

plt.title("各ODペアの時間帯別利用者数（標準化済）")
plt.xlabel("入場時間帯")
plt.ylabel("標準化スコア")
plt.xticks(rotation=45)
plt.legend(loc='upper right', bbox_to_anchor=(1.2, 1.0), fontsize=8)
plt.tight_layout()
plt.show()

標準化後の推移を確認すると以下のようになります．
まぁぐちゃぐちゃすぎて意味わからんです．

というわけで，参考サイトに従って，距離指標を動的時間伸縮法(Dynamic Time Warping: DTW)としたk-meansクラスタリングを実行します．
まずはエルボー法でクラスタ数を決定します．

clustering.py

distortions = [] 

for i in range(1,11): 
    ts_km = TimeSeriesKMeans(n_clusters=i,metric="dtw",random_state=42) 
    ts_km.fit_predict(pivot_df_OD_standardized) 
    distortions.append(ts_km.inertia_) 

plt.plot(range(1,11),distortions,marker="o") 
plt.xticks(range(1,11)) 
plt.xlabel("Number of clusters") 
plt.ylabel("Distortion") 
plt.show()

もう少しはっきり折れ曲がってくれると嬉しいのですが，今回はクラスタ数を5に設定しました．

clustering.py

# エルボー法の結果からクラスタ数を5にする
n_clusters = 5

# 決定したクラスタ数で再度クラスタリング実行
ts_km = TimeSeriesKMeans(n_clusters=n_clusters,metric="dtw",random_state=42) 
y_pred = ts_km.fit_predict(pivot_df_OD_standardized)

# 可視化
n_clusters = n_clusters
cluster_labels = np.unique(y_pred)

cmap = cm.get_cmap('tab10', n_clusters)  # カラーマップの取得
cluster_colors = [cmap(i) for i in range(n_clusters)]

n_cols = 3  # 3列で表示
n_rows = math.ceil(n_clusters / n_cols)

fig, axes = plt.subplots(n_rows, n_cols, figsize=(5 * n_cols, 4 * n_rows), sharex=True, sharey=True)
axes = axes.flatten()

for cluster_label in cluster_labels:
    ax = axes[cluster_label]
    cluster_indices = np.where(y_pred == cluster_label)[0]
    cluster_color = cluster_colors[cluster_label]
    
    for idx in cluster_indices:
        od_pair = pivot_df_OD_standardized.index[idx]
        ax.plot(
            pivot_df_OD_standardized.columns,
            pivot_df_OD_standardized.loc[od_pair],
            color=cluster_color,
            alpha=0.6
        )
    
    ax.set_title(f"クラスタ {cluster_label}")
    ax.set_xlabel("入場時間帯")
    ax.set_ylabel("標準化スコア")
    ax.tick_params(axis='x', rotation=45)

for i in range(len(cluster_labels), len(axes)):
    fig.delaxes(axes[i])

plt.tight_layout()
plt.show()

乗降客数の推移をクラスタごとに可視化した結果が下記になります．

クラスタ1は朝ラッシュ時間帯における利用者が相対的に多いクラスタでした．
他方，クラスタ4は夕ラッシュ時間帯における利用者が相対的に多いクラスタでした．

各ODペアの発駅の分布を，QGISでクラスタごとに色分けしてプロットすると，下記のようになりました．¹²

クラスタ0や1は都市間の地域に多く，住宅街が多いために朝ラッシュ帯の利用者が相対的に多くなったと考えられます．
他方，大阪駅周辺や大阪駅から離れた都心の駅はクラスタ4に振り分けられました．
前者は退勤後の娯楽/立ち寄り需要，後者は遠距離通勤者の乗り換えや立ち寄りなどによる出場記録ではないかと考えられます．

おわりに

今回は鉄道ICカードデータを用いた時系列クラスタリングを行ってみました．
もう少しきれいに分かれてくれるとより解釈しやすい結果になったかなという不満もありますが，朝ラッシュ型と夕ラッシュ型の違いが見られただけでも面白いかなと思っています．
様々なODペアやクラスタリング手法で試行してみるとまた異なった結果が得られるかと思います．

1. 「国土数値情報（鉄道データ）」³（国土交通省） (https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_3.html) をもとに著者作成 ↩

2. 背景地図はOpen Street Map (https://www.openstreetmap.org/copyright) をもとに作成 ↩

3. 出典：国土交通省国土数値情報ダウンロードサイト (https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_3.html) ↩