After having fun with the visuals, let’s get back to the graph itself and quantify it. Here, I’ll compute the whole degree of every node, measuring the variety of connections it has, and the unnormalized betweenness centrality of every node, counting the whole variety of shortest paths crossing each node.

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

Now, I actually have the importance scores of every crossing. Moreover, within the nodes table, we even have their location — now it’s time to go for the principal query. For this, I quantify the importance of every node falling into the admin boundaries of Rome. For this, I’ll need the admin boundaries of Rome, which is comparatively easy to get from OSMnx (note: Rome today might be different from Rome back within the day, but approximatively, it needs to be high quality).

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

The output of this cell:

Also, on the visuals, it’s pretty clear that Rome shouldn’t be there as a single node within the road network; as a substitute, many are nearby. So we ned some type of binning, spatial indexing, which helps us group all road network nodes and intersections belonging to Rome. Moreover, this aggregation could be desired to be comparable across the Empire. That is why, as a substitute of just mapping nodes into the admin area of Rome, I’ll go for Uber’s H3 hexagon binning and create hexagon grids. Then, map each node into the enclosing hexagon and compute the aggregated importance of that hexagon based on the centrality scores of the enclosed network nodes. Finally, I’ll discuss how essentially the most central hexagons overlay with Rome.

First, let’s get the admin area of the Roman Empire in an approximative way:

import alphashape # version:  1.1.0
from descartes import PolygonPatch
# take a random sample of the node points
sample = nodes.sample(1000)
sample.plot()
# create its concave hull
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'))

The output of this cell:

Let’s split the Empire’s polygon right into a hexagon grid:

import h3 # version: 3.7.3
from shapely.geometry import Polygon # version: 1.7.1
import numpy as np # version: 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()

Result:

Now, map the road network nodes into hexagons and fasten the centrality scores to every hexagon. Then. I aggregate the importance of every node inside each hexagon by summing up their variety of connections and the variety of shortest paths crossing them:

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)

Finally, mix the aggregated centrality scores with the hexagon map of the Empire:

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

And visualize it. On this visual, I also add the empty grid as a base map after which color each grid cell based on the whole importance of the road network nodes inside. This fashion, the coloring will highlight essentially the most critical cells in green. Moreover, I added the polygon of Rome in white. First, coloured by degree:

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)

Result: