PythonでGBFSを遊んでみる〜シェアサイクル到達圏の可視化〜

たまには本業っぽいことをしよう。

みなさんはシェアサイクルをご存じでしょうか？
日本国内だとHELLO CYCLINGやドコモ・バイクシェアが有名ですね。単位時間で自転車をレンタルして、任意のポート（またはステーションともいいますが、この記事ではポートに統一します）に返却できるというサービスです。

GBFS(General Bikeshare Feed Specification)というシェアサイクル（といいつつキックボードなども扱えます）の情報を扱うフォーマットが策定されており、日本でも先に挙げた2社は公共交通オープンデータセンターからこの形式でデータが公開されています。

GBFSとはなんぞやについては、@kumatira さんの記事に詳しいのでそちらを参照ください。

https://qiita.com/kumatira/items/f9229c21d9f0db3b5ae3

さて、この記事ですが、ラストワンマイルだなんだと前置きが長い記事なので、スパッとソースだけ見てえという方は、前置きは飛ばしてください。大したことは言ってないです。

:::note info
本記事では、OpenStreet株式会社（ハローサイクリング） / 公共交通オープンデータ協議会のバイクシェア関連情報(GBFS形式)をCC-BY4.0に従って利用しています。
https://ckan.odpt.org/dataset/c_bikeshare_gbfs-openstreet
:::

# やること is シェアサイクルはラストワンマイルの解決策たるや？を可視化する

ラストワンマイルとはなんぞやをつらつら書こうと思ったのですが冗長になるので省きます。
要は、既存の公共交通網、特に鉄道、バスから離れた箇所の交通手段としてシェアサイクルは有効かどうか？を可視化しますということで。

下記を可視化します。

1. 既存の鉄道駅、バス停の近傍（歩いて5分程度、1分80m/minとして４００m以内）のポート群Pを抽出
1. ポート群Pから自転車で30分以内に到達可能なポート群Pnを抽出
   1. ここで30分以内としているのは、現在日本国内でサービスしているシェアサイクルは、３０分過ぎると料金が増加する体系が多いからです
1. ポート群Pnから徒歩５分圏内をシェアサイクルを用いた場合の到達圏とする

![sharecycle_area.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/343307/1bdaff13-fe00-81be-64dd-b87460c0a64e.png)
赤線を鉄道、青線をバス、水色円をバス停とする。
×円はポートとし、緑は駅・バス停の近傍ポート群P、紫はPから３０分内のポート群Pn、水色はそれ以外のポートとする
破線円はポート群Pから３０分圏内、半透明な水色円はポート群Pnの徒歩５分圏内とする

この半透明な水色の円がシェアサイクルにおけるラストワンマイルの到達圏となります。
この到達圏を地図上に重ね合わせれば、シェアサイクルで最低料金で行ける範囲が明示化されます。

上の図は模式図なのでスカスカに見えますが、シェアサイクルのポートはなかなか密ですし、半径400mの円は0.5平方kmあるので、実際にはもっと広く見えるはずです。

# 実際やること

さて、ここで一つごめんなさいです。
上をやりたかったのですが、めんどくさいのでだいぶ端折ります。

というのもバス停の座標を探すのが面倒なのです。
私は、小平あたりでシビックテックというグループに一応属しているため、国分寺から小平あたりでオープンデータを用いて上の分析をしたいのですが、現状、このエリアを走る立川バス、西武バス、小平のコミュニティバスはオープンデータ化されていません。

なので上で定義した到達圏の分析で、バスはやめて、鉄道だけします。
~~鉄道は小平とその近辺の駅なら10駅行くか行かないかなので、OpenStreetMapあたりで拾います。~~(20駅もあったよ。後悔した)

# Python GBFSクライアントを入れる

Pythonでは[GBFSのクライアントモジュール](https://pypi.org/project/gbfs-client/)があるので、それを使います。
GBFSは単なるJSONなので、デフォルトのjsonモジュールでも扱えますが、日本語でこのクライアントを紹介する記事がなかったので、使います。

```Shell
pip install gbfs-client
```

これでインストールOKです。
使い方は公式サイトにある通りで、難しくありません。

```Python
from gbfs.services import SystemDiscoveryService

ds = SystemDiscoveryService()
ds.get_system_by_id('hellocycling')
```
```Shell
{'Country Code': 'JP',
 'Name': 'HELLO CYCLING\u3000-どこでも借りられて好きな場所で返せる自転車シェア',
 'Location': 'Tokyo, JP',
 'System ID': 'hellocycling',
 'URL': 'https://www.hellocycling.jp',
 'Auto-Discovery URL': 'https://api-public.odpt.org/api/v4/gbfs/hellocycling/gbfs.json',
 'Validation Report': 'https://gbfs-validator.netlify.app/?url=https%3A%2F%2Fapi-public.odpt.org%2Fapi%2Fv4%2Fgbfs%2Fhellocycling%2Fgbfs.json'}
```

GBFSは、GBFSで公開している事業者に一意のsystem_idを設定することを求めています。
get_system_by_id()には、そのsystem_idを指定します。
HELLO CYCLINGであれば`hellocycling`、ドコモ・バイクシェア（の東京エリアのみ）であれば、`docomo-cycle-tokyo`になります。
（余談ですが、@kumatiraさんのGBFSの記事のコメントでHELLO CYCLINGの中の人とシェアサイクルの運用エリアをどう分けるか？の議論がされていましたが、このidの付け方からすると、ドコモ・バイクシェアとしては仙台など別エリアはsystem_idごと分けるつもりのように思われますね）

system_idは、Githubの[GBFSのリポジトリ](https://github.com/NABSA/gbfs/blob/master/systems.csv)で調べられます。
また`ds.systems`でsystem.csv登録の事業者の情報は全て辞書で取れるので、そこから探してもOKです。

```Python
for x in ds.systems:
    if 'HELLO' in x['Name']:
        print(x)
```
```Shell
{'Country Code': 'JP', 'Name': 'HELLO CYCLING\u3000-どこでも借りられて好きな場所で返せる自転車シェア', 'Location': 'Tokyo, JP', 'System ID': 'hellocycling', 'URL': 'https://www.hellocycling.jp', 'Auto-Discovery URL': 'https://api-public.odpt.org/api/v4/gbfs/hellocycling/gbfs.json', 'Validation Report': 'https://gbfs-validator.netlify.app/?url=https%3A%2F%2Fapi-public.odpt.org%2Fapi%2Fv4%2Fgbfs%2Fhellocycling%2Fgbfs.json'}
```

サポートしている情報の一覧をとるには、下記のようにします。

```Python
client = ds.instantiate_client('hellocycling', 'jp')
client.feed_names
```
```Shell
['system_information',
 'station_information',
 'vehicle_types',
 'station_status']
```

system_information以下、GBFSの4つの情報がサポートされているのがわかります。

なおinstantiate_client()で'jp'を指定しないとエラーになります。
これはinstantiate_client()は、言語指定をしないと'en'が初期値になり、HELLO CYCLINGは現在日本語のデータしかなく、英語の情報が取れないためです。
これはドコモ・バイクシェアも同様です。

このclientを使ってstation_informationなどの情報を取って行きます。

# 小平近辺のHELLO CYCLINGのポート情報をとる

HELLO CYCLING全体をとると数が多いのと、先述の通り、小平あたりに絞るので、小平市と接している下記の市に存在するポートに絞ります。

* 小平市
* 東久留米市
* 西東京市
* 小金井市
* 国分寺市
* 東大和市
* 東村山市

東がつく市区町村多いな！？
武蔵野市は小平市に接してはいるものの、接しているのがほんの少しのため割愛します。

```Python
target_cities = [
    '小平市',
    '東久留米市',
    '西東京市',
    '小金井市',
    '国分寺市',
    '東大和市',
    '東村山市']

def is_target_cities(station_addr):
    for city_name in target_cities:
        if city_name in station_addr:
            return True
    return False

stations = [station for station in client.request_feed('station_information').get('data').get('stations') if is_target_cities(station['address'])]
```

```Shell
[{'lat': 35.758123,
  'lon': 139.467186,
  'name': 'サイクルスポット東村山店',
  'address': '東京都東村山市本町2-21-1',
  'station_id': '1556',
  'rental_uris': {'ios': 'https://www.hellocycling.jp/app/port/detail/1556?referrer=odpt',
   'web': 'https://www.hellocycling.jp/app/port/detail/1556?referrer=odpt',
   'android': 'https://www.hellocycling.jp/app/port/detail/1556?referrer=odpt'},
  'parking_hoop': False,
  'parking_type': 'street_parking',
  'contact_phone': '+8105038218282',
  'is_charging_station': False,
  'vehicle_type_capacity': {'num_bikes_now': 3,
   'num_bikes_limit': 4,
   'num_bikes_parkable': 1,
   'num_bikes_rentalable': 2}},
〜以下略〜
```

取れました。
GBFSのstation_informationにはaddressがあり、HELLO CYCLINGはサポートしてくれています。
なかったらlat/lonから逆ジオコーディングする必要がありましたが、楽でいいですね。

# 駅の座標を用意する

対象となる駅は小平近辺ということで、私の独断と偏見で下記に絞ります。

東伏見
西武柳沢
田無
花小金井
小平
久米川
東村山
国分寺
一橋学園
青梅街道
萩山
八坂
武蔵大和
多摩湖
恋ヶ窪
鷹の台
小川
東小金井
武蔵小金井
西国分寺

めんどいのでオープンデータで探したら、@uedayouさんが鉄道駅をオープンデータ化してくれていたので、ここから座標を取ることにします。
https://qiita.com/uedayou/items/b5131b5ca930fe0bef69

```Python
ekis = [
    ('西武鉄道/西武新宿線/東伏見'),
    ('西武鉄道/西武新宿線/西武柳沢'),
    ('西武鉄道/西武新宿線/田無'),
    ('西武鉄道/西武新宿線/花小金井'),
    ('西武鉄道/西武新宿線/小平'),
    ('西武鉄道/西武新宿線/久米川'),
    ('西武鉄道/西武新宿線/東村山'),
    ('西武鉄道/西武多摩湖線/国分寺'),
    ('西武鉄道/西武多摩湖線/一橋学園'),
    ('西武鉄道/西武多摩湖線/青梅街道'),
    ('西武鉄道/西武多摩湖線/萩山'),
    ('西武鉄道/西武多摩湖線/八坂'),
    ('西武鉄道/西武多摩湖線/武蔵大和'),
    ('西武鉄道/西武多摩湖線/西武遊園地'), #多摩湖
    ('西武鉄道/西武国分寺線/恋ヶ窪'),
    ('西武鉄道/西武国分寺線/鷹の台'),
    ('西武鉄道/西武国分寺線/小川'),
    ('東日本旅客鉄道/中央本線/東小金井'),
    ('東日本旅客鉄道/中央本線/武蔵小金井'),
    ('東日本旅客鉄道/中央本線/西国分寺'),
]

import requests

eki_pos = dict()
for eki in ekis:
    data = requests.get('https://uedayou.net/jrslod/{}.json'.format(eki)).json()
    eki_name = eki.split('/')[2]
    if eki_name == '西武遊園地':
        eki_name = '多摩湖'
    eki_pos.setdefault(eki_name, dict())
    eki_pos[eki_name]['lat'] = data['https://uedayou.net/jrslod/{}'.format(eki)]['http://www.w3.org/2003/01/geo/wgs84_pos#lat'][0]['value']
    eki_pos[eki_name]['lon'] = data['https://uedayou.net/jrslod/{}'.format(eki)]['http://www.w3.org/2003/01/geo/wgs84_pos#long'][0]['value']
eki_pos
```

2021年3月の改称が未反映らしく、多摩湖が西武遊園地で登録されているので、そこだけ名前を変えます。

```Shell
{'東伏見': {'lat': '35.72868', 'lon': '139.56435'},
 '西武柳沢': {'lat': '35.72863', 'lon': '139.55294'},
 '田無': {'lat': '35.72729', 'lon': '139.53987'},
 '花小金井': {'lat': '35.72641', 'lon': '139.5127'},
 '小平': {'lat': '35.73648', 'lon': '139.48944'},
 '久米川': {'lat': '35.75003', 'lon': '139.47209'},
 '東村山': {'lat': '35.7609', 'lon': '139.46591'},
 '国分寺': {'lat': '35.70105', 'lon': '139.47928'},
 '一橋学園': {'lat': '35.7217', 'lon': '139.4801'},
 '青梅街道': {'lat': '35.73134', 'lon': '139.4766'},
 '萩山': {'lat': '35.74094', 'lon': '139.47668'},
 '八坂': {'lat': '35.74459', 'lon': '139.46844'},
 '武蔵大和': {'lat': '35.75607', 'lon': '139.44427'},
 '多摩湖': {'lat': '35.76511', 'lon': '139.44288'},
 '恋ヶ窪': {'lat': '35.71174', 'lon': '139.46358'},
 '鷹の台': {'lat': '35.72344', 'lon': '139.46113'},
 '小川': {'lat': '35.73814', 'lon': '139.46359'},
 '東小金井': {'lat': '35.70152', 'lon': '139.52247'},
 '武蔵小金井': {'lat': '35.70097', 'lon': '139.50481'},
 '西国分寺': {'lat': '35.69972', 'lon': '139.46577'}}
```

# 徒歩5分の距離と自転車30分の距離をとる

この場合の距離は、単純なメートルではなく、楕円体をしている地球上において何度になるか？です。
今回、座標はdegree(度表記)なので、バッファを取るのに移動距離が何度になるかを指定する必要があります。

2点間の座標があって、その座標間の距離を取るのは割とよくある問題です。よくありすぎて、公式化されていますし、大体の言語で関数化もされていてどこかしらにソースもあります。

ところがある座標があって、そこからXメートル離れた場所の座標を求めるとなると途端に面倒臭くなります。
その辺りはこの記事に詳しいです。
https://qiita.com/kasugab3621/items/47a5077e31c5209e7f9c

今回欲しいのは、指定された座標（駅やシェアサイクルのステーション）からの半径Xメートルになるバッファです。
そのバッファに指定する距離を知りたいのです。
例えば徒歩5分は、80m/minとすると400m移動できます。駅から半径400mになるバッファを取るのに、指定する距離（座標平面上の距離）が欲しいわけです。

厳密に計算するとかなりややこしいので、端折ります。

今回、扱う座標はほぼ全て北緯35度に属します。北緯35度では、経度1度あたり91287.7885mです。
1mあたりの度を導いて、それに距離をかければ、大体の座標平面上の距離がでます。

```Python
#メートルを距離に変換(端折って経度1度長とする)
def get_distance(distance):
    x = distance / 91287.7885
    return x


#徒歩の距離
TOHO_DISTANCE = get_distance(5 * 80) #徒歩5分 * 80m/min

#自転車の距離
BIKE_DISTANCE = get_distance(30 * 200) #自転車30分 * 12km/h(= 200m/min)
```

念の為述べておきますが、今回みたいな遊びには十分な精度なものの、この関数、かなり雑なので仕事など精度を求める場合はきちんと作られた方がいいです。

なお下記の距離になります(徒歩5分以内、自転車30分以内)。
```Shell
0.0043817470723370635 0.06572620608505594
```

# shapelyで解析する

バッファに指定する距離を導いたところで、解析していきます。

まず駅とステーションの一覧をShapelyでPolygonやPointにしていきます。

```Python
import itertools
from shapely.geometry import Point

#鉄道駅を徒歩5分圏内のPolygonにする
eki_points = dict()
for eki in eki_pos:
    #駅から徒歩5分以内のバッファを作る
    eki_points[eki] = Point(float(eki_pos[eki]['lon']), float(eki_pos[eki]['lat'])).buffer(TOHO_DISTANCE)
    eki_points[eki].name = eki

#シェアサイクルステーションをPointにする
station_points = dict()
for station in stations:
    station_points[station['name']] = Point(station['lon'], station['lat'])
    station_points[station['name']].name = station['name']
```

ここから解析です。

1. 駅から徒歩5分圏内にあるステーションを抽出(P)
1. Pから自転車30分圏内にあるステーションを抽出(Pn)
1. Pnの徒歩5分圏内を示す

```Python
#駅から5分圏内のステーションを抽出し、自転車30分以内のポリゴンにする
near_stations = dict() #P
for pair in itertools.product(eki_points, station_points):
    if eki_points[pair[0]].contains(station_points[pair[1]]):
            near_stations[pair[1]] = station_points[pair[1]].buffer(BIKE_DISTANCE)

#駅から徒歩5分以内のステーションから自転車30分以内で行けるステーションから徒歩5分以内のポリゴンをとる
service_area = dict() #Pn
for pair in itertools.product(near_stations, station_points):
    if near_stations[pair[0]].contains(station_points[pair[1]]):
            service_area[pair[1]] = station_points[pair[1]].buffer(TOHO_DISTANCE)
```

やってみたら、抽出されたステーションの数は、GBFSからとった小平近辺のステーションの数と一致しました。
つまり小平近辺のステーションは、全て駅から自転車で30分以内で到達可能でした。

よくよく考えたら、自転車で12km/hで30分って、半径6kmあるので、これかなり広いやん…

# foliumでプロットする

解析結果をfoliumで可視化します。

foliumの使い方は、@Kumanuron-1910さんの記事と公式ドキュメントを参考にしました。
https://qiita.com/Kumanuron-1910/items/12ce7aa02922927de2f4

```Python
import folium

# 地図生成
folium_map = folium.Map(location = [35.73648, 139.48944], zoom_start=13) #小平駅の座標

#駅のマーカーを置く
for eki in eki_points:
    folium.Marker(
        location = [eki_pos[eki]['lat'], eki_pos[eki]['lon']],
        popup = eki
    ).add_to(folium_map)

#到達圏をプロットする
for area in service_area:
    folium.Circle(
        location=[station_points[area].y, station_points[area].x],
        popup = area,
        radius=400,
        color='#0000aa',
        fill_color='#0000ff'
    ).add_to(folium_map)

# 地図表示
folium_map
```
目印に駅のマーカーを置いて、到達圏は半径をメートル指定できるCircleを使います。
こちらができた図になります。

![スクリーンショット 2022-10-05 19.06.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/343307/abf1f0c8-45e2-9bcf-9596-703797883915.png)

なお駅近傍の自転車で30分以内の圏内は、半径6kmの円になるのでこのエリアがほぼ入ります。
実際は信号や道が曲がっていたり、人通りがあったりでもっと小さいと思います。
自転車ネットワーク上で30分以内のステーションを探すなど、そのあたりを加味すればもう少し精緻な到達圏になると思います。