After enjoying the images, let’s go back to the graph and quantify it. Here, I will calculate the total degree of each node, by measuring the number of connections it has and the unnormalized betweenness centrality of each node, by counting the total number of shortest paths that cross each node.
node_degrees = dict(G_clean2.degree)
node_betweenness = dict(nx.betweenness_centrality(G_clean2, normalized = False))
Now I have the importance scores for each junction. Furthermore, in the nodes table, we also have its location; Now it’s time to move on to the main question. To do this, I quantify the importance of each node falling within the administrative boundaries of Rome. For this, I’ll need the administrative boundaries of Rome, which are relatively easy to obtain from OSMnx (note: the Rome of today is probably different from the Rome of yesteryear, but roughly, it should be fine).
admin = ox.geocode_to_gdf('Rome, Italy')
admin.plot()
The output of this cell:
Furthermore, visually, it is quite clear that Rome is not a single node in the road network; On the other hand, many are close. Therefore, we need some kind of spatial indexing, to help us group all the nodes and intersections of the road network that belong to Rome. Furthermore, it would be desirable for this aggregation to be comparable throughout the Empire. That’s why, instead of just mapping nodes in the Rome admin area, I’ll go with Uber. H3 Hexagonal grouping and creating hexagonal grids. Then, map each node to the surrounding hexagon and calculate the aggregate importance of that hexagon based on the centrality scores of the nodes in the closed network. Finally, I will discuss how the most central hexes overlap with Rome.
First, let’s get the administration area of the Roman Empire roughly:
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 divide the Empire polygon into a hexagonal grid:
import h3 # version: 3.7.3
from shapely.geometry import Polygon # version: 1.7.1
import numpy as np # version: 1.22.4def 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 to hexagons and attach the centrality scores to each hexagon. So. I add the importance of each node within each hexagon by adding its number of connections and the number of shortest paths that cross 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, combine the aggregated centrality scores with the Empire’s hexagonal map:
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. In this visual, I also add the empty grid as a basemap and then color each grid cell based on the total importance of the road network nodes it contains. In this way, the coloring will highlight the most critical cells in green. Additionally, I added the Rome polygon in white. First, colored by grade:
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: