# ライブラリを読込
import geopandas as gpd
from pysal.lib import weights
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap # 配色用
import japanize_matplotlib # 日本語の描画用

データの読込

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

国土数値情報ダウンロードサイトのデータの場合

　この例では、国土数値情報ダウンロードサイトから2026年版(基準年)の大阪府の行政区域データをダウンロードして、大阪市のデータを利用します。

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

# ファイルパスを指定
DISTRICT_PATH = 'data/N03-20260101_27_GML/N03-20260101_27.shp' # ポリゴンデータ:大阪府(2026年版)

# データを読込
gdf_district = gpd.read_file(DISTRICT_PATH, encoding='UTF-8') # (2026年版の場合)

# 行政区域データを取得
gdf_district = gdf_district[['N03_004', 'N03_005', 'N03_007', 'geometry']] # (2026年版の場合)
gdf_district.columns = ['city1', 'city2', 'd_code', 'geometry']

# データを整形
dst_proj     = 6668
gdf_district = gdf_district.to_crs(epsg=dst_proj) # 空間座標系を再設定

　事前に保存しておいた行政区域のポリゴンデータを読み込みます。
　必要に応じて扱いやすいようにデータを整形しておきます。

　行政区域データをまとめます。

# データを統合
gdf_target = gdf_district.dissolve(by='d_code', as_index=False) # 飛び地を統合
gdf_target['centroids'] = gdf_target['geometry'].centroid # 重心座標

# データフレームを整形
gdf_target = gdf_target.reindex(
    columns=['city1', 'city2', 'd_code', 'geometry', 'centroids']
) # (目視での確認用)

　市区町村(行政区域コード)ごとに、飛び地などのポリゴンデータをまとめて、重心座標を求めます。(この処理だと警告文が出ますが、ここでの目的には影響しない誤差なので無視します。)

　大阪市のデータを取り出します。

# 地域を指定
city_name = '大阪市'

# データを抽出 
gdf_target = gdf_target[gdf_target['city1'] == city_name]
gdf_target.head()

	city1	city2	d_code	geometry	centroids
0	大阪市	都島区	27102	POLYGON ((135.51481 34.72057, 135.51489 34.720...	POINT (135.52721 34.71193)
1	大阪市	福島区	27103	POLYGON ((135.46247 34.6874, 135.46253 34.6879...	POINT (135.47534 34.6939)
2	大阪市	此花区	27104	MULTIPOLYGON (((135.34359 34.62112, 135.34357 ...	POINT (135.42011 34.668)
3	大阪市	西区	27106	POLYGON ((135.46667 34.67284, 135.4642 34.6747...	POINT (135.48375 34.67776)
4	大阪市	港区	27107	MULTIPOLYGON (((135.43006 34.64979, 135.43006 ...	POINT (135.45084 34.66032)

　大阪府全域だとゴチャゴチャして分かりにくいので、大阪市内の区のデータを利用することにします。

　以上で、空間重み行列の作成や作図に利用する地理空間データを用意できました。

スポンサードリンク

空間重み行列の作成

　次は、地理空間データから空間重み行列を作成します。
　詳しくは本の2.2節を参照してください。

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

# 隣接関係を作成
adj_obj = weights.Rook.from_dataframe(
    df=gdf_target, geom_col='geometry', use_index=False
) # ルーク型
'''
adj_obj = weights.Queen.from_dataframe(
    df=gdf_target, geom_col='geometry', use_index=False
) # クイーン型
'''
# 空間隣接行列を作成
adj_mat, area_ids = adj_obj.full()
adj_mat = adj_mat.astype(dtype=np.int8)
print(adj_mat[:10, :10])
print(adj_mat.shape)

[[0 0 0 0 0 0 0 0 0 1]
 [0 0 1 1 0 0 0 0 1 0]
 [0 1 0 1 1 0 0 0 1 0]
 [0 1 1 0 1 1 0 1 0 0]
 [0 0 1 1 0 1 0 0 0 0]
 [0 0 0 1 1 0 0 1 0 0]
 [0 0 0 0 0 0 0 1 0 0]
 [0 0 0 1 0 1 1 0 0 0]
 [0 1 1 0 0 0 0 0 0 0]
 [1 0 0 0 0 0 0 0 0 0]]
(24, 24)

　weights ライブラリを利用して、行政区域のポリゴンデータから隣接関係を作成して、空間隣接行列に変換します。
　空間隣接行列 $\tilde{\mathbf{W}}$ を adj_mat とします。ただし、Pythonではインデックスを0から数えるので、各要素 $\tilde{w}_{ij}$ は adj_mat[i-1, j-1] と対応します。

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

# 空間重み行列を作成
weight_mat = adj_mat / np.sum(adj_mat, axis=1, keepdims=True) # 正規化
np.nan_to_num(weight_mat, nan=0.0) # (隣接数0の対策用)　
print(weight_mat[:10, :10].round(2))
print(weight_mat.shape)

[[0.   0.   0.   0.   0.   0.   0.   0.   0.   0.2 ]
 [0.   0.   0.2  0.2  0.   0.   0.   0.   0.2  0.  ]
 [0.   0.25 0.   0.25 0.25 0.   0.   0.   0.25 0.  ]
 [0.   0.14 0.14 0.   0.14 0.14 0.   0.14 0.   0.  ]
 [0.   0.   0.33 0.33 0.   0.33 0.   0.   0.   0.  ]
 [0.   0.   0.   0.2  0.2  0.   0.   0.2  0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.17 0.   0.  ]
 [0.   0.   0.   0.2  0.   0.2  0.2  0.   0.   0.  ]
 [0.   0.33 0.33 0.   0.   0.   0.   0.   0.   0.  ]
 [0.25 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]]
(24, 24)

　空間隣接行列 $\tilde{\mathbf{W}}$ を正規化(規準化) $w_{ij} = \frac{\tilde{w}_{ij}}{\sum_{j'=1}^N \tilde{w}_{ij'}}$ して、空間重み行列 $\mathbf{W}$ を計算します。
　行和の計算時に keepdims=True を指定して、計算結果をN行1列の2次元配列にすることで、ブロードキャストにより2次元配列の全体を行和で正規化できます。
　空間重み行列 $\mathbf{W}$ を weight_mat とします。こちらも、各要素 $w_{ij}$ は weight_mat[i-1, j-1] と対応します。
　または、.transform 属性を変更して、空間重み行列に変換します。

# 空間重み行列を作成
adj_obj.transform = 'R' # 正規化
weight_mat, _ = adj_obj.full()
print(weight_mat[:10, :10].round(2))
print(weight_mat.shape)

[[0.   0.   0.   0.   0.   0.   0.   0.   0.   0.2 ]
 [0.   0.   0.2  0.2  0.   0.   0.   0.   0.2  0.  ]
 [0.   0.25 0.   0.25 0.25 0.   0.   0.   0.25 0.  ]
 [0.   0.14 0.14 0.   0.14 0.14 0.   0.14 0.   0.  ]
 [0.   0.   0.33 0.33 0.   0.33 0.   0.   0.   0.  ]
 [0.   0.   0.   0.2  0.2  0.   0.   0.2  0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.17 0.   0.  ]
 [0.   0.   0.   0.2  0.   0.2  0.2  0.   0.   0.  ]
 [0.   0.33 0.33 0.   0.   0.   0.   0.   0.   0.  ]
 [0.25 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]]
(24, 24)

　.transform 属性に 'O' を指定すると空間隣接行列、'R' を指定すると空間重み行列を取り出せます。

　空間重み行列の行和と列和を確認します。

# 行和を計算
sum_row = np.sum(weight_mat, axis=1)
print(sum_row[:10].round(2))

# 列和を計算
sum_col = np.sum(weight_mat, axis=0)
print(sum_col[:10].round(2))

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1.01 1.14 1.01 1.49 0.59 1.18 1.16 0.82 0.7  0.87]

　行ごとに正規化しているので、空間重み行列 $\mathbf{W}$ の行和 $\sum_{i=1}^N w_{ij}$ は1になり、列和 $\sum_{j=1}^N w_{ij}$ は(基本的に)1になりません。

　以上で、空間重み行列を得られました。

スポンサードリンク

グラフの作成

　続いて、空間重み行列の図を作成します。
　コロプレス図については本を参照してください。

ヒートマップ

　空間重み行列のヒートマップを作成します。

# カラーマップを作成
cmap = LinearSegmentedColormap.from_list(
    name='white_red',
    colors=['white', 'red']
)

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

# 配色を設定
w_min, w_max = 0.0, 1.0 # 最小値・最大値

# ラベルを作成
param_lbl = f'$N = {N}$'

# 空間重み行列を作図
fig, ax = plt.subplots(
    figsize=(7, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
pc = ax.pcolormesh(
    weight_mat, 
    cmap=cmap, vmin=w_min, vmax=w_max, 
    shading='auto'
) # 全区域
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=gdf_target['city2'], size=10, rotation=90) # 区域名
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=gdf_target['city2'], size=10) # 区域名
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title(param_lbl, loc='left')
fig.suptitle('spatial weight matrix: contiguity', fontsize=20)
fig.colorbar(
    mappable=pc, ax=ax, shrink=1.0, 
    label='spatial weight'
) # 重み軸
ax.grid()
ax.invert_yaxis() # (行番号との対応用)
ax.set_aspect(aspect='equal', adjustable='box')
plt.show()

　区域の組み合わせごとの重み $w_{ij}$ に応じて濃淡を付けて示します。

　2次元のNumPy配列と ax.pcolormesh() を使って、ヒートマップを描画します。
　非隣接区域(重みが0)のタイルを白色にするために、(最小値が白色のカラーマップがないため、)カラーマップを作成して cmap 引数に指定します。

　空間隣接行列とは異なり、空間重み行列が対称行列でないことを確認できます。
　各行が、各区域 $i$ に関する重み $\mathbf{w}_i = (w_1, \cdots, w_N)$ に対応します。区域 $i$ の隣接数(隣接する区域の数)を $N_i$ とすると、行ごとに重みは $0$ か $\frac{1}{N_1}$ の2値 $w_{ij} \in \{0, \frac{1}{N_i}\}$ をとるので、行ごとの重みは同じ値になります。

コロプレス図

　区域を指定します。

# 区域を指定
area_idx = 22

# 重みを格納
gdf_target['weight'] = weight_mat[area_idx]
print(gdf_target['weight'].to_numpy()[:10].round(2))

[0.17 0.17 0.   0.17 0.   0.   0.   0.   0.   0.17]

　可視化の対象とする区域番号 $i$ を指定します。ただし、Pythonではインデックスを0から数えるので、 $i-1$ の値を area_idx とします。
　対象の区域に関する重み(空間重み行列の $i$ 行)をデータフレームに格納します。

　指定した区域の重みのコロプレス図を作成します。

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

# 配色を設定
w_min, w_max = 0.0, 1.0 # 最小値・最大値

# ラベルを作成
param_lbl = f'$N = {N}, i = {area_idx+1}$'

# 重みを作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
gdf_target.boundary.plot(
    ax=ax, 
    edgecolor='black', linewidth=0.5
) # 行政区界
gdf_target.plot(
    ax=ax, column='weight', 
    cmap=cmap, vmin=w_min, vmax=w_max, 
    legend=True, legend_kwds={'label': 'adjacency weight', 'shrink': 1.0}
) # 重み
for i in range(N):
    adj_idx, = np.where(adj_mat[i] == 1) # 隣接区域のインデックス
    adj_idx  = adj_idx[adj_idx > i]      # 重複を除去
    for j in adj_idx:
        O_x, O_y = gdf_target.loc[i, 'centroids'].coords[0] # 各区域の座標
        P_x, P_y = gdf_target.loc[j, 'centroids'].coords[0] # 隣接区域の座標
        ax.plot(
            [O_x, P_x], 
            [O_y, P_y], 
            color='C0', linewidth=2.0
        ) # 各区域 - 隣接区域
for x, y, area_lbl in zip(gdf_target['centroids'].x, gdf_target['centroids'].y, gdf_target['city2']):
    ax.text(
        x=x, y=y, s=area_lbl, ha='center', va='center', 
        size=10
    ) # 区域名
ax.set_xlabel('longitude')
ax.set_ylabel('latitude')
ax.set_title(param_lbl, loc='left')
fig.suptitle('choropleth map: contiguity', fontsize=20)
ax.grid()
ax.set_aspect(aspect='equal', adjustable='box')
plt.show()

　可視化対象の区域 $i$ について、 $N$ 個の区域 $j = 1, \dots, N$ それぞれに対する重み $w_{ij}$ に応じて濃淡を付けて示します。この図は大阪市北区の例です。
　また、隣接区域( $w_{ij} = 1$ のときの区域 $i$ と区域 $j$ )を線分で示します。

　GeoDataFrameと gdf.boundary.plot() を使って、境界線図を描画します。
　GeoDataFrameと gdf.plot() を使って、コロプレス図(ポリゴン図・面図)を描画します。
　column 引数に重みの列を指定します。
　ネットワーク図については「空間隣接行列の作図」を参照してください。

　区域 $i$ の隣接数(隣接する区域の数)を $N_i$ とすると、区域 $i$ と隣接関係のある区域の重みは $w_{ij} = \frac{1}{N_i}$ 、隣接関係のない区域(区域 $i$ 自体の重みも含む)の重みは $w_{ij} = 0$ になるのを確認できます。
　1つの区域に関する重みの図なので、区域全体の重みは確認できません。

　以上で、空間重み行列の図を作成できました。

　この記事では、境界の共有による空間重み行列をプログラムと図で確認しました。次の記事では、図のアニメーションで確認します。

参考文献

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

作者:植杉　大
Independently published

Amazon

おわりに

2026.07.09：加筆修正の際に「空間重み行列の可視化」から記事を分割しました。

　図の種類を増やして、作図処理の解説も少し整えました。

　2冊目『はじめての地理空間データサイエンス』と合わせて1冊目『空間データサイエンス入門』の加筆修正を進めたことで、空間隣接行列の記事も含めた空間重み行列編の記事が充実しました。
　どれもコピペみたいな記事になっていますが、手法ごとに微妙に図が変えており、その差異を表現するための試行錯誤が思いの外に面倒臭いんです。そして、その細かな処理の違いが後々混同することになるので、備忘録としても書き残しておく必要が出てくるんです。
　折角なので、他の記事の図と見比べて、手法ごとの差分から(全ての手法を覚えるわけではなく)目的の手法の理解を深めてください。

　記事の投稿日に公開されたBEYOOOOONDSのライブ映像をどうぞ♪