オープンデータから作るシェアサイクルネットワークデータ
GBFSについて
本記事はマイクロモビリティの標準的なデータフォーマットとして規格化されている「General Bikeshare Feed Specification(GBFS)」を用いてネットワークを可視化しました。
日本でも、ドコモ・バイクシェア、HELLO CYCLING、コトバイクなどで、GBFSでオープンデータ化されています。
GBFSに関しては以下の記事を御覧ください。
ネットワーク構造の概要
自転車シェアリングにおけるステーション間の関係性を把握するために、ステーションを「ノード（頂点）」として、各ステーションの1.5 km以内に位置する別のステーションとの間に「エッジ（辺）」を引くことで、ネットワーク構造を構築しました。このプログラムでは、1.5 km以内という条件で接続されるステーションのペアを生成し、ステーション間の距離や接続のパターンを分析できるようにしています。
このネットワークの全体構造を通して、どのステーションが他のステーションと密に接続されているか、つまり「中心的な役割」を果たしているのかを視覚的に捉えられます。例えば、利用頻度が高いエリアや交通のハブとなるステーションがネットワークの「中心」として可視化されやすくなります。こうした接続構造の理解は、ステーションの配置やメンテナンスの最適化、利用パターンの分析に役立ちます。
プログラム
1. データの読み込み:
CSVファイルからステーションデータを読み込みます。データには、ステーションの緯度経度やIDなどが含まれています。
2. Haversine関数の定義:
Haversineの公式を使って、2点間の直線距離を計算します。この距離は、緯度経度の座標を使って地球上の表面距離を求めるためのものです。
世界的に様々な場面で使用されている計算手法で，球面三角法を用いて地球を理想的な球体として扱います．ざっくり説明すると，球面上に北極点と扱う２点A，Bで球面三角形をつくり，球面三角法の余弦定理を用いて２地点間の距離Dを求めています．
他にも、2点間の直線距離の計算にはHubeny（ヒュベニ）の公式などもありますが、詳しくはこちらの記事を参考にしてみてください。
記事の中では以下の手法が解説されています。
- Haversineの公式
- Hubenyの公式
- Vincenty法
- 測地線航海算法(Geodesic Sailing)
- 国土地理院の方法
- GeographicLibライブラリを用いる
3. ステーションのペアを生成:
全ステーションのIDを取得し、すべてのステーションペアの組み合わせを作成します。これにより、各ステーションが他のすべてのステーションと組み合わせられ、隣接関係の確認が行いやすくなります。
4. 緯度経度のラジアン変換とKDTreeの構築:
距離計算を効率化するために、ステーションの緯度経度をラジアンに変換し、KDTree（高速な空間検索構造）を構築します。これにより、特定の地点から一定距離以内にある他の地点を素早く検索できます。
-
ラジアン変換について
緯度経度の角度情報を毎回ラジアンに変換しながら計算する場合、1つのペアの距離を求める際に複数の三角関数の計算が必要です。実際の時間短縮効果はデータ規模によりますが、数千レコードのデータに対して数十％程度の時間短縮になることがあります。 -
KDTreeについて
もっとも簡単な探索方法は，全ての点と比較して最も距離の近い点を探す方法です．この方法は全探索，或いは力任せ探索(brute force search)と呼ばれています。しかし，今回のように点群の数が多く，最近傍探索を実施する点も多い場合，計算コストは非常に大きくなります．kd-treeでは，２分木のデータ構造を利用した効果的な探索を行うことで，この計算コストを抑えることができます。
出典: Wikipedia - kd木
5. 1.5km以内のステーションペアの探索:
KDTreeを用いて、各ステーションから1.5 km以内にある他のステーションを検索します。指定した距離範囲内にあるステーションのインデックスを取得し、Haversine関数で実際の距離を確認します。
6. GeoJSON形式での保存:
1.5 km以内の距離にあるステーションペアをGeoJSON形式のラインストリングとして保存します。これにより、視覚化や地図上での解析に役立つ形式で出力されます。
import pandas as pd
import numpy as np
from shapely.geometry import LineString
from scipy.spatial import KDTree
import json
import itertools
import math
import geojson
from geojson import Feature, FeatureCollection, LineString
import time
# CSVファイルを読み込む
df = pd.read_csv('in/hello-cycling_20240620_0600.csv')
# レコード数を出力する
record_count = df.shape[0]
print(f'CSVファイルの結合が完了しました。{record_count}レコードです。')
# Haversine formulaを使って直線距離を計算する関数
def haversine(lat1, lon1, lat2, lon2):
R = 6371 # 地球の半径（キロメートル）
phi1, phi2 = math.radians(lat1), math.radians(lat2)
delta_phi = math.radians(lat2 - lat1)
delta_lambda = math.radians(lon2 - lon1)
a = math.sin(delta_phi / 2.0) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(delta_lambda / 2.0) ** 2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
return R * c
# ステーションのペアを生成する
stations = df['station_id'].tolist()
pairs = list(itertools.combinations(stations, 2))
num_pairs = len(pairs)
print(f"ステーションペアが {num_pairs} ペア作成されました。")
# プロファイリング用のタイマー開始
start_time = time.time()
print(f"処理開始時間: {start_time}")
# 緯度経度をラジアンに変換
df['lat_rad'] = np.radians(df['lat'])
df['lon_rad'] = np.radians(df['lon'])
# ラジアン変換の時間を表示
lat_lon_conversion_time = time.time()
print(f"緯度経度のラジアン変換時間: {lat_lon_conversion_time - start_time:.2f}秒")
# KDTreeを構築
print("KDTreeの構築を開始...")
coords = df[['lat_rad', 'lon_rad']].values
tree = KDTree(coords)
print("KDTreeの構築が完了しました。")
# 1.5km以内のステーションペアを探索
pairs_with_coords = []
features = []
for i, (lat1, lon1) in enumerate(zip(df['lat'], df['lon'])):
idx = tree.query_ball_point([df.at[i, 'lat_rad'], df.at[i, 'lon_rad']], 1.5 / 6371) # 1.5km以内のインデックスを取得
for j in idx:
if i < j: # 重複を避けるためにi < jのペアのみを保持
lat2, lon2 = df.at[j, 'lat'], df.at[j, 'lon']
distance = haversine(lat1, lon1, lat2, lon2)
if distance <= 1:
pairs_with_coords.append([df.at[i, 'station_id'], df.at[j, 'station_id'], lat1, lon1, lat2, lon2, distance])
# GeoJSON LineString featureを作成
line = LineString([(lon1, lat1), (lon2, lat2)])
feature = Feature(geometry=line, properties={
"distance": float(distance), # floatに変換
"station_id_1": str(df.at[i, 'station_id']), # strに変換
"station_id_2": str(df.at[j, 'station_id']) # strに変換
})
features.append(feature)
# 一定数のステーションごとに進捗を出力
if (i + 1) % 100 == 0 or i == len(df) - 1:
print(f"{i + 1}/{len(df)} ステーションを処理しました。")
# ペアをDataFrameに変換する
pairs_df = pd.DataFrame(pairs_with_coords, columns=['station_id_1', 'station_id_2', 'lat1', 'lon1', 'lat2', 'lon2', 'distance'])
# GeoJSONとして保存、ペアの数を報告
feature_collection = FeatureCollection(features)
file_name = f'out/station_pairs_within_1500m_{num_pairs}_pairs.geojson'
with open(file_name, 'w') as f:
geojson.dump(feature_collection, f)
print(f"ファイル '{file_name}' にステーションペアが保存されました。")
ネットワークの可視化
Mapbox Studioを使ってシェアサイクルのネットワークを可視化しました。
みなさんも、QGIS,kepler などお好みのツールを使って可視化してみてください。
まとめ
今回は、シェアサイクルのオープンデータを用いてネットワークを可視化しました。
やはり、ステーション（ポート）のピンを可視化するだけでは見えてこないネットワークの世界が見えてきてとても興味深いと改めて認識しました。
ぜひみなさんも手を動かして、GBFSデータを触ってみてください。
Discussion