はじめに
『Pythonで学ぶ空間データサイエンス入門』の独学ノートです。本の内容から寄り道してアレコレ考えます。
本を読んだ上で補助的に読んでください。
この記事では、空間重み行列について、図を使って解説します。
【前の内容】
【他の内容】
【今回の内容】
2.2 空間重み行列の可視化
空間重み行列(spatial weight matrix)を数式とプログラムや図で確認します。
空間隣接行列(spatial adjacency matrix)については「【Python】2.1：空間隣接行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。
利用するライブラリを読み込みます。
# ライブラリを読込 import geopandas as gpd from pysal.lib import weights import numpy as np import matplotlib.pyplot as plt import japanize_matplotlib import matplotlib.cm as cm from matplotlib.animation import FuncAnimation
空間重み行列の数式表現
まずは、空間重み行列の定義や計算を数式で確認します。
の(
から
のインデックスで表す)
個の区域に関する空間隣接行列を
、空間重み行列を
とします。
番目と
番目の区域が隣接していない場合は
です。
番目の区域自体(
番目同士の組み合わせの場合)は隣接していない
とします。
番目と
番目の区域が隣接している場合は
です。
の各要素
は、
の行
ごとに各要素を総和(行和)で割った値で定義されす。
は、0から1の値で行和が1になる値に基準化されます。
列和 は、基本的に1になりません(たまたま1になることはあります)。
ここまでは、空間重み行列の定義や性質、空間隣接行列との関係を数式で表現しました。
データの読込
可視化の前に、例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。
「国土数値情報ダウンロードサイト」の行政区域データと地価公示データを読み込みます。
処理コード(クリックで展開)
# ファイルパスを指定 DISTRICT_PATH = 'N03-20230101_27_GML/N03-23_27_230101.shp' LAND_PATH = 'L01-23_27_GML/L01-23_27.shp' # データを読込 gdf_district = gpd.read_file(DISTRICT_PATH, encoding='shift-jis') # 行政区域 gdf_land = gpd.read_file(LAND_PATH) # 地価公示 # 対象データを取得 gdf_district = gdf_district[['N03_003', 'N03_004', 'N03_007', 'geometry']] gdf_district.columns = ['city1', 'city2', 'd_code', 'geometry'] gdf_land = gdf_land[['L01_024', 'L01_022', 'L01_006', 'geometry']] gdf_land.columns = ['address', 'd_code', 'price', 'geometry'] # データを整形 dst_proj = 6668 gdf_district = gdf_district.to_crs(epsg=dst_proj) gdf_land = gdf_land.to_crs(epsg=dst_proj) # 空間座標系を再設定 gdf_district['city2_label'] = gdf_district['city2'].str.replace('大阪市|堺市', '', regex=True) # 地名が重なる対策用 # データを加工 gdf_land['log_price'] = np.log(gdf_land['price'])
公示価格の対数をとった列を作成します。その他にも扱いやすいように整形しておきます。
行政区域データと地価公示データをまとめます。
# データを統合 gdf_target = gdf_district.merge( gdf_land.groupby('d_code')['log_price'].mean(), on='d_code' ) gdf_target = gdf_target.dissolve(by=['d_code'], as_index=False) gdf_target['centers'] = gdf_target['geometry'].centroid # データフレームを整形 gdf_target = gdf_target.reindex(columns=['city1', 'city2', 'd_code', 'geometry', 'centers', 'log_price', 'city2_label']) gdf_target
| city1 | city2 | d_code | geometry | centers | log_price | city2_label | |
|---|---|---|---|---|---|---|---|
| 0 | 大阪市 | 大阪市都島区 | 27102 | POLYGON ((135.51481 34.72057, 135.51489 34.720... | POINT (135.52721 34.71193) | 12.769969 | 都島区 |
| 1 | 大阪市 | 大阪市福島区 | 27103 | POLYGON ((135.46247 34.6874, 135.46253 34.6879... | POINT (135.47534 34.6939) | 13.214822 | 福島区 |
| 2 | 大阪市 | 大阪市此花区 | 27104 | MULTIPOLYGON (((135.36506 34.6325, 135.35885 3... | POINT (135.42012 34.668) | 12.191524 | 此花区 |
| 3 | 大阪市 | 大阪市西区 | 27106 | POLYGON ((135.46667 34.67284, 135.4642 34.6747... | POINT (135.48375 34.67776) | 13.746399 | 西区 |
| 4 | 大阪市 | 大阪市港区 | 27107 | MULTIPOLYGON (((135.43009 34.64977, 135.4301 3... | POINT (135.45084 34.66032) | 12.410058 | 港区 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 67 | 泉南郡 | 田尻町 | 27362 | MULTIPOLYGON (((135.29259 34.40094, 135.29293 ... | POINT (135.25775 34.417) | 10.923984 | 田尻町 |
| 68 | 泉南郡 | 岬町 | 27366 | MULTIPOLYGON (((135.18163 34.34026, 135.18187 ... | POINT (135.15203 34.30321) | 9.918761 | 岬町 |
| 69 | 南河内郡 | 太子町 | 27381 | POLYGON ((135.66132 34.53314, 135.66139 34.532... | POINT (135.65269 34.51765) | 10.697891 | 太子町 |
| 70 | 南河内郡 | 河南町 | 27382 | MULTIPOLYGON (((135.64123 34.45932, 135.64105 ... | POINT (135.64897 34.48455) | 10.327751 | 河南町 |
| 71 | 南河内郡 | 千早赤阪村 | 27383 | POLYGON ((135.62578 34.47662, 135.62579 34.476... | POINT (135.64812 34.43844) | 10.122623 | 千早赤阪村 |
72 rows × 7 columns
市区町村(標準地行政区域コード)ごとに、対数公示価格の平均をとった列を作成して、ポリゴンデータと結合します。また、ポリゴンデータをまとめて、中心座標を求めます。
以上で、作図に利用する空間データを用意できました。
空間重み行列の作成
次は、空間データから空間重み行列を作成します。
空間隣接行列を作成します。
# 隣接関係を作成 wr = weights.Rook.from_dataframe(gdf_target) # 空間隣接行列を作成 adj_mat = wr.full()[0].astype(int) print(adj_mat) print(adj_mat.shape)
[[0 0 0 ... 0 0 0]
[0 0 1 ... 0 0 0]
[0 1 0 ... 0 0 0]
...
[0 0 0 ... 0 1 0]
[0 0 0 ... 1 0 1]
[0 0 0 ... 0 1 0]]
(72, 72)
weights ライブラリを利用して、空間データから隣接関係を作成して、空間隣接行列に変換します。
空間隣接行列 を
adj_mat とします。ただし、Pythonはインデックスを0から数えるので、各要素は と
adj_mat[i-1, j-1] が対応します。
空間隣接行列の列和と行和の処理を確認します。
# 列和を計算 sum_col = adj_mat.sum(axis=0, keepdims=True) print(sum_col[:, :5]) print(sum_col.shape) # 行和を計算 sum_row = adj_mat.sum(axis=1, keepdims=True) print(sum_row[:5, :]) print(sum_row.shape)
[[5 5 4 7 3]]
(1, 72)
[[5]
[5]
[4]
[7]
[3]]
(72, 1)
np.sum() の axis 引数に 0 を指定すると0軸方向の和 、
1 を指定すると1軸方向の和 を計算します。基準化の計算のために、
keepdims=True を指定して2次元配列のままにしておきます。空間隣接行列の列和と行和は同じ値になりますが、形状が異なります。
空間重み行列を作成します。
# 空間重み行列を作成 weight_mat = adj_mat / adj_mat.sum(axis=1, keepdims=True) np.nan_to_num(weight_mat, nan=0.0) # 列要素が全て0の場合用 print(weight_mat.round(2)) print(weight_mat.shape)
[[0. 0. 0. ... 0. 0. 0. ]
[0. 0. 0.2 ... 0. 0. 0. ]
[0. 0.25 0. ... 0. 0. 0. ]
...
[0. 0. 0. ... 0. 0.33 0. ]
[0. 0. 0. ... 0.33 0. 0.33]
[0. 0. 0. ... 0. 0.33 0. ]]
(72, 72)
空間隣接行列を行ごとに総和で割って基準化 します。
空間重み行列 を
weight_mat とします。こちらも、各要素は と
weight_mat[i-1, j-1] が対応します。
空間重み行列の列和と行和を確認します。
# 列和を計算 sum_col = weight_mat.sum(axis=0) print(sum_col[:5].round(2)) # 行和を計算 sum_row = weight_mat.sum(axis=1) print(sum_row[:5].round(2))
[0.82 1.06 1.01 1.49 0.59]
[1. 1. 1. 1. 1.]
行ごとに基準化しているので、行和が1になり、列和は(基本的に)1にならないのを確認できます。
以上で、空間重み行列行列を得られました。
空間重み行列の図
続いて、空間重み行列の定義や性質を可視化します。コロプレス図については本を参照してください。
空間重み行列のヒートマップを作成します。
# 地区数を取得 N = len(gdf_target) # 色の調整用 w_max = 1.0 # 空間重み行列を作図 fig, ax = plt.subplots(nrows=1, ncols=1, constrained_layout=True, figsize=(10, 10), dpi=100, facecolor='white') ax.pcolor(weight_mat, cmap='seismic', vmin=-w_max, vmax=w_max) # 行列 ax.set_xticks(ticks=np.arange(N)+0.5) ax.set_xticklabels(labels=gdf_target.city2, size=9, rotation=90) ax.set_yticks(ticks=np.arange(N)+0.5) ax.set_yticklabels(labels=gdf_target.city2, size=9) ax.invert_yaxis() # 軸の反転 ax.set_xlabel('city ( $j$ )') ax.set_ylabel('city ( $i$ )') ax.set_title(f'Osaka: $N = {len(gdf_target)}$', loc='left') fig.suptitle('spatial weight matrix', fontsize=20) plt.grid() ax.set_aspect('equal', adjustable='box') plt.show()
グラデーションの範囲の引数 vmin, vmax の中心が0になるように値を指定することで、隣接しない組み合わせのタイルに白色が配色されるように設定します。
隣接する区域の組み合わせごとに色付けし、重み の値の大小を色の濃淡で示します。空間重み行列は対称行列ではないのを確認できます。
空間重み行列と隣接ネットワークの対応をアニメーションで確認します。
作図コード(クリックで展開)
# 区域数を取得 N = len(gdf_target) # 色の調整用:(固定) w_max = 1.0 # 隣接区域以外を非表示化 weight_mat_masked = np.ma.masked_where(adj_mat==0, weight_mat) # グラフオブジェクトを初期化 fig, axes = plt.subplots(nrows=1, ncols=2, constrained_layout=True, figsize=(20, 10), width_ratios=[1, 1], dpi=100, facecolor='white') fig.suptitle('spatial weight matrix', fontsize=20) # 作図処理を定義 def update(n): # 前フレームのグラフを初期化 [ax.cla() for ax in axes] # 重みを格納 gdf_target['weight'] = weight_mat[:, n] # 隣接ネットワークを作図 ax = axes[0] gdf_district.boundary.plot(ax=ax, linewidth=0.5, edgecolor='black') # 行政区界 gdf_target.plot(ax=ax, column='weight', cmap='seismic', vmin=-w_max, vmax=w_max) # 各区域の値 for j in range(N): # 隣接区域のインデックスを抽出 adj_idx, = np.where(adj_mat[:, j] == 1) adj_idx = adj_idx[adj_idx > j] # 重複を除去 for i in adj_idx: ax.plot([gdf_target.centers.x[j], gdf_target.centers.x[i]], [gdf_target.centers.y[j], gdf_target.centers.y[i]], color='C0', linewidth=1.0) # 各区域-隣接区域 # 隣接区域のインデックスを抽出 adj_idx, = np.where(adj_mat[:, n] == 1) for i in adj_idx: ax.plot([gdf_target.centers.x[n], gdf_target.centers.x[i]], [gdf_target.centers.y[n], gdf_target.centers.y[i]], color='orange', linewidth=2.5) # 対象区域-隣接区域 for x, y, label in zip(gdf_target.centers.x, gdf_target.centers.y, gdf_target.city2_label): ax.annotate(text=label, xy=(x-0.015, y-0.005), size=9) # 区域名 ax.set_xlabel('longitude') ax.set_ylabel('latitude') ax.set_title(f'Osaka: $N = {N}$', loc='left') ax.grid() ax.set_aspect('equal', adjustable='box') # 対象区域以外・隣接区域以外を非表示化 target_mat_bool = np.tile(True, reps=adj_mat.shape) target_mat_bool[:, n] = False target_mat_masked = np.ma.masked_array(weight_mat_masked, target_mat_bool) # 空間重み行列を作図 ax = axes[1] ax.pcolor(weight_mat, cmap='seismic', vmin=-w_max, vmax=w_max) # 全区域 ax.pcolor(target_mat_masked, cmap='seismic', vmin=-w_max, vmax=w_max, color='gray') # 対象区域の隣接区域 ax.vlines(x=n+0.5, ymin=0, ymax=N, color='gray', linestyle='dashed') # 対象区域 ax.set_xticks(ticks=np.arange(N)+0.5) ax.set_xticklabels(labels=gdf_target.city2, size=9, rotation=90) ax.set_yticks(ticks=np.arange(N)+0.5) ax.set_yticklabels(labels=gdf_target.city2, size=9) ax.invert_yaxis() # 軸の反転 ax.set_xlabel('city ( $j$ )') ax.set_ylabel('city ( $i$ )') ax.set_title(f'{gdf_target.city2[n]}: $\\sum_{{i=1}}^N \\tilde{{w}}_{{ij}} = {np.sum(adj_mat[:, n])}$', loc='left') # 対象区域の隣接数 plt.grid() ax.set_aspect('equal', adjustable='box') # 動画を作成 ani = FuncAnimation(fig=fig, func=update, frames=N, interval=1000) # 動画を書出 ani.save( filename='ch2_2_weight_matrix.mp4', progress_callback = lambda i, n: print(f'frame: {i} / {n}') )
隣接区域のネットワーク図と空間重み行列のヒートマップ上で、区域ごとに、隣接区域をオレンジ色の線分、重みの大小を塗りつぶし色の濃淡で示します。隣接行列の列とネットワークが対応しているのを確認できます。
重みと隣接区域の隣接数の関係をアニメーションで確認します。
作図コード(クリックで展開)
# 隣接区域数を格納 gdf_target['count'] = np.sum(adj_mat, axis=1) # 区域数を取得 N = len(gdf_target) # 色の調整用 cnt_min, cnt_max = 0.0, gdf_target['count'].max() # 隣接区域以外を非表示化 weight_mat_masked = np.ma.masked_where(adj_mat==0, weight_mat) # グラフオブジェクトを初期化 fig, ax = plt.subplots(nrows=1, ncols=1, constrained_layout=True, figsize=(10, 10), dpi=100, facecolor='white') fig.suptitle('spatial weight matrix', fontsize=20) gdf_target.plot(ax=ax, column='count', cmap='jet', vmin=cnt_min, vmax=cnt_max, legend=True, legend_kwds={'label': 'count'}) # カラーバー表示用のダミー # 作図処理を定義 def update(n): # 前フレームのグラフを初期化 ax.cla() # 隣接ネットワークを作図 gdf_district.boundary.plot(ax=ax, linewidth=0.5, edgecolor='white') # 行政区界 gdf_target.plot(ax=ax, column='count', cmap='jet', vmin=cnt_min, vmax=cnt_max) # 隣接数 for i in np.where(adj_mat[:, n] == 1)[0]: ax.plot([gdf_target.centers.x[n], gdf_target.centers.x[i]], [gdf_target.centers.y[n], gdf_target.centers.y[i]], color='white', linewidth=6.0) # 白抜き for j in range(N): # 隣接区域のインデックスを抽出 adj_idx, = np.where(adj_mat[:, j] == 1) adj_idx = adj_idx[adj_idx > j] # 重複を除去 for i in adj_idx: ax.plot([gdf_target.centers.x[j], gdf_target.centers.x[i]], [gdf_target.centers.y[j], gdf_target.centers.y[i]], color='C0', linewidth=1.0) # 各区域-隣接区域 for i in np.where(adj_mat[:, n] == 1)[0]: ax.plot([gdf_target.centers.x[n], gdf_target.centers.x[i]], [gdf_target.centers.y[n], gdf_target.centers.y[i]], color=cm.seismic(weight_mat[i, n]*0.5+0.5), linewidth=2.5) # 対象区域-隣接区域 for x, y, label in zip(gdf_target.centers.x, gdf_target.centers.y, gdf_target.city2_label): ax.annotate(text=label, xy=(x-0.015, y-0.005), size=9) # 区域名 ax.set_xlabel('longitude') ax.set_ylabel('latitude') ax.set_title(f'Osaka: $N = {N}$', loc='left') ax.grid() ax.set_aspect('equal', adjustable='box') # 動画を作成 ani = FuncAnimation(fig=fig, func=update, frames=N, interval=1000) # 動画を書出 ani.save( filename='ch2_2_weight_matrix_net.mp4', progress_callback = lambda i, n: print(f'frame: {i} / {n}') )
隣接重み行列 weight_mat の各要素は 0 から 1 の値をとります。pcolor() の vmin, vmax 引数に -1, 1 を指定します。cm.seismic() は 0 から 1 の範囲で色付けされるので、weight_mat の値を 0.5 から 1 の値に変換して配色に使います。
線のグラデーションが見えにくいので、重みの線分の後に太めの白色の線を描画します。
区域ごとに隣接数で色付けし、重みの大小を線分の色の濃淡で示します。
隣接区域の隣接数が少ないほど重みが大きく、多いほど重みが小さくなるのを確認できます。
ここまでの記事では、空間重み行列を確認しました。次からの記事では、空間的自己相関を確認します。
参考文献
おわりに
図の試作をTwitterに流したときや1つ目の記事を投稿した段階で思いの外ウケてまして、正直なんで！？となってます。多少目を引く図かなとは思いますが、これまでの記事を比べて特別優れているとは思わないのですが。たまたまレコメンドアルゴリズムに乗ったんですかね。それとも私が知らないだけで空間統計って流行ってるんですか？
同様のクオリティ感で色々と書いているので他の記事も覗いてみてください。
【次の内容】
