▲可視化の練習をしよう～Basemap編～

こんにちは。
テービーテックの村松です。

SIGNATEではCompetitionsにて練習問題が公開されています。
本ブログでも深津パイセンがチャレンジしておりましたね。
ds-blog.tbtech.co.jp

今回は練習問題の紹介ではなく、コンペで公開されているチュートリアルをやってみて覚えた可視化を覚書ついでにご紹介いたします。

因みに私が可視化の練習がてらチャレンジしたチュートリアルはこちら↓↓
signate.jp

コンペ自体はもう終了していますが、データはダウンロードできるのでぜひ皆さん遊びがてらチャレンジしてみてください。

「Basemap」って？

さて、今回は上記のチュートリアル内で使用されている「Basemap」についてご紹介していきます。
「Basemap」はざっくり説明するとPythonで地図を表示するライブラリです。
経緯度情報を地図上に落とし込むことが可能です。
なお、可視化で普段よくお世話になっているmatplotlibの一部となります。

「Basemap」インストール

「Basemap」のインストール方法は環境により色々な方法があるようです。
今回はGoogleColabにて。

!apt-get install libgeos-3.5.0
!apt-get install libgeos-dev
!pip install https://github.com/matplotlib/basemap/archive/master.zip

地図を表示してみよう

範囲を指定しない状態だと世界地図が表示されますが、
今回は日本辺りを指定しておきます。

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();

※北と南がちょん切れていて申し訳ない。日本が縦に長いから全部入れようと思うととても縦長になってしまうので今回は勘弁してください・・・。
カクカクな日本地図ですね。
もうちょっと詳細に滑らかな線で表示して欲しいところ。

調べたところBasemap()内にオプションで解像度を指定する「resolution=」がありました。
デフォルトだと「resolution='c'」で、粗い解像度になるようですね。
’c’以外に指定できる解像度は4種類。

'l'→low、低解像度
'i’→intermediate、中間解像度
'h'→high、高解像度
'f'→full、最高解像度

最高解像度はとても綺麗にですが、ちょっと重かったです。
ちょろっと表示させようかなってときは低解像度が早いし、見栄えも十分な印象でした。
表示範囲や用途によって変えていきたいですね。

それを踏まえてもうちょっと綺麗に表示させてみましょう。

# 地図の表示
north = 46.
south = 30.
east = 147.
west = 128.
m = Basemap(projection='merc', 
            llcrnrlat=south,
            urcrnrlat=north, 
            llcrnrlon=west,
            urcrnrlon=east, 
            resolution='l')

# 陸地を緑に
m.fillcontinents(color='g')
# 海を青に
m.drawlsmask(ocean_color='lightcyan')

# 海岸線を引く
m.drawcoastlines();

先ほどよりもそれっぽくなりましたね。
陸地と海もそれっぽい色に変えてみました。
ついでにやや地図の形が変わったことにお気づきでしょうか？
Basemap()内のオプションで「projection='merc'」を追加しました。
地図の表示の仕方をメルカトル図法にしたよって意味です。
この見え方がおそらく私たちが一番よく目にする日本地図の形でしょうかね。
他にも色々選べますので地図を正しく使わないと怒られる方は用途に合わせた投影図法を選んでくださいね。

ぶっちゃけ本ブログの趣旨に関係ないので見なくても良いのですが、
メルカトル図法ってなんぞ？って思ったら↓↓をどうぞ。
メルカトル図法の地図に騙されてる？本当の国のサイズを教えてくれる The Truth Size - PhotoshopVIP

経緯度情報を地図上に表示してみよう

まずは簡単に一点を表示してみましょう。
今回の経緯度は豊田市駅のものを使用しています。

豊田市駅

# 地図の表示(愛知県)
north = 35.5
south = 34.4
east = 138
west = 136.5
fig = plt.figure(figsize=(8, 9))
m = Basemap(projection='merc', 
            llcrnrlat=south,
            urcrnrlat=north, 
            llcrnrlon=west,
            urcrnrlon=east, 
            resolution='h')

#豊田市駅の経緯度を指定・表示
toyota_lon = 137.09228
toyota_lat = 35.0514
x,y = m(toyota_lon, toyota_lat)
m.plot(x, y,color='r', marker='o', markersize=10)
plt.text(x + 1500, y + 1500, u'TOYOTASHI')

m.drawcoastlines();

市の境界線がないのでわかりずらいですが、おそらく間違ってはなさそう。
「Basemap」では日本の都道府県や市町村の境界を表示する機能は標準装備されていないようです。
南北アメリカ、オーストラリアのみ「m.drawstates()」で州の境界線を描くことができるそうです。
参照：気象データ解析のためのmatplotlibの使い⽅

地価

続きまして国土数値情報ダウンロードサービスからcsvファイルを取得してもっとたくさんの点をプロットしてみましょう。
今回は「都道府県地価調査」を取得しました。

まずは下準備。

# 地価の情報をファイルから取得
import pandas as pd
tika = pd.read_csv('ファイル名.csv')

#平米単価から平米単価スコアを計算
tika['平米単価'] = tika['Ｒ１価格'] / tika['地積']
tika['平米単価スコア'] = np.fmin(1.0, tika['平米単価'] / (tika['平米単価'].median() * 3))

# 秒で緯度経度が入っているので変換
tika['緯度'] = tika['緯度'] / 3600
tika['経度'] = tika['経度'] / 3600

あとはさっきとさほど変わりません。

# 地図の表示(愛知県)
north = 35.5
south = 34.4
east = 138
west = 136.5
fig = plt.figure(figsize=(8, 9))
m = Basemap(projection='merc', 
            llcrnrlat=south,
            urcrnrlat=north, 
            llcrnrlon=west,
            urcrnrlon=east, 
            resolution='h')
m.drawcoastlines()

# 地価の表示
lon, lat = m(tika['経度'].values, tika['緯度'].values)
plt.scatter(lon, lat, s=15, c=tika['平米単価スコア'], cmap=plt.get_cmap('jet') , alpha=0.4);

赤が平米単価の高い場所、
青が平米単価の低い場所となっています。

・・・
名古屋やっべえ、ということがよくわかりますね。

以上、「Basemap」の覚書でした。

可視化は色々なツールがあって面白いですね。
また覚えた可視化があったら覚書ついでにご紹介いたします。