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

はじめに

pythonを使って地図上にデータをプロットする。
地図はgoogle mapsではなく、ライセンス的に使いやすいOpenStreetMapを使用する。

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

© OpenStreetMap contributors

Requirements

• python3
• pandas
• matplotlib
• basemap

環境準備

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

basemapのインストール

basemapの公式サイトの手順に従って進める。まずは、依存ライブラリのインストール。

$ brew install geos
$ pip3 install pillow html5lib BeautifulSoup4

ここからbasemapソースコードをダウンロード。

basemapをインストール。

$ python3 setup.py install

データ取得

駅の緯度経度、乗降客数はそれぞれ、山手線停車駅の座標一覧と各駅の乗車人員 2014年度を利用させてもらう。

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

マップに描画

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')

結果

© OpenStreetMap contributors

参考

• matplotlib color example code