Python geopandas + Bokeh で地理情報をプロットしたい

数日前、pandas を利用して地理情報をプロットするという非常によいエントリが翻訳されていた。

上のエントリ、前処理が手間に見えるが pd.read_html や .str アクセサを使えばもっと簡単に書けると思う、、、がそれは本題でない。

pandas で地理情報を扱う場合、geopandas という拡張パッケージを利用すると便利なため、その使い方を書きたい。また、処理を Python で完結させるため、QGIS ではなく Bokeh でプロットしたい。

`geopandas` のインストール

pip で。

$ pip install geopandas geopy

このエントリでは依存パッケージである shapely、geopy の機能も利用する。shapely は自動的にインストールされるはずだが、geopy については上のように別途インストールが必要。

地理情報の読み込み

以降の操作は IPython Notebook で行う。まずは必要なパッケージをインポートする。

%matplotlib inline
import numpy as np
import pandas as pd
import geopandas as gpd

例として、東京都の地図を読み込んで利用したい。geopandas は、依存パッケージである fiona がサポートする地理情報ファイル形式を読み込むことができる。GitHub に地球地図日本のデータから作成した GeoJSON ファイルがあったため、こちらを利用させていただく。

dataofjapan/land · GitHub
- "出典: 地球地図日本(国土地理院技術資料D・1-No.576)" を以下記載のスクリプトで加工して利用。

ファイルの読み込みは、geopandas.read_file で行う。データは pandas.DataFrame を継承した geopandas.GeoDataFrame インスタンスとして読み込まれる。

df = gpd.read_file('tokyo.geojson')

type(df)
# geopandas.geodataframe.GeoDataFrame

isinstance(df, pd.DataFrame)
# True

df.head()
#   area_en area_ja    code                                           geometry  
# 0  Tokubu     都区部  131211  POLYGON ((139.821051 35.815077, 139.821684 35....   
# 1  Tokubu     都区部  131059  POLYGON ((139.760933 35.732206, 139.761002 35....   
# 2  Tokubu     都区部  131016  POLYGON ((139.770135 35.705352, 139.770172 35....   
# 3  Tokubu     都区部  131067  POLYGON ((139.809714 35.728135, 139.809705 35....   
# 4  Tokubu     都区部  131091  (POLYGON ((139.719199 35.641847, 139.719346 35...

GeoDataFrame.plot を使うと、データを matplotlib を利用してプロットできる。

df.plot()

f:id:sinhrks:20150718210549p:plain

島嶼部があると細部がわからないため、都区部だけをフィルタしてプロットする。

df[df['area_en'] == 'Tokubu'].plot()

f:id:sinhrks:20150718210600p:plain

地理情報の処理

GeoDataFrame の各行はそれぞれひとつの領域 (ジオメトリオブジェクト) に対応している。各領域のジオメトリオブジェクトの実体は geometry キー中に shapely の Polygon もしくは MultiPolygon インスタンスとして保存されている。

Polygon には、領域各点の経度/緯度からなる tuple がリストとして保存されている。以下の処理で領域各点の経度/緯度をそれぞれ別のリストとして取り出すことができる。この処理は、続けて記載する Bokeh でのプロット時に必要になる。

# 1行目 (1領域目) を取得
s = df.iloc[0, :]
s
# area_en                                                Tokubu
# area_ja                                                   都区部
# code                                                   131211
# geometry    POLYGON ((139.821051 35.815077, 139.821684 35....
# ward_en                                             Adachi Ku
# ward_ja                                                   足立区
# Name: 0, dtype: object

type(s)
# <class 'pandas.core.series.Series'>

# 領域のジオメトリオブジェクトの型を表示
type(s['geometry'])
# <class 'shapely.geometry.polygon.Polygon'>

# 領域の経度/緯度を取得
list(s['geometry'].exterior.coords[:5])
# [(139.821051, 35.815077), (139.821684, 35.814887), (139.822599, 35.814509),
#  (139.8229070000001, 35.81437), (139.822975, 35.814339)]

# 経度、緯度別のリストに分割
lons, lats = zip(*list(s['geometry'].exterior.coords))

# 経度を表示
lons[:5]
# (139.821051, 139.821684, 139.822599, 139.8229070000001, 139.822975, 139.822977)

地理情報のプロット

このデータを Bokeh でプロットしたい。Bokeh についてはこちらを。

Chainer で Deep Learning: Bokeh で Live Monitoring したい - StatsFragments

Bokeh で地理情報をプロットするサンプルとして、ギャラリーにコロプレス図を描く例がある。Bokeh では描画する各領域の x 座標 (経度)、y 座標 (緯度) をリストで指定する。各領域は複数座標からなるため、"各領域の x/y 座標のリスト" のリスト ( 例: [[境界1の点1のx座標, 境界1の点2のx座標, ...], [領域2の点1のx座標, ...], ...] ) を引数として用意する必要がある。

import shapely

# "各領域の x 座標 (経度) のリスト" のリスト
countries_xs = []

# "各領域の y 座標 (緯度) のリスト" のリスト
countries_ys = []

for i, row in df.iterrows():
    if row['area_en'] != 'Tokubu':
        # 都区部のみ描画
        continue

    polygons = row['geometry']
    if isinstance(polygons, shapely.geometry.multipolygon.MultiPolygon):
        # MultiPolygon に含まれる Polygon を取り出す
        polygons = [x for x in polygons]
    elif isinstance(polygons, shapely.geometry.polygon.Polygon):
        polygons = [polygons]
    else:
        raise ValueError
        
    for p in polygons: 
        # 境界の各点を経度、緯度に分割
        lons, lats = zip(*list(p.exterior.coords))
        countries_xs.append(lons)
        countries_ys.append(lats)

作成したデータを Bokeh に渡す。出力先は IPython Notebook にした。

from bokeh import plotting

# 出力先を指定
plotting.output_notebook()

# figure を作成
p = plotting.figure(toolbar_location="left", plot_width=900, plot_height=700)

# 各領域を描画 (塗り分けはしない)
p.patches(countries_xs, countries_ys, fill_color='white',
          line_color="black", line_width=0.5)
plotting.show(p)

f:id:sinhrks:20150718210612p:plain

Geocoding の利用

geopandas は geopy を利用して Geocoding (住所から緯度経度への変換) を行うための API geopandas.geocode.geocode を持っている。が、geopandas v0.1.1, geopy v1.10.0 では、これらのインターフェースに不整合がありそのままでは利用できない。詳細と回避策は以下 Stack Overflow を。

python - Geocoding error with geopandas and geopy - Stack Overflow

from geopandas.geocode import geocode

# geopandas.geocode の API をそのまま使うとエラー
geocode(['千代田区千代田1-1'])
# AttributeError: 'module' object has no attribute 'MapQuest'

# 回避策として geopy の以下プロパティを明示的に None にする
import geopy
geopy.geocoders.MapQuest = None
geocode(['千代田区千代田1-1'])
#                                              address                        geometry
# 0  1-1 Chiyoda, Chiyoda-ku, Tōkyō-to 100-0001, Japan  POINT (139.7539454 35.6838012)

Geocoding を利用したデータの追加

上で作成した地図に、道路交通センサスにある各地点の交通量をプロットしたい。もっと細かいデータが見られるとよいのだが、公開情報ではこれくらいしかなさそうだ。以下のリンクから東京都の箇所別基本表 kasyo13.xls をダウンロードし、利用させていただく。

平成22年度道路交通センサス
- 出典：「平成22年度　全国道路・街路交通情勢調査（道路交通センサス）一般交通量調査　集計表」（国土交通省）

pdf = pd.read_excel('kasyo13.xls', header=8, index_col=False)
pdf.head()
#           (km)  Unnamed: 1  Unnamed: 2   Unnamed: 3  Unnamed: 4   Unnamed: 5  
# 0  13110100010           1        1010       東名高速道路           1  14110100090   
# 1  13111010010           1        1101  中央自動車道富士吉田線           4  13200400170   
# 2  13111010020           1        1101  中央自動車道富士吉田線           6          NaN   
# 3  13111010030           1        1101  中央自動車道富士吉田線           6          NaN   
# 4  13111010040           1        1101  中央自動車道富士吉田線           6          NaN

# ヘッダがセル結合されており正しく読めないため、連番で上書き
pdf.columns = range(len(pdf.columns))

# 計測地点の住所が入力されている列を表示
pdf.iloc[:, 28].head()
# 0             東京～東名川崎
# 1    高井戸ＩＣ～世田谷区・東京都境間
# 2    高井戸ＩＣ～世田谷区・東京都境間
# 3     世田谷区・東京都境～調布ＩＣ間
# 4    高井戸ＩＣ～世田谷区・東京都境間
# Name: 28, dtype: object

# NaN, 正しい住所でないデータをフィルタ
pdf = pdf.dropna(subset=[28])
pdf = pdf[pdf[28].str.startswith('東京都')]

# 計測地点住所と 24時間交通量の列のみにフィルタ
pdf = pdf.iloc[:, [28, 36]]
pdf.columns = ['計測地点', '交通量合計']
pdf.head()
#             計測地点   交通量合計
# 22  東京都中央区八丁堀３丁目  107885
# 23  東京都中央区八丁堀１丁目   98482
# 24   東京都中央区新富１丁目   89089
# 25   東京都中央区銀座１丁目   94143
# 26   東京都中央区築地１丁目   93844

上のデータから、計測地点 カラムの住所をもとに、Geocoding で緯度経度を取得する。

codes = geocode(pdf['計測地点'])

type(codes)
# <class 'geopandas.geodataframe.GeoDataFrame'>

codes.head()
#                                               address                        geometry
# 22  3 Chome Hatchōbori, Chūō-ku, Tōkyō-to 104-0032...  POINT (139.7756944 35.6755445) 
# 23  1 Chome Hatchōbori, Chūō-ku, Tōkyō-to 104-0032...  POINT (139.7771192 35.6774165)  
# 24  1 Chome Shintomi, Chūō-ku, Tōkyō-to 104-0041, ...  POINT (139.7725002 35.6729586) 
# 25  Ginza Itchome Station, 1 Chome-6 Ginza, Chūō-k...     POINT (139.767045 35.67435) 
# 26  1 Chome Tsukiji, Chūō-ku, Tōkyō-to 104-0045, J...  POINT (139.7716809 35.6702089)

取得したデータは上記のような GeoDataFrame となる。geometry カラムに経度/緯度を示す Point インスタンスが含まれる。この座標と交通量を利用してバブルチャートが描きたい。Bokeh のギャラリーにバブルチャートをプロットする例があるため、それに倣うと、

pdf = pdf.join(codes)

# 各点の x 座標
points_x = []

# 各点の y 座標
points_y = []

# 各点の円の半径
points_v = []

for i, row in pdf.iterrows():
    points_x.append(row['geometry'].x)
    points_y.append(row['geometry'].y)
    points_v.append(row['交通量合計'])
# 円の半径を適当に調整
points_v = np.sqrt(np.array(points_v)) / 50000

# Bokeh での描画
plotting.output_notebook()
p = plotting.figure(toolbar_location="left", plot_width=900, plot_height=700)
p.patches(countries_xs, countries_ys, fill_color='white',
          line_color="black", line_width=0.5)
# バブルチャートを描画
p.scatter(points_x, points_y, radius=points_v, fill_alpha=0.1)
plotting.show(p)