駅圏の可視化



## 駅圏の可視化とは？

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

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

https://gigazine.net/news/20140516-world-voronoi-airport/

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

## 山手線の座標

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

https://qiita.com/butchi_y/items/3a6b70b38e13dc56ef13

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

## コード

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

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


```py
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]]
```

![no_station_name.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88357/cc0946a2-e2d9-8e3e-a936-29e5468ff1d7.png)

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

ここに駅名を載せる.


```py
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]))
```

![f1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88357/c954f140-c7aa-1b07-ee78-89647a2e533a.png)

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

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

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

![f4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88357/583571cd-38fd-df33-6cf7-f12b78fd482d.png)


```py
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）とは、分割時の垂直二等分線どうしの交点であり、以下のようにアクセスできる。

```py
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を含むものの場合は、ボロノイ要素が閉じられてないことを示す。

```py
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`の順番になっている。これを利用し、各駅とボロノイ分割によりできた要素の面積を算出してみる。


```py
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は、ボロノイ要素が頂点で閉じられてないために、面積の算出ができなかったもの。



## 参考

https://qiita.com/kwi/items/4b7162fbada390501d06