「すべての道はローマに通じるのですか？」

「すべての美容とファッションの道はローマに通じるのですか？」

Python、ネットワークサイエンス、および地理空間データによる古代の問いを量化する

最近、ハーバードのデータベースである ローマ帝国の歴史的道路ネットワーク（2008年版） という興奮するデータセットに出会いました。 また、公共交通ネットワークについてのプロジェクトに取り組んでおり、このようなネットワークサイエンスのホットスポットとボトルネックを特定する方法を探っています。 その後、これらをすべて組み合わせることによって、この古代の問いに迅速に答えることができ、当時のローマの中心地がいかに重要であったかを見ることができることに気付きました。

この記事では、すべての画像は著者によって作成されました。

1. データの読み込みと可視化

まず、GeoPandasとMatplotlibを使用して、ローマ道路ネットワークデータを迅速に読み込み、探索します。

• 「ジェンAIができることとできないことの5つのこと」
• イメージセグメンテーション：詳細ガイド
• 『日常のデザイン（AI）』

import geopandas as gpd # version: 0.9.0import matplotlib.pyplot as plt # version: 3.7.1gdf = gpd.read_file('dataverse_files-2')gdf = gdf.to_crs(4326)print(len(gdf))gdf.head(3)

このセルの出力:

それを可視化する:

f、ax = plt.subplots(1,1,figsize=(15,10))gdf.plot(column = 'CERTAINTY', ax=ax)

2. 道路ネットワークをグラフオブジェクトに変換する

前の図は道路ネットワークを一連の線ポリゴンとして表示しています。 ただし、例えばローマの重要性を量化するためには、いくつかのグラフ計算を行う必要があります。 これは、これらのラインストリングをグラフに変換する必要があるということです。

OSMNxパッケージは、地理空間データツールと有名なグラフ解析ライブラリであるNetworkXとの交差点に位置しています。 特に、元のデータセットからノードテーブルとエッジテーブルを導出するために このスレッド に従いました。

# create an edge tableedges = gdf.copy()edges['u'] = [str(g.coords[0][0]) + '_' + str(g.coords[0][1]) for g in edges.geometry.to_list()]edges['v'] = [str(g.coords[-1][0]) + '_' + str(g.coords[-1][1]) for g in edges.geometry.to_list()]edges_copy = edges.copy()edges['key'] = range(len(edges))edges = edges.set_index(['u', 'v', 'key'])edges.head(3)

このセルの結果:

import pandas as pd # バージョン：1.4.2
from shapely.geometry import Point # バージョン：1.7.1

# ノードテーブルを作成
nodes = pd.DataFrame(edges_copy['u'].append(edges_copy['v']), columns=['osmid'])
nodes['geometry'] = [Point(float(n.split('_')[0]), float(n.split('_')[1])) for n in nodes.osmid.to_list()]
nodes['x'] = [float(n.split('_')[0]) for n in nodes.osmid.to_list()]
nodes['y'] = [float(n.split('_')[1]) for n in nodes.osmid.to_list()]
nodes = gpd.GeoDataFrame(nodes)
nodes = nodes.set_index('osmid')
nodes.head(3)

このセルの結果:

グラフを作成:

import osmnx as ox # バージョン：1.0.1

# グラフを構築
graph_attrs = {'crs': 'epsg:4326', 'simplified': True}
G = ox.graph_from_gdfs(nodes, edges[['geometry']], graph_attrs)
print(type(G))
print(G.number_of_nodes()), print(G.number_of_edges())

ここで、5122のノードと7154のエッジを持つネットワークオブジェクトにGISデータファイルを変換することができました。さて、見てみましょう。NetworkXでもネットワークを可視化することができますが、私はオープンソースのソフトウェアGephiをおすすめします。より柔軟性があり、視覚の微調整オプションも優れています。GをGephi互換のファイルに変換してエクスポートしましょう。このバージョンでは、重みのない無向グラフで作業します。

# グラフを変換してエクスポート
import networkx as nx # バージョン：2.5
G_clean = nx.Graph()
for u, v, data in G.edges(data=True):
    G_clean.add_edge(u, v)
G_clean2 = nx.Graph()
G_clean2.add_edges_from(G_clean.edges(data=True))
nx.write_gexf(G_clean2, 'roman_empire_network.gexf')

さらに、道路ネットワークの各ノード（交差点）の座標を保存するためのデータテーブル「coordinates.csv」を作成します。

nodes2 = nodes[nodes.index.isin(set(G.nodes))].drop(columns=['geometry'])
nodes2.index.name = 'Id'
nodes2.to_csv('coordinates.csv')

3. Gephiでネットワークを可視化する

Gephiでネットワークを可視化する方法についての詳細な方法は独立したセッションで説明する価値がありますが、ここでは結果を示します。

この可視化では、各ノードは交差点に対応し、色はネットワークのコミュニティ（密結合なサブグラフ）をエンコードし、ノードのサイズはその媒介中心性に応じて変化します。媒介中心性は、ノードの架け橋の役割を数量化するネットワークの中心性尺度です。したがって、ノードが大きいほど中心的です。

ビジュアルでは、地理がクラスターを形成する様子や、驚くべきことにイタリアが独自の位置に立っていることも興味深いです。おそらく、より密な内部道路ネットワークのためです。

4. ネットワークの中心性

ビジュアルを楽しんだ後は、グラフ自体に戻って数量化します。ここでは、各ノードの総次数（接続数）と各ノードの非正規化媒介中心性（ノードを横断する最短経路の総数を数える）を計算します。

node_degrees = dict(G_clean2.degree)
node_betweenness = dict(nx.betweenness_centrality(G_clean2, normalized=False))

今、それぞれの交差点の重要度スコアを持っています。ノードのテーブルには、それらの位置もあります。そして、これからメインの質問に進みます。それには、ローマの行政区域に含まれる各ノードの重要度を数量化します。そのために、ローマの行政区域を取得する必要があります。OSMnxから比較的簡単に取得できます（注意：現在のローマは昔とは異なる可能性がありますが、大体は問題ありません）。

admin = ox.geocode_to_gdf('Rome, Italy')
admin.plot()

このセルの出力：

また、視覚的には、ローマは道路ネットワークの単一のノードとして存在していないことがわかります。代わりに、多くのノードが近くにあります。そのため、私たちはヒンビン（空間インデックス）が必要です。これにより、ローマに属するすべての道路ネットワークのノードと交差点をグループ化できます。さらに、この集計は帝国全体で比較可能であることが望まれます。そのため、ただノードをローマの行政区域にマッピングするだけではなく、UberのH3のヘキサゴンビンニングを使用してヘキサゴングリッドを作成します。そして、各ノードを囲むヘキサゴンにマッピングし、そのヘキサゴンの中に含まれるネットワークノードの中心性スコアに基づいて集計された重要度を計算します。最後に、最も中心的なヘキサゴンがローマとオーバーレイする方法について説明します。

まず、ローマ帝国の行政区域をおおまかに取得しましょう：

import alphashape # バージョン：1.1.0
from descartes import PolygonPatch
# ノードポイントのランダムなサンプルを取得
sample = nodes.sample(1000)
sample.plot()
# 凹型の凸包を作成
points = [(point.x, point.y) for point in sample.geometry]
alpha = 0.95 * alphashape.optimize_alpha(points)
hull = alphashape.alphashape(points, alpha)
hull_pts = hull.exterior.coords.xy
fig, ax = plt.subplots()
ax.scatter(hull_pts[0], hull_pts[1], color='red')
ax.add_patch(PolygonPatch(hull, fill=False, color='green'))

このセルの出力：

次に、ローマ帝国のポリゴンをヘキサゴングリッドで分割しましょう：

import h3 # バージョン：3.7.3
from shapely.geometry import Polygon # バージョン：1.7.1
import numpy as np # バージョン：1.22.4
def split_admin_boundary_to_hexagons(polygon, resolution):
    coords = list(polygon.exterior.coords)
    admin_geojson = {"type": "Polygon", "coordinates": [coords]}
    hexagons = h3.polyfill(admin_geojson, resolution, geo_json_conformant=True)
    hexagon_geometries = {hex_id: Polygon(h3.h3_to_geo_boundary(hex_id, geo_json=True)) for hex_id in hexagons}
    return gpd.GeoDataFrame(hexagon_geometries.items(), columns=['hex_id', 'geometry'])
roman_empire = split_admin_boundary_to_hexagons(hull, 3)
roman_empire.plot()

結果：

次に、道路ネットワークのノードをヘキサゴンにマッピングし、各ヘキサゴンに中心性スコアを付けます。そして、各ヘキサゴン内の各ノードの重要度を集計し、接続数や最短経路のクロス数を合計します：

gdf_merged = gpd.sjoin(roman_empire, nodes[['geometry']])gdf_merged['degree'] = gdf_merged.index_right.map(node_degrees)gdf_merged['betweenness'] = gdf_merged.index_right.map(node_betweenness)gdf_merged = gdf_merged.groupby(by = 'hex_id')[['degree', 'betweenness']].sum()gdf_merged.head(3)

最後に、集計された中心性スコアを帝国の六角形地図に組み合わせます：

roman_empire = roman_empire.merge(gdf_merged, left_on = 'hex_id', right_index = True, how = 'outer')roman_empire = roman_empire.fillna(0)

そして、可視化します。このビジュアルでは、空のグリッドをベースマップとして追加し、その後、グリッドセルの総合的な道路ネットワークノードの重要度に基づいて各グリッドセルの色を変えます。これにより、色付けによって最も重要なセルを緑で強調します。また、ローマのポリゴンを白い色に追加しました。まず、次のように度数によって色分けします：

f, ax = plt.subplots(1,1,figsize=(15,15))gpd.GeoDataFrame([hull], columns = ['geometry']).plot(ax=ax, color = 'grey', edgecolor = 'k', linewidth = 3, alpha = 0.1)roman_empire.plot(column = 'degree', cmap = 'RdYlGn', ax = ax)gdf.plot(ax=ax, color = 'k', linewidth = 0.5, alpha = 0.5)admin.plot(ax=ax, color = 'w', linewidth = 3, edgecolor = 'w')ax.axis('off')plt.savefig('degree.png', dpi = 200)

結果：

次に、以下のように媒介中心性に基づいて色分けします：

f, ax = plt.subplots(1,1,figsize=(15,15))gpd.GeoDataFrame([hull], columns = ['geometry']).plot(ax=ax, color = 'grey', edgecolor = 'k', linewidth = 3, alpha = 0.1)roman_empire.plot(column = 'betweenness', cmap = 'RdYlGn', ax = ax)gdf.plot(ax=ax, color = 'k', linewidth = 0.5, alpha = 0.5)admin.plot(ax=ax, color = 'w', linewidth = 3, edgecolor = 'w')ax.axis('off')plt.savefig('betweenness.png', dpi = 200, bbox_inches = 'tight')

最後に、安心感を持って結論付けられます。累積度数に基づいて六角形セルを色分けすると、ローマの領域が圧倒的に勝ちます。六角形を媒介中心性に基づいて色分けすると、同様の結果となります。ここでの追加要素は、ローマと中東を結ぶ高速道路も重要であり、媒介中心性が高いセグメントとして浮かび上がります。

要約すると、ネットワーク科学によれば、全ての道はローマに至ると言えます！

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!