basemapやcartopyで地図上でplot

概要

地図と絡むことが増えてきたのでGoogleMapのAPI使おうかと思ったのだけど、今ひとつ要件を満たせないことがあったので、mpl_toolkits.basemapを使って緯度経度指定で地図にあれこれplotする方法を確認する。

インストール（Ubuntu）

matplotlib/basementのgithubのページからソースを落としてきてインストールする。

$ wget https://github.com/matplotlib/basemap/archive/v1.0.7rel.tar.gz
$ tar xvf v1.0.7rel.tar.gz
$ cd basemap-1.0.7rel/geos-3.3.3
$ configure
$ make
$ sudo make install

1つ上の階層にsetup.pyが置かれてるのでpipで叩く。

$ cd ../
$ sudo pip install .

ldconfigで /usr/local/lib/libgeos-3.3.3.so がいることを確認しておく。

$ sudo ldconfig
$ ldconfig -p | grep /usr/local/lib

これでインストール完了。下記のようにインストールを実行してエラーにならなければ成功。

$ python -c "from mpl_toolkits.basemap import Basemap"

mpl_toolkits.basemapで日本地図を描画

地図用のライブラリが入ったので、これで日本地図を描画してみる。

とりあえず何も指定せずに海岸線だけ引いてみる。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# 緯度経度で範囲を指定する
north = 46.
south = 30.
east = 147.
west = 128.

# 地図の表示
m = Basemap(llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east)
# 海岸線を引く
m.drawcoastlines() 

デフォルトのままだとちょっとガッカリな地図。

引数を指定すればもっと詳細な地図を出すこともできる。

resolutionの指定

コンストラクタを呼ぶ際にresolutionを指定すると、もう少し良い感じの地図が描ける。

具体的には、c (crude), l (low), i (intermediate), h (high), f (full) が指定できるらしい。地図を詳細にするとより綺麗な地図が描けるけど、それだけ描画には時間がかかるようになる。

試しにlowを指定してみる。

m = Basemap(llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east, resolution='l')
m.drawcoastlines() 

low, intermediate, high, fullを並べて表示してみる。

resolutions = ['l', 'i', 'h', 'f']
f, axarr = plt.subplots(2, 2)
for idx, resolution in enumerate(resolutions):
    ax = axarr[idx / 2][idx % 2]
    m = Basemap(llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east, resolution=resolution, ax=ax)
    m.drawcoastlines() 
    ax.set_title(resolution.upper())

デフォルトと比べるとlowでもそれなりに綺麗に見える。highやfullは描画に時間がかかるので、普段はもっぱらlowを利用する。

陸地と海を色分けする

白地図のままだと寂しいので色を付けてみる。

m = Basemap(projection='merc', llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east, resolution='h')

# 陸地を茶色に, 湖を水色に
m.fillcontinents(color='#8B4513', lake_color='#90FEFF')

# 海を濃い青に
m.drawlsmask(ocean_color='#00008b')

関東だけ拡大表示してみる

north = 36.5
south = 35.
east = 141.
west = 139.

# 地図の表示
m = Basemap(projection='merc', llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east, resolution='f')
m.drawcoastlines()
m.drawstates()

なお、県境等を表示する機能はデフォルトではない模様。

参考: Not defined inter-state lines in Japan? Any other countries?

上記によるとデータ落としてきてcartopy使えばいけるということだったのでやってみる。

とりあえず下記国土交通省の出してるデータから「行政区域」のデータを落としてくる。

国土数値情報ダウンロードサービス

落としてきたら、コードを実行するディレクトリにファイルを解凍しておく。

次にcartopyのインストール。condaを入れてれば簡単に入るけど（conda install -c scitools cartopy）、うちは入れてないのでソースから。

事前にproj.4 4.9.0がいるらしいので事前に入れておく。

$ git clone https://github.com/OSGeo/proj.4.git
$ cd proj.4
$ git checkout -b 4.9.0 tags/4.9.0
$ ./configure
$ make
$ sudo make install

libgeosも入れておく。

$ sudo apt-get install libgeos-dev

cloneしてタグ確認。

$ git clone https://github.com/SciTools/cartopy.git
$ cd cartopy
$ git tag

最新バージョンをチェックアウトしてpipでインストール。あとimportできるかチェック。

$ git checkout -b v0.14.2 tags/v0.14.2
$ sudo pip install .
$ python -c "import cartopy"

あとはStackoverflowを参考に国土交通省のファイルを読み込みつつ表示をしてみる。

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

# 落としてきた行政区域のshpファイルを指定
fname = 'N03-140401_GML/N03-14_140401.shp'
shapes = list(shpreader.Reader(fname).geometries())

# 東京あたりを描画
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_geometries(shapes, ccrs.PlateCarree(), edgecolor='black', facecolor='gray', alpha=0.3)
ax.set_extent([139, 141, 35, 36], ccrs.PlateCarree())
plt.show()

かなり細かく区分けされたデータが出た。

ax.stock_img()を使うと平地や海などが色分けされた状態で出力されたりする。

fname = 'N03-140401_GML/N03-14_140401.shp'
shapes = list(shpreader.Reader(fname).geometries())

ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_geometries(shapes, ccrs.PlateCarree(), edgecolor='black', facecolor='gray', alpha=0.3)
ax.set_extent([137, 141, 35, 39], ccrs.PlateCarree())
ax.stock_img()
plt.show()

駅の位置にラベルを付けてみる(basemap)

以前の記事でも使った駅データ.jpの情報から各駅の緯度経度を取得する。

試しに山手線の駅名を描画してみる。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# 日本語フォントの設定（takaoフォントが入ってる場合）
font = {'family' : 'TakaoGothic'}
matplotlib.rc('font', **font)

# 地図の表示
north = 35.8
south = 35.5
east = 140
west = 139.5
m = Basemap(projection='merc', llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east, resolution='h')
m.drawcoastlines()
m.drawstates()

# 駅情報取得
stations = pd.read_csv('station20150414free.csv')
lines = pd.read_csv('line20150414free.csv')
stations = stations[['station_cd', 'station_name', 'line_cd', 'lon', 'lat']].merge(lines[['line_cd', 'line_name']])

# 山手線の駅名表示
for idx, station in stations[stations.line_name == 'JR山手線'].iterrows():
        x, y = m(station.lon, station.lat)
        plt.text(x, y, station.station_name)

ちょっと上の方が被ってしまったけど、とりあえず地図上に表示できた。

駅の位置にラベルを付けてみる(cartopy)

Cartopyの方で大江戸線を出してみる。

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

# 落としてきた行政区域のshpファイルを指定
fname = 'N03-140401_GML/N03-14_140401.shp'
shapes = list(shpreader.Reader(fname).geometries())

# 東京あたりを描画
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_geometries(shapes, ccrs.PlateCarree(), edgecolor='black', facecolor='gray', alpha=0.3)
ax.set_extent([139.6, 139.9, 35.6, 35.8], ccrs.PlateCarree())
plt.show()

for idx, station in stations[stations.line_name == '都営大江戸線'].iterrows():
    plt.text(station.lon, station.lat, station.station_name)

東京の地価を表示してみる(basemap)

国土数値情報のページから地価のデータをcsvで落としてきて、UTF-8に変換したファイルをtika.csvという名前で保存しておく。

これを利用して平米単価が安いところを青で、高いところを赤でscatterでplotしてみる。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# 地図の表示
north = 35.85
south = 35.53
east = 140
west = 139
m = Basemap(projection='merc', llcrnrlat=south,urcrnrlat=north, llcrnrlon=west,urcrnrlon=east, resolution='f')
m.drawcoastlines()
m.drawstates()

# 地価の情報をファイルから取得
tika = pd.read_csv('tika.csv')
tika['平米単価'] = tika['Ｈ２７価格'] / tika['地積']

# 秒で緯度経度が入っているので変換
tika['緯度'] = tika['緯度'] / 3600
tika['経度'] = tika['経度'] / 3600
tika['平米単価スコア'] = np.fmin(1.0, tika['平米単価'] / (tika['平米単価'].median() * 3))

# 座標取得
lon, lat = m(tika['経度'].values, tika['緯度'].values)

plt.scatter(lon, lat, c=tika['平米単価スコア'], cmap=plt.get_cmap('jet') , alpha=0.7)

山手線の内側が高いのは当然として、中央線などもけっこう高めになっているのが目で見てわかる結果が表示される。

東京の地価を表示してみる(cartopy)

続いてcartopyでも同じものを出してみる。

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

# 落としてきた行政区域のshpファイルを指定
fname = 'N03-140401_GML/N03-14_140401.shp'
shapes = list(shpreader.Reader(fname).geometries())

# 東京あたりを描画
north = 35.85
south = 35.53
east = 140
west = 139
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_geometries(shapes, ccrs.PlateCarree(), edgecolor='black', facecolor='gray', alpha=0.3)
ax.set_extent([west, east, south, north], ccrs.PlateCarree())
plt.show()

# 地価の情報をファイルから取得
tika = pd.read_csv('tika.csv')
tika['平米単価'] = tika['Ｈ２７価格'] / tika['地積']

# 秒で緯度経度が入っているので変換
tika['緯度'] = tika['緯度'] / 3600
tika['経度'] = tika['経度'] / 3600
tika['平米単価スコア'] = np.fmin(1.0, tika['平米単価'] / (tika['平米単価'].median() * 3))

# 座標取得
plt.scatter(tika['経度'], tika['緯度'], c=tika['平米単価スコア'], cmap=plt.get_cmap('jet') , alpha=0.7)

ポイントが多過ぎてどこに区の境があるかわかりづらくなっている。区の名前を入れて地域を絞って見ればもう少しわかりやすい図になりそう。

感想

オープンなデータである程度やりくりできるということで、GoogleMapのAPIを使うよりは自由度があるものの、描画の時間もかかるし面倒さもいろいろあるなと。