Tutte le strade portano a Roma?

La bellezza e la moda si incontrano a Roma il tuo viaggio verso lo stile perfetto

Quantificare l’antica domanda con Python, Network Science e Dati Geospaziali

Di recente sono venuto a conoscenza di un eccitante dataset intitolato Roman Road Network (versione 2008) su Harvard’s Dataverse: le reti stradali storiche dell’Impero Romano in un perfetto formato GIS! Inoltre, sto lavorando su un progetto sulle reti di trasporto pubblico e su come identificare i punti caldi e i punti di strozzatura di una tale rete tramite la Network Science. Poi mi sono reso conto che mettendo tutto insieme, avrei potuto rispondere rapidamente a questa antica domanda e vedere quanto fosse centrale l’area di Roma ai tempi.

In questo articolo, tutte le immagini sono state create dall’autore.

1. Leggi e visualizza i dati

Prima di tutto, utilizziamo GeoPandas e Matplotlib per caricare ed esplorare rapidamente i dati della rete stradale romana.

• Cinque cose che GenAI può e non può fare
• Classificazione equa con debiasing avversario
• Ricercatori di Google DeepMind presentano Promptbreeder un sistema di intelligenza artificiale auto-referenziale e auto-migliorante che può automaticamente evolvere promemoria efficaci specifici del dominio in un dato dominio.

import geopandas as gpd # version: 0.9.0
import matplotlib.pyplot as plt # version: 3.7.1
gdf = gpd.read_file('dataverse_files-2')
gdf = gdf.to_crs(4326)
print(len(gdf))
gdf.head(3)

L’output di questa cella:

Ora visualizziamolo:

f, ax = plt.subplots(1,1,figsize=(15,10))
gdf.plot(column = 'CERTAINTY', ax=ax)

2. Trasforma la rete stradale in un oggetto graph

La figura precedente mostra la rete stradale come un insieme di poligoni di linee. Tuttavia, per poter quantificare l’importanza di, per esempio, Roma, ho intenzione di effettuare alcuni calcoli grafici. Ciò significa che devo trasformare queste linee in un grafico.

Il pacchetto OSMNx è perfetto per questo – all’intersezione degli strumenti di dati geospaziali e della famosa libreria di analisi dei grafi, NetworkX. In particolare, ho seguito questo thread per derivare una tabella di nodi e una tabella di edge dal dataset originale.

# crea una tabella edge
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)

Risultato di questa cella:

import pandas as pd # versione: 1.4.2
from shapely.geometry import Point # versione: 1.7.1
# crea una tabella dei nodi
nodes = pd.DataFrame(edges_copy ['u'] .append (edges_copy ['v']), colonne = ['osmid'])
nodes ['geometry'] = [Point (float (n.split ('_') [0]), float (n.split ('_') [1])) per n in nodi.osmid.to_list ()]
nodes ['x'] = [float (n.split ('_') [0]) per n in nodi.osmid.to_list ()]
nodes ['y'] = [float (n.split ('_') [1]) per n in nodi.osmid.to_list ()]
nodes = gpd.GeoDataFrame(nodi)
nodes = nodi.set_index ('osmid')
nodes.head (3)

Risultato di questa cella:

Crea il grafico:

import osmnx as ox # versione: 1.0.1 # Ora crea il grafografo_attributi = {'crs': 'epsg: 4326', 'semplificato': vero} G = ox.graph_from_gdfs(nodi, archi[['geometria']], grafico_attributi) print (type (G)) print (G.number_of_nodes ()), print (G.number_of_edges ()) 

Quindi qui, sono riuscito a trasformare il file di dati GIS in un oggetto di rete con 5122 nodi e 7154 archi. Ora, vorrei dare un’occhiata. È possibile visualizzare le reti anche con NetworkX. Tuttavia, preferisco optare per il software open-source Gephi. Offre maggiore flessibilità e migliori opzioni di regolazione visiva. Trasformiamo G in un file compatibile con Gephi e esportiamolo, in questa versione lavorerò con un grafico non pesato e non orientato.

# Trasforma ed esporta il graficonetworkx come nx # versione: 2.5G_pulito = nx.Graph () per u, v, dati in G.edges (data = True): G_pulito.add_edge (u, v) G_clean2 = nx.Graph () G_clean2.add_edges_from (G_pulito.edges (data = True)) nx.write_gexf (G_clean2, 'roman_empire_network.gexf') 

Inoltre, creo una tabella dati chiamata coordiantes.csv in cui salvo le coordinate di ogni nodo (incrocio) della rete stradale.

nodi2 = nodi[nodi.index.isin(set(G.nodes))].drop(columns = ['geometria'])nodi2.index.name = 'Id'nodi2.to_csv('coordinate.csv') 

3. Visualizzazione della rete in Gephi

Il modo esatto per visualizzare le reti in Gephi merita una sessione a sé stante, quindi qui mostrerò invece il risultato.

In questa visualizzazione, ogni nodo corrisponde a un incrocio, il colore codifica le cosiddette comunità di rete (sottografi densamente interconnessi), mentre i nodi sono dimensionati in base alla loro centralità tra nodi. La centralità tra nodi è una misura della centralità di rete che quantifica il ruolo di ponte di un nodo. Pertanto, maggiore è un nodo, più centrale è.

Sulla visualizzazione, è anche interessante vedere come la geografia influisce sui cluster e come sorprendentemente l’Italia si distingue, probabilmente a causa della rete stradale interna molto più densa.

4. Centralità della rete

Dopo aver apprezzato le immagini, torniamo al grafico stesso e lo quantifichiamo. Qui, calcolerò il grado totale di ciascun nodo, misurando il numero di connessioni che ha, e la centralità non normalizzata tra nodi di ciascun nodo, conteggiando il numero totale di percorsi più brevi che attraversano ciascun nodo.

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

Ora ho i punteggi di importanza di ogni incrocio. Inoltre, nella tabella dei nodi, abbiamo anche la loro posizione – adesso è il momento di affrontare la domanda principale. Per fare questo, quantifico l’importanza di ogni nodo che ricade all’interno dei confini amministrativi di Roma. Per questo, avrò bisogno dei confini amministrativi di Roma, che è relativamente facile da ottenere da OSMnx (nota: Roma oggi è probabilmente diversa da Roma nel passato, ma approssimativamente dovrebbe essere sufficiente).

admin = ox.geocode_to_gdf('Roma, Italia')
admin.plot()

L’output di questa cella:

Inoltre, dalle immagini, è abbastanza evidente che Roma non è presente come singolo nodo nella rete stradale; invece, ce ne sono molti nelle vicinanze. Quindi abbiamo bisogno di una sorta di raggruppamento, di indicizzazione spaziale, che ci aiuti a raggruppare tutti i nodi della rete stradale e gli incroci appartenenti a Roma. Inoltre, sarebbe auspicabile che questa aggregazione sia comparabile in tutto l’Impero. Ecco perché, invece di semplicemente mappare i nodi nell’area amministrativa di Roma, opterò per il raggruppamento esagonale di Uber’s H3 e creerò griglie esagonali. Quindi, mappo ogni nodo nell’esagono che lo racchiude e calcolo l’importanza aggregata di quell’esagono in base ai punteggi di centralità dei nodi della rete inclusi. Infine, discuterò come i centroidi esagonali più centrali si sovrappongono a Roma.

Prima di tutto, otteniamo l’area amministrativa dell’Impero Romano in modo approssimativo:

import alphashape # versione: 1.1.0
from descartes import PolygonPatch

# prendi un campione casuale dei punti nodo
sample = nodes.sample(1000)
sample.plot()

# crea il suo inviluppo concavo
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'))

L’output di questa cella:

Dividiamo il poligono dell’Impero in una griglia esagonale:

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

Risultato:

Ora, mappa i nodi della rete stradale negli esagoni e associ i punteggi di centralità a ciascun esagono. Successivamente, aggrego l’importanza di ciascun nodo all’interno di ciascun esagono sommando il numero delle loro connessioni e il numero di percorsi più brevi che li attraversano:

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)

Infine, combinare i punteggi di centralità aggregati con la mappa esagonale dell’Impero:

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

E visualizzarlo. In questa visualizzazione, aggiungo anche la griglia vuota come base mappa e coloro ogni cella della griglia in base all’importanza complessiva dei nodi della rete stradale al suo interno. In questo modo, la colorazione evidenzierà le celle più critiche in verde. Inoltre, ho aggiunto il poligono di Roma in bianco. Prima, colorato per gradi:

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)

Risultato:

Ora, colorato per intermediarietà:

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

Infine, arriviamo a una conclusione rassicurante. Se coloriamo le celle esagonali in base al loro grado cumulativo, l’area di Roma prevale di gran lunga. Se coloriamo gli esagoni in base all’intermediarietà, l’immagine è simile – Roma prevale nuovamente. Un punto aggiuntivo qui è che l’autostrada che collega Roma con il Medio Oriente emerge anche come un segmento critico ad alta intermediarietà.

tl;dr anche according to la scienza delle reti, tutte le strade portavano a Roma!