からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

トップ > 攻略ノート > 【Python】2.4：空間的自己相関(Local Moran's I)の可視化【空間データサイエンス入門のノート】
最終更新日

【Python】2.4：空間的自己相関(Local Moran's I)の可視化【空間データサイエンス入門のノート】

攻略ノート 攻略ノート-空間データサイエンス入門 空間分析 Python

はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学ノートです。本の内容から寄り道してアレコレ考えます。
　本を読んだ上で補助的に読んでください。 　

　この記事では、ローカル・モランのIについて、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

2.4 空間的自己相関(Local Moran's I)の可視化

　空間的自己相関(spatial autocorrelation)の指標であるローカル・モランのI(LMI・local Moran's I)を図で確認します。
　LMIの計算については「【Python】2.4：空間的自己相関(Local Moran's I)の実装【空間データサイエンス入門のノート】 - からっぽのしょこ」、空間重み行列(spatial weight matrix)については「【Python】2.2：空間重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import geopandas as gpd
from pysal.lib import weights
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib # 日本語の表示用
from mpl_toolkits.axes_grid1 import make_axes_locatable # カラーバーのサイズ調整用
import matplotlib.cm as cm # 配色の統一用
from matplotlib.animation import FuncAnimation


ローカル・モランのIの数式表現

　まずは、LMIの定義や計算を数式で確認します。

　 N 個の区域に関して、 i 番目の区域のデータ(変数)を x_i として、 N 個のデータをまとめて N 次元ベクトル \mathbf{x}、空間重み行列を (N \times N) の行列 \mathbf{W} とします。

\displaystyle \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix} ,\ \mathbf{W} = \begin{pmatrix} w_{11} & w_{12} & \cdots & w_{1N} \\ w_{21} & w_{22} & \cdots & w_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ w_{N1} & w_{N2} & \cdots & w_{NN} \end{pmatrix}

　LMIは、次の式で定義されます。

\displaystyle \mathrm{LMI}_i = \frac{ (x_i - \bar{x}) \sum_{j=1}^N w_{ij} (x_j - \bar{x}) }{ \frac{1}{N} \sum_{j=1}^N (x_j - \bar{x})^2 }

　ベクトルや行列を用いた計算については「2.3-4：空間的自己相関(Moran's I)の計算式の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

ローカル・モランのIの実装

　続いて、LMIの計算を行う関数を作成します。

　LMIの計算を自作関数として定義します。

# LMIの計算を実装
def culc_LMI(x, W):
    
    # データ数を取得
    N = len(x)
    
    # 偏差を計算
    x_dev = x - np.mean(x)

    # LMIを計算
    LMI = (x_dev * (W @ x_dev)) * (N-1) / (x_dev.T @ x_dev)
    return LMI.flatten()


　ここまでは、他の記事で確認しました。

データの読込

　可視化の前に、例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。

　「国土数値情報ダウンロードサイト」の行政区域データと地価公示データを読み込みます。

処理コード(クリックで展開)

# ファイルパスを指定
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

　市区町村(標準地行政区域コード)ごとに、対数公示価格の平均をとった列を作成して、ポリゴンデータと結合します。また、ポリゴンデータをまとめて、中心座標を求めます。


　以上で、作図に利用する空間データを用意できました。

ローカル・モランのIの計算

　次は、空間データからLMIを計算します。

　変数を設定します。

# 変数を抽出
x = gdf_target['log_price'].values
print(x[:5].round(2))
print(x.shape)

# 区域数を取得
N = len(x)
print(N)

[12.77 13.21 12.19 13.75 12.41]
(72,)
72

　この例では、区域ごとの平均対数公示価格を変数として用います。

　空間隣接行列を作成して、空間重み行列に変換します。

# 隣接関係を作成
wr = weights.Rook.from_dataframe(gdf_target)

# 空間隣接行列を作成
adj_mat = wr.full()[0].astype(int)
print(adj_mat)
print(adj_mat.shape)

# 空間重み行列を作成
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 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)
[[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)

　weights ライブラリを利用して、空間データから空間隣接行列を作成して、空間重み行列に変換します。

　LMIを計算します。

# 重みを設定
W = weight_mat.copy()

# LMIを計算
LMI = culc_LMI(x, W)
print(LMI[:5].round(3))

[1.494 1.865 0.439 2.9   0.562]

　自作関数を使ってLMIを計算します。

ローカル・モランのIの図

　続いて、LMIの定義や性質を可視化します。コロプレス図については本を参照してください。
　変数と偏差の関係については「【Python】2.3：空間ラグの可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　LMIのコロプレス図とヒートマップを作成します。

# LMIを格納
gdf_target['LMI'] = LMI

# 色の調整用
v_max = np.ceil(LMI.max())

# LMIを作図
fig, axes = plt.subplots(nrows=1, ncols=2, constrained_layout=True, 
                       figsize=(12, 10), width_ratios=[8, 2], dpi=100, facecolor='white')
fig.suptitle("Local Moran's I", fontsize=20)

# コロプレスマップを作図
ax = axes[0]
gdf_district.boundary.plot(ax=ax, linewidth=0.5, edgecolor='white') # 行政区界
gdf_target.plot(ax=ax, column='LMI', cmap='jet', vmin=-v_max, vmax=v_max) # 各区域の値
for label_x, label_y, label_str in zip(gdf_target.centers.x, gdf_target.centers.y, gdf_target.city2_label):
    ax.annotate(text=label_str, xy=(label_x-0.015, label_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')

# ヒートマップを作図
ax = axes[1]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='200%', pad=0.25)
c = ax.pcolor(LMI[:, np.newaxis], cmap='jet', vmin=-v_max, vmax=v_max) # ベクトル
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=gdf_target.city2)
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('city ( $i$ )')
ax.grid()
fig.colorbar(c, cax=cax, label='${LMI}_i$')
ax.set_aspect('equal', adjustable='box')

plt.show()

ローカル・モランのIのコロプレス図

　左のコロプレス図は、LMIの値に応じて区域ごとを塗りつぶして示します。右のヒートマップは、全ての区域を1列に並べた図です。

　各区域に隣接する区域の変数によってLMIが決まります。

　LMIの計算のヒートマップを作成します。

作図コード(クリックで展開)

# 偏差を計算
x_dev = (x - np.mean(x))[:, np.newaxis]

# 偏差を複製
X_dev = np.ones((N, 1)) @ x_dev.T

# SLを計算
sl = np.sum(W*X_dev, axis=1, keepdims=True)

# LMIを計算
tilde_lmi = x_dev * sl
lmi = tilde_lmi * (N-1) / np.sum(x_dev**2)

# 色の設定用
x_max  = np.ceil(max(x_dev.max(), abs(x_dev.min())))
w_max  = 1.0
wx_max = max((W*X_dev).max(), abs((W*X_dev).min()))
sl_max = max(sl.max(), abs(sl.min()))
tilde_lmi_max = max(tilde_lmi.max(), abs(tilde_lmi.min()))
lmi_max = max(lmi.max(), abs(lmi.min()))

# LMIの計算を作図
fig, axes = plt.subplots(nrows=2, ncols=5, constrained_layout=True, 
                       figsize=(31, 20), width_ratios=[3, 3, 3, 11, 11], 
                       dpi=100, facecolor='white')
fig.suptitle("Local Moran's I", fontsize=20)

# 偏差を作図
ax = axes[0, 2]
ax.pcolor(x_dev, cmap='jet', vmin=-x_max, vmax=x_max) # 偏差
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=gdf_target.city2, size=9) # 区域名
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('city ( $i$ )')
ax.set_title('$x - \\bar{x}$', fontsize=20)
ax.grid()
ax.set_aspect('equal', adjustable='box')

# 複製した偏差を作図
ax = axes[0, 3]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(X_dev, cmap='jet', vmin=-x_max, vmax=x_max) # 偏差
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels='')
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_title('$\\mathbf{1} (x-\\bar{x})^T$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$(x_j-\\bar{x})$')
ax.set_aspect('equal', adjustable='box')

# 空間重み行列を作図
ax = axes[0, 4]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(W, cmap='seismic', vmin=-w_max, vmax=w_max) # 重み
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title('$W$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$w_{ij}$')
ax.set_aspect('equal', adjustable='box')

# 重み付け偏差を作図
ax = axes[1, 3]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(W*X_dev, cmap='RdYlBu_r', vmin=-wx_max, vmax=wx_max) # 重み付け偏差
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title('$W \\odot \\mathbf{1} (x-\\bar{x})^T$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$w_{ij} (x_j-\\bar{x})$')
ax.set_aspect('equal', adjustable='box')

# 空間ラグを作図
ax = axes[1, 2]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='300%', pad=0.25)
c = ax.pcolor(sl, cmap='RdYlBu_r', vmin=-sl_max, vmax=sl_max) # 重み付け偏差の和
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('$i$')
ax.set_title('$W (x-\\bar{x})$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='${SL}_i = \\sum_{j=1}^N w_{ij} (x_j-\\bar{x})$')
ax.set_aspect('equal', adjustable='box')

# 空間ラグと偏差の積を作図
ax = axes[1, 1]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='300%', pad=0.25)
c = ax.pcolor(tilde_lmi, cmap='RdYlBu_r', vmin=-tilde_lmi_max, vmax=tilde_lmi_max) # 重み付け偏差の和と偏差の積
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('$i$')
ax.set_title('$(x-\\bar{x}) \\odot W (x-\\bar{x})$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$(x_i-\\bar{x}) \\sum_{j=1}^N w_{ij} (x_j-\\bar{x})$')
ax.set_aspect('equal', adjustable='box')

# ラベル用の文字列を作成
fnc_label  = '$(x-\\bar{x}) \\odot W (x-\\bar{x}) '
fnc_label += '\\frac{N-1}{\\sum_{j=1}^N (x_j-\\bar{x})^2}$'
elem_label  = '${LMI}_i = (x_i-\\bar{x}) \\sum_{j=1}^N w_{ij} (x_i-\\bar{x}) '
elem_label += '\\frac{N-1}{\\sum_{j=1}^N (x_j-\\bar{x})^2}$'

# LMIを作図
ax = axes[1, 0]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='300%', pad=0.25)
c = ax.pcolor(lmi, cmap='RdYlBu_r', vmin=-lmi_max, vmax=lmi_max) # 調整済み重み付け偏差の和と偏差の積
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('$i$')
ax.set_title(fnc_label, fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label=elem_label)
ax.set_aspect('equal', adjustable='box')

# 不要な図を非表示
axes[0, 0].axis('off')
axes[0, 1].axis('off')
axes[1, 4].axis('off')

plt.show()

ローカル・モランのIの計算

　それぞれ値の意味合いが異なるのでカラーマップ(配色)を変えて、値の大小をグラデーションで示します。
　偏差ベクトル(ベクトルとスカラの差)は、正確には \mathbf{x} - \bar{x} \mathbf{1} ですが、この図では(これ以上ゴチャゴチャしてもなんなので) \mathbf{x} - \bar{x} と表記しています。

　左上図は、縦ベクトルの偏差です。
　中上図は、 N 個(行)に複製した横ベクトル(転置した)の偏差です。
　右上図は、隣接重み行列です。
　中下図は、隣接重み行列によって「重み付けした偏差」(中上図と右上図の要素ごとの積)です。隣接しない組み合わせは0になりLMIに影響せず、隣接数が少ない区域との組み合わせほど値が大きく(割引率が小さく)なります。
　左下の右図は、区域ごとの「重み付け偏差の和」(中下図の行ごとの和)であり、区域ごとの空間ラグ(SL)です。
　中図は、区域ごとの「重み付け偏差の和」と「偏差」の積(右図と左上図の要素ごとの積)です。
　左図は、不偏分散で調整した区域ごとの「重み付け偏差の和」であり、区域ごとのLMIです。

隣接関係とLMIの関係

　簡単な例として、格子状に隣接する区域における変数とLMIの関係を図で確認します。

作図コード(クリックで展開)

　区域数を設定します。

# 一辺の区域数を指定
n = 5

# 全体の区域数を計算
N = n**2
print(N)

25

　格子状に並ぶ区域の一辺の区域数 n を指定して、全体の区域数 N を計算します。

　空間隣接行列を作成します。

# 空間隣接行列を作成
tilde_W = np.zeros((N, N), dtype='int')
adj_i = np.arange(n, N)
adj_j = np.arange(0, N-n)
tilde_W[adj_i, adj_j] = 1 # 上:(i, i-n)
tilde_W[adj_j, adj_i] = 1 # 下:(i, i+n)
adj_i = np.hstack([np.arange(n*i+1, n*(i+1)) for i in range(n)])
adj_j = np.hstack([np.arange(n*i, n*(i+1)-1) for i in range(n)])
tilde_W[adj_i, adj_j] = 1 # 左:(i, i-1)
tilde_W[adj_j, adj_i] = 1 # 右:(i, i+1)
print(tilde_W[:10, :10])

[[0 1 0 0 0 1 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0]
 [0 1 0 1 0 0 0 1 0 0]
 [0 0 1 0 1 0 0 0 1 0]
 [0 0 0 1 0 0 0 0 0 1]
 [1 0 0 0 0 0 1 0 0 0]
 [0 1 0 0 0 1 0 1 0 0]
 [0 0 1 0 0 0 1 0 1 0]
 [0 0 0 1 0 0 0 1 0 1]
 [0 0 0 0 1 0 0 0 1 0]]

　格子状に並ぶ区域に対応する空間隣接行列を作成します。
　一辺の区域数を n とすると、区域 i の上の区域番号は i-n、下の区域番号は i+n、左の区域番号は i-1、右の区域番号は i+1 です。ただし、区域 i の場所によって、上下左右のどの区域が存在するかは異なります。

　空間重み行列を作成します。

# 隣接重み行列に変換
W = tilde_W / tilde_W.sum(axis=1, keepdims=True)
np.nan_to_num(W, 0.0) # 列要素が全て0の場合用
print(W[:10, :10].round(2))

[[0.   0.5  0.   0.   0.   0.5  0.   0.   0.   0.  ]
 [0.33 0.   0.33 0.   0.   0.   0.33 0.   0.   0.  ]
 [0.   0.33 0.   0.33 0.   0.   0.   0.33 0.   0.  ]
 [0.   0.   0.33 0.   0.33 0.   0.   0.   0.33 0.  ]
 [0.   0.   0.   0.5  0.   0.   0.   0.   0.   0.5 ]
 [0.33 0.   0.   0.   0.   0.   0.33 0.   0.   0.  ]
 [0.   0.25 0.   0.   0.   0.25 0.   0.25 0.   0.  ]
 [0.   0.   0.25 0.   0.   0.   0.25 0.   0.25 0.  ]
 [0.   0.   0.   0.25 0.   0.   0.   0.25 0.   0.25]
 [0.   0.   0.   0.   0.33 0.   0.   0.   0.33 0.  ]]

　空間隣接行列を行ごとに基準化(行和で割る)して、空間重み行列を作成します。

　フレームごとの変数の変化を設定します。

# 変数の値を指定
x_val = 10.0

# 追加する区域番号を指定
if n%2 == 1: # 区域数が奇数の場合
    var_add_i  = [np.arange(0, N, step=2)] # 互い違いのインデックス
    var_add_i += [np.arange(1, N, step=2)] # 隙間のインデックス
    var_add_i  = np.hstack(var_add_i)
else: # 区域数が偶数の場合
    var_add_i  = [(n*i)+np.arange(i%2, n, step=2) for i in range(n)] # 互い違いのインデックス
    var_add_i += [(n*i)+np.arange((i+1)%2, n, step=2) for i in range(n)] # 隙間のインデックス
    var_add_i  = np.array(var_add_i)
print(var_add_i)

# 除去する区域番号を指定
var_del_i = np.arange(N)
print(var_del_i)

[ 0  2  4  6  8 10 12 14 16 18 20 22 24  1  3  5  7  9 11 13 15 17 19 21
 23]
[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 24]

　この例では、全ての区域が値を持たない(値が 0 の)状態を初期値として、値を持つ区域を隣接しないように増やしていき、続いて全ての区域が値を持つようにします。さらに、順番に値を消していきます。
　各区域の変数の値を x_val として指定しておきます。

　LMIの計算のヒートマップのアニメーションを作成します。

# フレーム数を取得
frame_num = len(np.hstack([var_add_i, var_del_i]))

# 色の設定用
x_max     = abs(x_val)
x_dev_max = abs(x_val)
w_max   = 1.0
wx_max  = 4.5
sl_max  = 8.5
lmi_max = 5.0
_wx_max_lt, _sl_max_lt, _lmi_max_lt = [], [], [] # 確認用

# ラベル用の文字列を作成
fnc_label  = '$(x-\\bar{x}) \\odot W (x-\\bar{x}) '
fnc_label += '\\frac{N-1}{\\sum_{j=1}^N (x_j-\\bar{x})^2}$'
elem_label  = '${LMI}_i = (x_i-\\bar{x}) \\sum_{j=1}^N w_{ij} (x_i-\\bar{x}) '
elem_label += '\\frac{N-1}{\\sum_{j=1}^N (x_j-\\bar{x})^2}$'

# サイズ調整用の値を指定
d = 0.6

# 変数を初期化
x = np.zeros(N)
city_mat = np.tile(np.nan, reps=(n+2, n+2))

# グラフオブジェクトを初期化
fig, axes = plt.subplots(nrows=2, ncols=4, constrained_layout=True, 
                         figsize=((12+4+4+12)*d, (10+10)*d), width_ratios=[10, 3, 3, 10], 
                         dpi=100, facecolor='white')
fig.suptitle("Local Moran's I", fontsize=20)

# 作図処理を定義
def update(frame_i):
    
    # 前フレームのグラフを初期化
    [ax.cla() for ax in axes.flatten()]

    # 0除算の回避用
    eps = np.random.normal(loc=0.0, scale=0.001, size=1).item()

    # 変数を変更
    if frame_i < N:
        # 値を追加
        var_i     = var_add_i[frame_i]
        x[var_i]  = x_val
        x[var_i] += eps
    else:
        # 値を除去
        var_i    = var_del_i[frame_i-N]
        x[var_i] = eps

    # 格子状の区域に変形
    x_mat = x.reshape((n, n))
    city_mat[1:(n+1), 1:(n+1)] = x_mat
    
    # 偏差を計算
    x_dev = (x - np.mean(x))[:, np.newaxis]
    
    # 差を複製
    X_dev = np.ones((N, 1)) @ x_dev.T
    
    # SLを計算
    sl = np.sum(W*X_dev, axis=1, keepdims=True)
    
    # LMIを計算
    lmi = x_dev * sl * (N-1) / np.sum(x_dev**2)
    
    # 色設定用の最大値の確認用
    _wx_max_lt.append(abs((W*X_dev)).max())
    _sl_max_lt.append(abs(sl).max())
    _lmi_max_lt.append(abs(lmi).max())
    print('weight x dev max:', max(_wx_max_lt))
    print('SL max:          ', max(_sl_max_lt))
    print('LMI max:         ', max(_lmi_max_lt))
    
    # 区域を作図
    ax = axes[0, 0]
    ax.pcolor(city_mat, cmap='jet', vmin=-x_max, vmax=x_max) # 変数
    ax.text(x=1, y=1+1, s='${LMI}_i = $', 
            size=10, ha='right', va='bottom') # 枠内に収まらなかった対策用
    city_i = 0
    for label_i in range(1, n+1):
        for label_j in range(1, n+1):
            lmi_val = lmi[city_i, 0]
            city_i += 1
            ax.text(x=label_j+0.5, y=label_i+0.5, s=f'$i = {city_i}$', 
                    size=10, ha='center', va='center') # 区域番号
            ax.text(x=label_j, y=label_i+1.0, s=f'${lmi_val:.3f}$', 
                    size=10, ha='left', va='bottom', 
                    bbox={'facecolor': cm.RdYlBu_r(lmi_val/lmi_max*0.5+0.5), 
                          'edgecolor': 'white', 'pad': 0.0}) # LMI
    ax.set_xticks(ticks=np.arange(1, n+1)+0.5)
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(1, n+1)+0.5)
    ax.set_yticklabels(labels='')
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('longitude')
    ax.set_ylabel('latitude')
    ax.set_title(f'$\\bar{{x}} = {np.mean(x):.2f}$', loc='left')
    ax.grid()
    ax.set_aspect('equal', adjustable='box')
    
    # 変数を作図
    ax = axes[0, 1]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='100%', pad=0.25)
    c = ax.pcolor(x[:, np.newaxis], cmap='jet', vmin=-x_max, vmax=x_max) # 変数
    ax.set_xticks(ticks=[0.5])
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_ylabel('$i$')
    ax.set_title('$x$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$x_i$')
    ax.set_aspect('equal', adjustable='box')
    
    # 偏差を作図
    ax = axes[0, 2]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='100%', pad=0.25)
    c = ax.pcolor(x_dev, cmap='jet', vmin=-x_dev_max, vmax=x_dev_max) # 変数
    ax.set_xticks(ticks=[0.5])
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_ylabel('$i$')
    ax.set_title('$x - \\bar{x}$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$x_i - \\bar{x}$')
    ax.set_aspect('equal', adjustable='box')
    
    # 複製した偏差を作図
    ax = axes[0, 3]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size=f'{1/N*100}%', pad=0.25)
    c = ax.pcolor(X_dev, cmap='jet', vmin=-x_dev_max, vmax=x_dev_max) # 偏差
    ax.set_xticks(ticks=np.arange(N)+0.5)
    ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels='')
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('$j$')
    ax.set_title('$\\mathbf{1} (x-\\bar{x})^T$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$(x_j-\\bar{x})$')
    ax.set_aspect('equal', adjustable='box')
    
    # 空間重み行列を作図
    ax = axes[1, 0]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size=f'{1/N*100}%', pad=0.25)
    c = ax.pcolor(W, cmap='seismic', vmin=-w_max, vmax=w_max) # 重み
    ax.set_xticks(ticks=np.arange(N)+0.5)
    ax.set_xticklabels(labels=np.arange(N)+1)
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('$j$')
    ax.set_ylabel('$i$')
    ax.set_title('$W$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$(x_i-\\bar{x}) (x_j-\\bar{x})$')
    ax.set_aspect('equal', adjustable='box')
    
    # 重み付け偏差を作図
    ax = axes[1, 3]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size=f'{1/N*100}%', pad=0.25)
    c = ax.pcolor(W*X_dev, cmap='RdYlBu_r', vmin=-wx_max, vmax=wx_max) # 重み付け偏差
    ax.set_xticks(ticks=np.arange(N)+0.5)
    ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1, size=8)
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('$j$')
    ax.set_ylabel('$i$')
    ax.set_title('$W \\odot \\mathbf{1} (x-\\bar{x})^T$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$w_{ij} (x_j-\\bar{x})$')
    ax.set_aspect('equal', adjustable='box')

    # 空間ラグを作図
    ax = axes[1, 2]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='100%', pad=0.25)
    c = ax.pcolor(sl, cmap='RdYlBu_r', vmin=-sl_max, vmax=sl_max) # 重み付け偏差の和
    ax.set_xticks(ticks=[0.5])
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1, size=8)
    ax.invert_yaxis() # 軸の反転
    ax.set_ylabel('$i$')
    ax.set_title('$W (x-\\bar{x})$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='${SL}_i = \\sum_{j=1}^N w_{ij} (x_j-\\bar{x})$')
    ax.set_aspect('equal', adjustable='box')
    
    # LMIを作図
    ax = axes[1, 1]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='100%', pad=0.25)
    c = ax.pcolor(lmi, cmap='RdYlBu_r', vmin=-lmi_max, vmax=lmi_max) # 調整済み重み付け偏差の和と偏差の積
    ax.set_xticks(ticks=[0.5])
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1, size=8)
    ax.invert_yaxis() # 軸の反転
    ax.set_ylabel('$i$')
    ax.set_title(fnc_label, fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label=elem_label)
    ax.set_aspect('equal', adjustable='box')

# 動画を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=500)

# 動画を書出
ani.save(
    filename='ch2_4_lmi_adj.mp4', 
    progress_callback = lambda i, n: print(f'frame: {i} / {n}')
)

　一定の値を持つ区域の数や位置を変化させて、LMIへの影響を確認します。

　左上図は、簡易的な表現のコロプレス図(地図)です(行列を可視化した図ではありません)。上下左右の区域に隣接しています。区域ごとにLMIも表示しており、中下の右図と対応しています。
　中上図は、変数と偏差のベクトルです。変数の値によって平均や偏差の値に影響します。
　右上図は、複製した横ベクトルの(転置した)偏差です。
　左下図は、隣接重み行列です。
　右下図は、隣接重み行列によって重み付けした横ベクトルの偏差です。
　中下の右図は、重み付け偏差の和であり、区域ごとの空間ラグ(SL)です。
　左図は、調整した重み付け偏差の和であり、区域ごとのLMIです。

　この章では、空間的自己相関を確認しました。次の章では、線形回帰モデルを確認します。

参考文献

Pythonで学ぶ 空間データサイエンス入門: 地域の特徴を発見する方法

Pythonで学ぶ 空間データサイエンス入門: 地域の特徴を発見する方法

Amazon

おわりに

　書く前はいつも1つ書けば他はコピペでサクッとできるだろうと思いながら始めるのですが、中々そうもいかずそれ相応に大変な目に合うんですよねぇ。
　これが表現力の限界……でも言葉も合わせるとギリギリ伝わる範囲であ……れ。

【次の内容】

www.anarchive-beta.com