¿Todos los caminos conducen a Roma?  Cuantificando la antigua cuestión con… |  de Milan Janosov |  octubre de 2023

Después de disfrutar de las imágenes, volvamos al gráfico y cuantifíquelo. Aquí, calcularé el grado total de cada nodo, midiendo el número de conexiones que tiene y la centralidad de intermediación no normalizada de cada nodo, contando el número total de caminos más cortos que cruzan cada nodo.

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

Ahora tengo las puntuaciones de importancia de cada cruce. Además, en el nodos tabla, también tenemos su ubicación; ahora es el momento de pasar a la pregunta principal. Para ello, cuantifico la importancia de que cada nodo caiga dentro de los límites administrativos de Roma. Para esto, necesitaré los límites administrativos de Roma, que son relativamente fáciles de obtener desde OSMnx (nota: la Roma actual probablemente sea diferente de la Roma de antaño, pero aproximadamente, debería estar bien).

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

La salida de esta celda:

Los límites administrativos de Roma.

Además, en lo visual, está bastante claro que Roma no es un nodo único en la red de carreteras; en cambio, muchos están cerca. Por lo tanto, necesitamos algún tipo de indexación espacial, que nos ayude a agrupar todos los nodos e intersecciones de la red de carreteras que pertenecen a Roma. Además, sería deseable que esta agregación fuera comparable en todo el Imperio. Es por eso que, en lugar de simplemente mapear nodos en el área de administración de Roma, optaré por Uber. H3 agrupación hexagonal y crear cuadrículas hexagonales. Luego, asigne cada nodo al hexágono circundante y calcule la importancia agregada de ese hexágono en función de las puntuaciones de centralidad de los nodos de la red cerrada. Finalmente, discutiré cómo los hexágonos más centrales se superponen con Roma.

Primero, obtengamos el área de administración del Imperio Romano de forma aproximada:

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'))

La salida de esta celda:

Un subconjunto de nodos de red y el casco cóncavo circundante.

Dividamos el polígono del Imperio en una cuadrícula hexagonal:

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()

Resultado:

La rejilla hexagonal del Imperio Romano.

Ahora, asigne los nodos de la red de carreteras a hexágonos y adjunte las puntuaciones de centralidad a cada hexágono. Entonces. Agrego la importancia de cada nodo dentro de cada hexágono sumando su número de conexiones y el número de caminos más cortos que las cruzan:

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)
Vista previa de la tabla de cuadrícula hexagonal agregada.

Finalmente, combine las puntuaciones de centralidad agregadas con el mapa hexagonal del Imperio:

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

Y visualizarlo. En este objeto visual, también agrego la cuadrícula vacía como mapa base y luego coloreo cada celda de la cuadrícula según la importancia total de los nodos de la red de carreteras que contiene. De esta forma, la coloración resaltará en verde las celdas más críticas. Además, agregué el polígono de Roma en blanco. Primero, coloreado por grado:

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)

Resultado: