駅圏の可視化とは？

ある地点からみた時に、どの駅が最寄り駅か知りたい時がある。
その駅が最寄りとなる領域を「駅圏」とし、可視化してみた。

領域の分割方法としては、ボロノイ分割を用いた。
ちょうど、空港での例を見つけたので、一度読んでみてほしい。

「ボロノイ図」というのは「ランダムに並んだ複数の点が存在している時、それぞれの点から距離の等しい位置を通る線(垂直二等分線)を引き、点ごとの領域を示した図」です。
文字で書くとわかりづらいですが、小学校の学区割りで、近隣の学校との距離を計算して極端に通学距離が長くならないようにエリア分けしたりするときに使われます。

山手線の座標

山手線停車駅の座標はこちらを利用した。

csvに変換して、 yamanote.csvというファイル名とした。

コード

scipyの中に Voronoiというクラスがある。これを利用することで、ボロノイ分割を簡単にできる。

参考: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Voronoi.html

from scipy.spatial import Voronoi, voronoi_plot_2d
import matplotlib.pyplot as plt
import pandas as pd

station = pd.read_csv("./yamanote.csv")

# 緯度経度のみ取得する
station_location = station[["longitude", "latitude"]].values

fig, ax = plt.subplots(figsize=(10, 10))

vor = Voronoi(station_location)

# 頂点は非表示
voronoi_plot_2d(vor, ax, show_vertices=False, point_size=5)

# 駅名を表示
station_x = station_location[:, [0]]
station_y = station_location[:, [1]]

駅名を表示させてないため、わかりづらい。

ここに駅名を載せる.

station_location = station[["longitude", "latitude"]].values

fig, ax = plt.subplots(figsize=(10, 10))

vor = Voronoi(station_location)

# vertexは非表示
voronoi_plot_2d(vor, ax, show_vertices=False, point_size=5)

# 駅名を取得する
station_names = station[["station"]].values.flatten() 

# 駅名を表示
station_x = station_location[:, [0]]
station_y = station_location[:, [1]]
for i, name in enumerate(station_names):
    ax.annotate(name, (station_x[i], station_y[i]))

適当にプロットして最寄り駅がどこか確認してみる

ここに、緯度 35.68, 経度 139.74をプロットしてみる。

ここからは、有楽町駅が近いことがわかる。

station_location = station[["longitude", "latitude"]].values

x = 139.74
y = 35.68

fig, ax = plt.subplots(figsize=(10, 10))

vor = Voronoi(station_location)

plt.scatter(x, y)

# vertexは非表示
voronoi_plot_2d(vor, ax, show_vertices=False, point_size=5)

# 駅名を取得する
station_names = station[["station"]].values.flatten() 

# 駅名を表示
station_x = station_location[:, [0]]
station_y = station_location[:, [1]]
for i, name in enumerate(station_names):
    ax.annotate(name, (station_x[i], station_y[i]))

詳細にデータを見る

ボロノイ分割の情報は、 Voronoi インスタンスの各種属性よりみることができる。

頂点

頂点（vertex）とは、分割時の垂直二等分線どうしの交点であり、以下のようにアクセスできる。

vor.vertices

# array([[139.70137049,  35.8642353 ],
#       [139.73438581,  35.6945346 ],
#       [139.73911137,  35.70021051],
#       [139.73373498,  35.66774947],
#      [139.96865369,  35.59334462],
#      [139.63225487,  35.67244304], ...

各領域

分割された各領域（ボロノイ要素）がどの頂点で構成されているかは、以下のように取得できる。配列の要素は、verticesのindexを指す。
-1を含むものの場合は、ボロノイ要素が閉じられてないことを示す。

vor.regions
# [[-1, 0, 9],
#  [16, 14, 13, 15],
#  [16, 11, 12, 4, 2, 15],
#  [20, 18, 17, 19],
#  [18, 1, 6, 7, 17],
#  [19, 3, 8, 7, 17],
#  [25, 21, 23, 22, 24],
#  [25, 1, 18, 20, 4, 2, 21],
#  [23, 9, 0, 13, 15, 2, 21], ....

vor.regions で出力される配列は、 dataframeのstationの順番になっている。これを利用し、各駅とボロノイ分割によりできた要素の面積を算出してみる。

from scipy.spatial import ConvexHull

# 座標に基づき、面積を算出
def voronoi_volumes(points):
    v = Voronoi(points)
    vol = np.zeros(v.npoints)
    for i, reg_num in enumerate(v.point_region):
        indices = v.regions[reg_num]
        if -1 in indices:
            vol[i] = np.inf
        else:
            vol[i] = ConvexHull(v.vertices[indices]).volume
    return vol

station["station_area"] = voronoi_volumes(station_location)

駅名	緯度	経度	緯度経度を数値としたときの面積
東京	35.681382	139.766084	0.001638
有楽町	35.675069	139.763328	0.001080
新橋	35.665498	139.759640	0.000976
浜松町	35.655646	139.756749	inf
田町	35.645736	139.747575	0.000785

infは、ボロノイ要素が頂点で閉じられてないために、面積の算出ができなかったもの。

参考