python + OpenStreetMapで地図にデータをプロットする

# はじめに
pythonを使って地図上にデータをプロットする。
地図はgoogle mapsではなく、ライセンス的に使いやすい[OpenStreetMap](http://www.openstreetmap.org/)を使用する。

サンプルデータとして、今回は山手線の乗降客数を用いる。
完成図は以下。

![ダウンロード.png](https://qiita-image-store.s3.amazonaws.com/0/25443/4cb70cfc-8503-1ca8-9dd3-fb3c5aca9c48.png)
[© OpenStreetMap contributors](http://www.openstreetmap.org/copyright)

## Requirements
- python3
- pandas
- matplotlib
- basemap


#  環境準備
環境はMac、また、python3, pandas, matplotlibのインストールは完了しているものとする。

## basemapのインストール
[basemapの公式サイト](http://matplotlib.org/basemap/users/installing.html)の手順に従って進める。まずは、依存ライブラリのインストール。

```shell-session
$ brew install geos
$ pip3 install pillow html5lib BeautifulSoup4
```

[ここ](https://sourceforge.net/projects/matplotlib/files/matplotlib-toolkits/)からbasemapソースコードをダウンロード。

basemapをインストール。

```
$ python3 setup.py install
```


# データ取得
駅の緯度経度、乗降客数はそれぞれ、[山手線停車駅の座標一覧](http://qiita.com/butchi_y/items/3a6b70b38e13dc56ef13)と[各駅の乗車人員 2014年度](http://www.jreast.co.jp/passenger/index.html)を利用させてもらう。

```py3
import sys
import pandas as pd
from http.client import IncompleteRead
import time

# Retrieve station's lat/lon
def fetch_latlon():
    url = 'http://qiita.com/butchi_y/items/3a6b70b38e13dc56ef13'
    data = pd.read_html(url)[0]
    data.columns = ['name', 'lat', 'lon']
    return data

# Retrieve the number of passenger
def fetch_passenger():
    # fetch data
    urls = ['http://www.jreast.co.jp/passenger/index.html', 'http://www.jreast.co.jp/passenger/2014_01.html']
    data = pd.DataFrame()
    for url in urls:
        df = None
        for i in range(5):
            try:
                df = pd.read_html(url)
                break
            except IncompleteRead as ir:
                time.sleep(0.1)
        # Sometimes fetching passenger fails
        if df is None:
            print('Failed to fetch passenger.')
            sys.exit(1)
        
        for i in range(2):
            # Delete rank column in the first table
            if url == urls[0]:
                df[i] = df[i].iloc[:, 1:]
                df[i].columns = [col for col in range(len(df[i].columns))]
            data = data.append(df[i])
    # arrange data
    data = data.iloc[2:, [0, 3]]
    data.columns = ['name', 'passenger']
    data.reset_index(inplace=True, drop=True)
    return data

# Retrieve lat/lon and passenger info and merge them into one data
def fetch_data():
    latlon = fetch_latlon()
    passenger = fetch_passenger()
    data = latlon.merge(passenger, how='inner', on='name')
    return data

data = fetch_data()
```

取得結果は以下。
`data.iloc[:3]`

|  | name | lat | lon | passenger |
|:---|:--------|------------:|------------:|----------:|
| 0 | 東京 | 35.681382 | 139.766084 | 418184 |
| 1 | 有楽町 | 35.675069 | 139.763328 | 165450 |
| 2 | 新橋 | 35.665498 | 139.759640 | 253874 |


# マップに描画
```py3
import requests
from io import BytesIO
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from PIL import Image
%matplotlib inline

# min/max circle size of plot
min_size = 200
max_size = 2000

# Retrive static OpenStreetMap
def get_osm_img(minlat, minlon, maxlat, maxlon, scale=60000, img_format='png'):
    url = 'http://www.openstreetmap.org/export/finish'
    payload = {
        'mapnik_format': img_format, 
        'mapnik_scale': scale, 
        'minlon': minlon, 
        'minlat': minlat, 
        'maxlon': maxlon, 
        'maxlat': maxlat, 
        'format': 'mapnik'
    }
    response = requests.post(url, payload)
    return Image.open(BytesIO(response.content))

fig = plt.figure(figsize=(15, 15))

minlat, minlon, maxlat, maxlon = 35.61, 139.67, 35.75, 139.80
bmap = Basemap(projection='merc', llcrnrlat=minlat, urcrnrlat=maxlat, llcrnrlon=minlon, urcrnrlon=maxlon, lat_ts=0, resolution='l')

x, y = bmap(data['lon'].values, data['lat'].values)

file_name = 'osm.png'
bg_img = None
try:
    bg_img = Image.open(file_name)
except FileNotFoundError as fnfe:
    bg_img = get_osm_img(minlat=minlat, minlon=minlon, maxlat=maxlat, maxlon=maxlon, scale=60000)
    bg_img.save(file_name)
bmap.imshow(bg_img, origin='upper')
bmap.scatter(x, y, c=data['passenger'], cmap=plt.cm.get_cmap('seismic'), alpha=0.5, s=data['passenger'].map(lambda x: (x - data['passenger'].min()) / (data['passenger'].max() - data['passenger'].min()) * (max_size - min_size) + min_size))

# plt.colorbar()
plt.show()
# plt.savefig('out.png')
```

# 結果
![ダウンロード.png](https://qiita-image-store.s3.amazonaws.com/0/25443/4cb70cfc-8503-1ca8-9dd3-fb3c5aca9c48.png)
[© OpenStreetMap contributors](http://www.openstreetmap.org/copyright)

# 参考
- [matplotlib color example code](http://matplotlib.org/examples/color/colormaps_reference.html)