Do All the Roads Lead to Rome?. Quantifying the Ancient Question with… | by Milan Janosov | Oct, 2023

After having fun with the visuals, let’s get again to the graph itself and quantify it. Right here, I’ll compute the whole diploma 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 every node.

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

Now, I’ve the significance scores of every crossing. Moreover, within the nodes desk, we even have their location — now it’s time to go for the primary query. For this, I quantify the significance of every node falling into the admin boundaries of Rome. For this, I’ll want the admin boundaries of Rome, which is comparatively straightforward to get from OSMnx (word: Rome immediately might be totally different from Rome again within the day, however approximatively, it needs to be superb).

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

The output of this cell:

Additionally, on the visuals, it’s fairly clear that Rome isn’t there as a single node within the highway community; as an alternative, many are close by. So we ned some form of binning, spatial indexing, which helps us group all highway community nodes and intersections belonging to Rome. Moreover, this aggregation can be desired to be comparable throughout the Empire. Because of this, as an alternative of simply mapping nodes into the admin space of Rome, I’ll go for Uber’s H3 hexagon binning and create hexagon grids. Then, map every node into the enclosing hexagon and compute the aggregated significance of that hexagon primarily based on the centrality scores of the enclosed community nodes. Lastly, I’ll focus on how essentially the most central hexagons overlay with Rome.

First, let’s get the admin space of the Roman Empire in an approximative approach:

import alphashape # model:  1.1.0
from descartes import PolygonPatch
# take a random pattern of the node factors
pattern = nodes.pattern(1000)
pattern.plot()
# create its concave hull
factors = [(point.x, point.y) for point in sample.geometry]
alpha = 0.95 * alphashape.optimizealpha(factors)
hull = alphashape.alphashape(factors, alpha)
hull_pts = hull.exterior.coords.xy
fig, ax = plt.subplots()
ax.scatter(hull_pts[0], hull_pts[1], coloration='purple')
ax.add_patch(PolygonPatch(hull, fill=False, coloration='inexperienced'))

The output of this cell:

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

import h3 # model: 3.7.3
from shapely.geometry import Polygon # model: 1.7.1
import numpy as np # model: 1.22.4
def split_admin_boundary_to_hexagons(polygon, decision):
coords = record(polygon.exterior.coords)
admin_geojson = {"kind": "Polygon",  "coordinates": [coords]}
hexagons = h3.polyfill(admin_geojson, decision, 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.gadgets(), columns = ['hex_id', 'geometry'])
roman_empire = split_admin_boundary_to_hexagons(hull, 3)
roman_empire.plot()

End result:

Now, map the highway community nodes into hexagons and fasten the centrality scores to every hexagon. Then. I combination the significance of every node inside every 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)

Lastly, 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 visible, I additionally add the empty grid as a base map after which coloration every grid cell primarily based on the whole significance of the highway community nodes inside. This manner, the coloring will spotlight essentially the most vital cells in inexperienced. Moreover, I added the polygon of Rome in white. First, coloured by diploma:

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

End result: