Pythonを用いた地理空間情報の可視化方法は数多くありますが、今回はBasemap Matplotlib Toolkitを使って位置情報を可視化してみます。
例えば地理空間情報可視化に関して、過去に紹介されている記事はこちら
１. Basemapについて
Basemapは簡単に言うと、様々な地図投影法や地図タイル、あるいは海岸線や河川などを描きつつ、matplotlibのプロット機能を追加していくことができます。主に、地球科学者の間で使われています。
The matplotlib basemap toolkit is a library for plotting 2D data on maps in Python.
Basemap does not do any plotting on it’s own, but provides the facilities to transform coordinates to one of 25 different map projections (using the PROJ.4 C library). Matplotlib is then used to plot contours, images, vectors, lines or points in the transformed coordinates. Shoreline, river and political boundary datasets (from Generic Mapping Tools) are provided, along with methods for plotting them. The GEOS library is used internally to clip the coastline and polticial boundary features to the desired map projection region.
Basemap is geared toward the needs of earth scientists, particularly oceanographers and meteorologists.
(https://matplotlib.org/basemap/users/intro.htmlより引用)
１–１. Basemapの基本
まずはcondaを使用して、basemapをインストールします。
$ conda install -c anaconda basemap
$ conda install -c conda-forge basemap-data-hires
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
fig = plt.figure()
# 地図を描く範囲を指定する。
m = Basemap(llcrnrlat=30, urcrnrlat=50, llcrnrlon=125, urcrnrlon=150) # Basemapインスタンス作成
# 緯度・経度を10度毎に引く。二つ目のオプションではラベルを上下左右どこに付けるかを設定している。
m.drawparallels(np.arange(-90, 90, 10), labels=[ True,False, True, False]) # 緯度線
m.drawmeridians(np.arange(0, 360, 10), labels=[True, False, False, True]) # 経度線
m.drawcoastlines() # 国境線
plt.show()
１–２. 地図の詳細設定
地図の解像度や投影法なども柔軟に変更できます。（細かい設定はたくさんあります、、、）
class mpl_toolkits.basemap.Basemap(llcrnrlon=None, llcrnrlat=None, urcrnrlon=None, urcrnrlat=None, llcrnrx=None, llcrnry=None, urcrnrx=None, urcrnry=None, width=None, height=None, projection='cyl', resolution='c', area_thresh=None, rsphere=6370997.0, ellps=None, lat_ts=None, lat_1=None, lat_2=None, lat_0=None, lon_0=None, lon_1=None, lon_2=None, o_lon_p=None, o_lat_p=None, k_0=None, no_rot=False, suppress_ticks=True, satellite_height=35786000, boundinglat=None, fix_aspect=True, anchor='C', celestial=False, round=False, epsg=None, ax=None)
| projection\resolution | 低解像度 | 中解像度 | 高解像度 |
|---|---|---|---|
| 正距円筒図法 |  |
 |
 |
| メルカトル図法 |  |
 |
 |
| ランベルト図法 |  |
 |
 |
ちなみに何も設定しなかった場合のデフォルト値は
projection='cyl'（正距円筒図法）
resolution='c'（粗い解像度）
２. 位置情報可視化方法
２–１. データソース
今回は無料公開されている擬似人流データをサンプルとして使わせてもらおうと思います。
２–２. データ構造
2013年9月16日の関西圏擬似人流データを読み込みんで、中身を確認します。
import pandas as pd
df = pd.read_csv('./Kansai/2013_09_16.csv')
df
'''
1列目：ユーザー ID
2列目：性別推定値（1:男性、2:女性、0・3：不明、NA：未推定）
3列目：日付・時刻（5分毎の24時間分）
4列目：緯度
5列目：経度
6列目：滞在者カテゴリ(大分類) ※文字列型
7列目：滞在者カテゴリ(小分類) ※文字列型
8列目：状態(滞在or移動) ※文字列型
9列目：滞在者カテゴリID(6、7行目に対応)
'''
今回は移動しているユーザ（STAY_MOVE=='MOVE'）を対象にします。
またtimestampを時分に分割しつつ、ユーザ毎・時間毎にランクを付与して扱いやすい形に整形しておきます。
from dfply import *
df = df >> filter_by(X.STAY_MOVE=='MOVE') >> select(columns_to(X.lon, inclusive=True))
df = df >> separate(X.timestamp, ['col1','hour','col2','minute','col3'], sep=[10,13,14,16],convert=True) >> select(~X.col1, ~X.col2, ~X.col3)
df = df >> group_by(X.uid,X.hour) >> mutate(rk=row_number(X.minute))
df
２–３. Basemapの初期値設定と位置情報の可視化
ではこれから位置情報を可視化していこうと思いますが、先に今回ベースとして考えていく地図を設定します。
ArcGISの背景地図や県境・市町村境を描画できるメソッドもあるため、追加でそれらも適用させていこうと思います。
県境や市町村境については__こちら__からShapefileをダウンロードしてください。
ファイル構造は以下の通り。
-gadm
-gadm36_JPN_1.cpg
-gadm36_JPN_1.shp
-gadm36_JPN_2.dbf
-gadm36_JPN_2.shx
-gadm36_JPN_1.dbf
-gadm36_JPN_1.shx
-gadm36_JPN_2.prj
-gadm36_JPN_1.prj
-gadm36_JPN_2.cpg
-gadm36_JPN_2.shp
地図の初期値設定
def basemap():
fig = plt.figure(dpi=300)
m = Basemap(projection="cyl", resolution="i", llcrnrlat=33.5,urcrnrlat=36, llcrnrlon=134, urcrnrlon=137)
m.drawparallels(np.arange(-90, 90, 0.5), labels=[True, False, True, False],linewidth=0.0, fontsize=8)
m.drawmeridians(np.arange(0, 360, 0.5), labels=[True, False, False, True],linewidth=0.0, rotation=45, fontsize=8)
m.drawcoastlines(color='k')
m.readshapefile('./gadm/gadm36_JPN_1', 'prefectural_bound1', color='k', linewidth=.8) # 県境
m.readshapefile('./gadm/gadm36_JPN_2', 'prefectural_bound2', color='lightgray', linewidth=.5) # 市町村境
m.arcgisimage(service='World_Street_Map', verbose=True, xpixels=1000, dpi=300)
ここまでやれば、あとはmatplotlibの感覚でプロットできます。
scatter(x, y, *args, **kwargs)
https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#scatter
tmp1=df[(df['gender']==1) & (df['hour']==9) & (df['rk']==1)]
tmp2=df[(df['gender']==2) & (df['hour']==9) & (df['rk']==1)]
basemap()
plt.scatter(tmp1['lon'],tmp1['lat'],color='b',s=0.5) # 男性を青で
plt.scatter(tmp2['lon'],tmp2['lat'],color='r',s=0.5) # 女性を赤で
ここではユーザ毎に、9時台の最初のログだけをプロットしています。
hexbin(x, y, **kwargs)
https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#hexbin
tmp=df[(df['hour']==9) & (df['rk']==1)]
basemap()
plt.hexbin(tmp['lon'],tmp['lat'], gridsize=50, cmap='rainbow', mincnt=1, bins='log',linewidths=0.3, edgecolors='k')
gridsizeでヘキサグリッドの大きさも変えられます。
３. おまけ
頻度分布に降水分布を重ね合わせて可視化する
例えば人の動きを左右しそうな気象データを掛け合わせることも可能。
ちなみに2013年9月16日は台風が通過した日でした。
 |
 |
|---|---|
| 気象庁より | １kmメッシュレーダーエコー |
人の移動速度と移動方向を可視化する
始点と終点の位置情報から速度と方位角を割り出してベクトル表示も可能。(使い道があるかは置いといて)
以上になります！
最後までお読みいただきありがとうございました！