からっぽのしょこ

はてなブログの仕様変更の影響で、数式の表示が崩壊しています。読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【Python】2.1:空間隣接行列の作図:境界の共有【空間データサイエンス入門のノート】

はじめに

 『Pythonで学ぶ空間データサイエンス入門』の独学時のまとめノートです。「導出編」「実装編」「可視化編」の三部構成でモデルやアルゴリズムの理解を目指します。
 本の内容から寄り道・回り道しながら進めます。本を読んだ上で補助的に読んでください。

 この記事では、境界の共有による空間隣接行列の定義について、プログラムと図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

2.1 空間隣接行列の作図:境界の共有

 空間隣接行列(spatial adjacency matrix)の定義をプログラムと図で確認します。この記事では、境界を共有している(実際に隣接している・contiguity)かにより隣接関係(adjacency relationship)を設定します。
 空間隣接行列については「2.1:空間隣接行列の定義式【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

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

# ライブラリを読込
import geopandas as gpd
from pysal.lib import weights
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap # 配色用
import matplotlib.patches as mpatches # 凡例用
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']]
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.1節を参照してください。

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

# 隣接関係を作成
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)
print(adj_mat.shape)
[[0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 1 1]
 [0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0]
 [0 1 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [0 1 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1]
 [0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 1 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 1]
 [0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1]
 [0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0]
 [1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0]
 [0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0]
 [1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0]
 [1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1]
 [0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 1 0 0]
 [0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0]
 [0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0]
 [1 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1]
 [1 0 0 1 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0]]
(24, 24)

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

 空間隣接行列の行和と列和を確認します。

# 行和を計算
sum_row = np.sum(adj_mat, axis=1)
print(sum_row)

# 列和を計算
sum_col = np.sum(adj_mat, axis=0)
print(sum_col)
[5 5 4 7 3 5 6 5 3 4 4 5 4 5 5 4 4 6 4 2 3 2 6 7]
[5 5 4 7 3 5 6 5 3 4 4 5 4 5 5 4 4 6 4 2 3 2 6 7]

 空間隣接行列  \mathbf{W} の行和  \sum_{j=1}^N w_{ij} と列和  \sum_{i=1}^N w_{ij} は一致します。

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

スポンサードリンク

グラフの作成

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

ヒートマップ

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

# カラーマップを作成
cmap = ListedColormap(colors=['orange']) # 単色を指定
# 区域数を取得
N = len(gdf_target)

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

# 非隣接区域を非表示化
adj_mat_masked = np.ma.masked_where(condition=adj_mat==0, a=adj_mat)

# 空間隣接行列を作図
fig, ax = plt.subplots(
    figsize=(7, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
ax.pcolormesh(
    adj_mat_masked, 
    cmap=cmap, #vmin=0.0, vmax=1.0, # (単色の設定用)
    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 adjacency matrix: contiguity', fontsize=20)
ax.grid()
ax.invert_yaxis() # (行番号との対応用)
ax.set_aspect(aspect='equal', adjustable='box')

# 凡例を追加
adj_patch = mpatches.Patch(
    color=cmap(1.0), 
    label='adjacent'
) # 凡例用のダミー
ax.legend(
    handles=[adj_patch], 
    bbox_to_anchor=(1.0, 1.0), loc='upper left'
) # 隣接関係

plt.show()

空間隣接行列のヒートマップ

 非隣接区域の要素をマスクしてタイルを非表示にすることで、隣接区域のみを色付けしています。
 非隣接区域のタイルを白色にすることで、2値の配色を設定する場合は、次のように処理します。

# カラーマップを作成
cmap = ListedColormap(colors=['white', 'orange']) # 2値色を指定
# 区域数を取得
N = len(gdf_target)

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

# 空間隣接行列を作図
fig, ax = plt.subplots(
    figsize=(7, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
ax.pcolormesh(
    adj_mat, 
    cmap=cmap, #vmin=0.0, vmax=1.0, # (2値色の設定用)
    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 adjacency matrix: contiguity', fontsize=20)
ax.grid()
ax.invert_yaxis() # (行番号との対応用)
ax.set_aspect(aspect='equal', adjustable='box')

# 凡例を追加
legend_handles = [
    mpatches.Patch(
        facecolor='white', #edgecolor='black', 
        label='not adjacent'
    ),
    mpatches.Patch(
        facecolor=cmap(1.0), 
        label='adjacent'
    )
] # 凡例用のダミー
ax.legend(
    handles=legend_handles, 
    bbox_to_anchor=(1.0, 1.0), loc='upper left'
) # 隣接関係

plt.show()

空間隣接行列のヒートマップ

 隣接関係のある区域(  w_{ij} = 1 のときの  i 行目・ j 列目のタイル)を色付けして示します。

 2次元のNumPy配列と ax.pcolormesh() を使って、ヒートマップを描画します。
 タイルはカラーコードなどで単色を指定できないので、カラーマップを作成して cmap 引数に指定します。カラーマップが1色で配列が1値の場合やカラーマップが2色で配列が2値の場合は vmin, vmax 引数の設定は不要です。
 各行が、各区域  i に関する隣接関係  \mathbf{w}_i = (w_1, \cdots, w_N) に対応します。

 空間隣接行列が対称行列であることを確認できます。

ネットワーク図

 隣接区域のネットワーク図を作成します。

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

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

# 隣接ネットワークを作図
fig, ax = plt.subplots(
    figsize=(8, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
gdf_target.boundary.plot(
    ax=ax, 
    edgecolor='black', linewidth=0.5
) # 行政区界
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('spatial adjacency graph: contiguity', fontsize=20)
ax.grid()
ax.set_aspect(aspect='equal', adjustable='box')
plt.show()

空間隣接行列のネットワーク図

 隣接区域(  w_{ij} = 1 のときの区域  i と区域  j )を線分で示します。

 GeoDataFrameと gdf.boundary.plot() を使って、境界線図を描画します。
 空間隣接行列と plt.plot() を使って、ネットワーク図を描画します。
 空間隣接行列 adj_mat から各行(対象の区域 i )の要素を取り出して、値が 1 の要素のインデックス(隣接区域の区域番号) jnp.where() で抽出します。
 地理空間データ gdf_target の重心座標の列 centroids から対象の区域 i と隣接区域 j のxy座標を取り出して、2点を線分で結びます。

 区域全体の隣接関係を確認できます。

コロプレス図

 隣接数を集計します。

# 隣接数を格納
#gdf_target['count'] = [len(v) for v in adj_obj.neighbors.values()] # 隣接関係の場合 
gdf_target['count'] = np.sum(adj_mat, axis=1) # 空間隣接行列の場合
print(gdf_target['count'].to_numpy())
[5 5 4 7 3 5 6 5 3 4 4 5 4 5 5 4 4 6 4 2 3 2 6 7]

 区域ごとの隣接数(空間隣接行列の行和)を集計して、データフレームに格納します。

 隣接数のコロプレス図を作成します。

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

# 配色を設定
u = 5.0
cnt_min = 0.0
cnt_max = gdf_target['count'].max()
cnt_max = np.ceil(cnt_max /u)*u # u単位で切り上げ

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

# 隣接数を作図
fig, ax = plt.subplots(
    figsize=(9, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
gdf_target.boundary.plot(
    ax=ax, 
    edgecolor='white', linewidth=0.5
) # 行政区界
gdf_target.plot(
    ax=ax, column='count', 
    cmap='jet', vmin=cnt_min, vmax=cnt_max, 
    legend=True, legend_kwds={'label': 'adjacency count', '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()

隣接数のコロプレス図

 各区域の隣接数を塗りつぶしのグラデーションで示します。

 GeoDataFrameと gdf.plot() を使って、コロプレス図(ポリゴン図・面図)を描画します。
 column 引数に隣接数の列を指定します。

 区域全体における隣接数の偏りを確認できます。各区域で交わる線分の数が隣接数に対応します。

ヒストグラム

 隣接数の度数を集計します。

# 隣接数を集計
adj_cnts = np.sum(adj_mat, axis=1)
print(adj_cnts)

# 度数軸の範囲を設定
u = 5.0
cnt_max = adj_cnts.max()
cnt_max = np.ceil(cnt_max /u)*u # u単位で切り上げ

# 度数を集計
cnt_vec  = np.arange(cnt_max+1, dtype=np.int16)
freq_vec = np.array([np.sum(adj_cnts == cnt) for cnt in cnt_vec])
print(cnt_vec)
print(freq_vec)
[5 5 4 7 3 5 6 5 3 4 4 5 4 5 5 4 4 6 4 2 3 2 6 7]
[ 0  1  2  3  4  5  6  7  8  9 10]
[0 0 2 3 7 7 3 2 0 0 0]

 区域ごとの隣接数を集計して adj_cnts とします。
 True1False0 として扱える(計算できる)ので、隣接数が cnt の度数(要素数・区域数)は np.sum(adj_cnts == cnt) で得られます。リスト内包表記により、cnt を1ずつ大きくしていき度数を求めて、配列に格納します。

 隣接数のヒストグラムを作成します。

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

# 隣接数を作図
fig, ax = plt.subplots(
    figsize=(8, 6), dpi=100, facecolor='white', 
    constrained_layout=True
)
ax.bar(
    x=cnt_vec, height=freq_vec, 
    width=1.0
) # 度数
ax.set_xlabel('adjacency count') # 隣接数軸
ax.set_ylabel('frequency') # 度数軸
ax.set_title(param_lbl, loc='left')
fig.suptitle('histogram: contiguity', fontsize=20)
ax.grid()
plt.show()

隣接数のヒストグラム

 横軸を隣接数、縦軸を度数のヒストグラムで隣接数の分布を示します。

 度数の配列と plt.bar() を使って、ヒストグラムを描画します。
 height 引数に度数の配列を指定します。

 隣接数の偏りを確認できます。

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

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

参考文献

おわりに

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

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

 この本を読んで、このブログを参考にして、わざわざこの図を手元で再現しようとするとは思わないのですが(なので初稿では作図の解説記事は書かなかったのですが)、詰まり所がいくつかあったので備忘録として独立した記事を書いておくことにしました。
 この後も作図編の記事がありますが、本の補足としては可視化編を読んでもらえば十分だと思います。

 最後に、つばきファクトリーのライブ映像をどうぞ♪


【次の内容】

 境界の共有による空間隣接行列を図のアニメーションで確認します。

www.anarchive-beta.com

 境界の共有による空間重み行列をプログラムと図で確認します。

https://www.anarchive-beta.com/entry/2026/07/10/210000www.anarchive-beta.com