“所有的路是否都通往罗马？”
美丽时尚的魅力：所有的路都通向罗马吗？
使用Python、网络科学和地理空间数据量化古老问题
我最近偶然发现了一个令人兴奋的数据集，标题为罗马道路网络（2008版），位于哈佛大学的Dataverse上：完美的GIS格式下的罗马帝国历史道路网络！此外，我一直在研究公共交通网络及如何识别这种网络科学中的热点和瓶颈。然后我很快意识到，通过将这些都结合起来，我可以迅速回答这个古老的问题，并看到罗马地区在过去有多么重要。
本文中的所有图片均由作者创建。
1. 读取和可视化数据
首先，让我们使用GeoPandas和Matplotlib快速加载和探索罗马道路网络数据。
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的交集。特别是，我按照该帖子中的方法从原始数据集中派生出了一个节点和边表。
# 创建边表edges = 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())
所以，在这里，我成功将GIS数据文件转换为了一个拥有5122个节点和7154条边的网络对象。现在，我想要查看一下。我们也可以使用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('意大利罗马')
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.optimizealpha(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')
最后，我们得出了一个令人欣慰的结论。如果按累计度数着色六边形单元格，罗马所占面积远远超过其他地区。如果按介数着色六边形，则情况相似 —— 罗马再次占主导地位。一个额外的发现是，连接罗马与中东的高速公路也显示为一个关键的、介数较高的部分。
简而言之，网络科学也证明了所有的道路都通往罗马！
