Questo articolo presenta un approfondimento tecnico su un modello computazionale, implementato in Python, che inferisce i flussi di mobilità partendo proprio da questi dati di presenza aggregata. Il nostro dataset di riferimento è un’estrazione relativa alla giornata del 31 Maggio 2025 per l’area di Ravenna.
Il percorso che descriviamo è un’evoluzione. Un primo approccio intuitivo, basato su un confronto “brute-force” di ogni cella con tutte le altre, si è rivelato computazionalmente insostenibile, con tempi di esecuzione stimati in giorni o settimane. Questo ha reso necessario un cambio di paradigma. Il codice e la metodologia che seguono rappresentano la soluzione ottimizzata e metodologicamente robusta, basata su tre pilastri:
- Efficienza algoritmica: Sostituzione della ricerca “brute-force” con un approccio basato su indici spaziali, riducendo la complessità da quadratica (O(N²)) a quasi lineare (O(N*k)).
- Robustezza statistica: Introduzione di filtri per eliminare il rumore di fondo e concentrarsi solo su eventi di mobilità significativi.
- Accuratezza geospaziale: Utilizzo di librerie geospaziali (
GeoPandas,scikit-mobility) e di una “tassellazione” precisa per ancorare l’analisi alla geografia reale (passaggio dai tiles del db a tasselli per scikit-mob).
Vediamo ora in dettaglio il codice che implementa questa soluzione avanzata.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')
# Import scikit-mobility components
try:
import skmob
from skmob import FlowDataFrame
print("✅ scikit-mobility importato con successo")
print(f"Versione: {skmob.__version__}")
# Verifichiamo le misure collettive disponibili
try:
from skmob.measures.collective import origin_destination_flows
print("✅ Misure collettive disponibili")
except ImportError:
print("⚠️ Misure collettive non disponibili, implementeremo manualmente")
origin_destination_flows = None
# Verifichiamo i modelli disponibili
try:
from skmob.models.gravity import Gravity
from skmob.models.radiation import Radiation
print("✅ Modelli Gravity e Radiation disponibili")
except ImportError:
print("⚠️ Modelli Gravity e Radiation non disponibili, continueremo senza")
Gravity = None
Radiation = None
except ImportError as e:
print(f"❌ Errore nell'importazione di scikit-mobility: {e}")
exit(1)
# Import per tessellazione
try:
import geopandas as gpd
from shapely.geometry import Polygon
from shapely import wkt
print("✅ GeoPandas e Shapely disponibili per tessellazione")
except ImportError:
print("❌ GeoPandas non disponibile, installare con: pip install geopandas")
exit(1)
# ==================================
# --- PARAMETRI DEL MODELLO ---
# ==================================
# Filtri sui dati
MIN_PRESENCE_THRESHOLD = 5 # Ignora celle con meno di 5 persone
MIN_FLOW_THRESHOLD = 10 # Inferisce flussi solo per movimenti > 10 persone
# Parametri di inferenza del flusso
# 1 grado di latitudine/longitudine è circa 111 km. 150m sono circa 0.00135 gradi.
# Raggio per pedoni (1.25 km in 15 min) -> 1250m / (111000 m/grado) = ~0.011 gradi
# Raggio per auto (5 km in 15 min) -> 5000m / (111000 m/grado) = ~0.045 gradi
SEARCH_RADIUS_DEGREES = 0.045 # Raggio di ricerca in gradi decimali (scegli quello adatto al tuo caso d'uso)
DISTANCE_PENALTY_FACTOR = 100 # Aumenta la "penalità" per la distanza nella distribuzione del flusso
def load_and_prepare_traffic_data():
"""
Carica e prepara i dati di traffico per analisi di flusso ottimizzata
"""
print("Caricamento dati di traffico per analisi di flusso ottimizzata...")
try:
# Carichiamo i dati delle tessere
tiles_df = pd.read_csv('superset_data/datasets/pixel/tiles_d.csv')
# Carichiamo i dati di traffico
traffic_df = pd.read_csv('superset_data/datasets/pixel/daily_cells/2025-05-31.csv')
# Prepariamo i dati
traffic_df['TimeStamp'] = pd.to_datetime(traffic_df['TimeStamp'], errors='coerce')
traffic_df['dateInsert'] = pd.to_datetime(traffic_df['dateInsert'], errors='coerce')
traffic_df = traffic_df.dropna(subset=['TimeStamp', 'dateInsert'])
traffic_df['timestamp_corretto'] = traffic_df['dateInsert'] - timedelta(hours=1)
# Uniamo con le coordinate
traffic_with_coords = pd.merge(traffic_df, tiles_df[['TileX', 'TileY', 'Lat', 'Lon']],
on=['TileX', 'TileY'], how='inner')
print(f"Dati di mobilità caricati: {len(traffic_with_coords)} righe")
return traffic_with_coords
except Exception as e:
print(f"Errore nel caricamento dei dati di mobilità: {e}")
return None
def create_temporal_snapshots_optimized(traffic_data):
"""
Crea snapshot temporali ottimizzati per inferire i flussi
"""
print("Creazione snapshot temporali ottimizzati...")
# Raggruppiamo per timestamp e cella
snapshots = traffic_data.groupby(['timestamp_corretto', 'TileX', 'TileY', 'Lat', 'Lon']).agg({
'P': 'sum'
}).reset_index()
# Creiamo fasce orarie per analisi
snapshots['hour'] = snapshots['timestamp_corretto'].dt.hour
snapshots['time_slot'] = snapshots['timestamp_corretto'].dt.floor('H')
# Filtriamo solo le celle con presenza significativa
snapshots = snapshots[snapshots['P'] > MIN_PRESENCE_THRESHOLD]
print(f"Snapshot ottimizzati: {len(snapshots)} punti temporali")
print(f"Fasce orarie: {snapshots['hour'].min()}-{snapshots['hour'].max()}")
print(f"Celle uniche: {snapshots[['TileX', 'TileY']].drop_duplicates().shape[0]}")
return snapshots
def find_nearby_cells(cell_x, cell_y, all_cells_df, max_distance=SEARCH_RADIUS_DEGREES):
"""
Trova celle vicine a una cella data (ottimizzazione geografica)
"""
# Coordinate della cella di riferimento
cell_data = all_cells_df[(all_cells_df['TileX'] == cell_x) & (all_cells_df['TileY'] == cell_y)]
if len(cell_data) == 0:
return []
ref_lat = cell_data['Lat'].iloc[0]
ref_lon = cell_data['Lon'].iloc[0]
# Calcoliamo le distanze per tutte le altre celle
distances = np.sqrt((all_cells_df['Lat'] - ref_lat)**2 + (all_cells_df['Lon'] - ref_lon)**2)
# Filtriamo celle vicine
nearby_mask = distances <= max_distance
nearby_cells = all_cells_df[nearby_mask]
return nearby_cells
def infer_flows_optimized_indexed(snapshots):
"""
Inferisce i flussi con complessità ridotta usando un INDICE SPAZIALE per la ricerca geografica.
"""
print("Inferenza ottimizzata dei flussi con INDICE SPAZIALE...")
flows = []
# Analizziamo ogni fascia oraria
for hour in range(snapshots['hour'].min(), snapshots['hour'].max()):
print(f" Analizzando ora {hour}:00...")
current_snapshot = snapshots[snapshots['hour'] == hour]
next_snapshot = snapshots[snapshots['hour'] == hour + 1]
if len(current_snapshot) == 0 or len(next_snapshot) == 0:
continue
# --- PREPARAZIONE INDICE SPAZIALE (una volta per ora) ---
all_cells_data = pd.concat([current_snapshot, next_snapshot]).drop_duplicates(subset=['TileX', 'TileY'])
gdf = gpd.GeoDataFrame(
all_cells_data,
geometry=gpd.points_from_xy(all_cells_data.Lon, all_cells_data.Lat),
crs="EPSG:4326"
)
# L'indice spaziale (sindex) viene creato automaticamente da GeoPandas
current_presence = dict(zip(zip(current_snapshot['TileX'], current_snapshot['TileY']), current_snapshot['P']))
next_presence = dict(zip(zip(next_snapshot['TileX'], next_snapshot['TileY']), next_snapshot['P']))
all_cells_tuples = set(current_presence.keys()) | set(next_presence.keys())
for cell_tuple in all_cells_tuples:
current_p = current_presence.get(cell_tuple, 0)
next_p = next_presence.get(cell_tuple, 0)
if current_p > next_p and (current_p - next_p) > MIN_FLOW_THRESHOLD:
outflow = current_p - next_p
cell_x, cell_y = cell_tuple
# --- RICERCA OTTIMIZZATA CON INDICE SPAZIALE ---
source_geom = gdf[(gdf.TileX == cell_x) & (gdf.TileY == cell_y)].iloc[0].geometry
search_area = source_geom.buffer(SEARCH_RADIUS_DEGREES)
# 1. Trova gli indici delle celle vicine in modo ultra-veloce
possible_matches_indices = list(gdf.sindex.intersection(search_area.bounds))
nearby_gdf = gdf.iloc[possible_matches_indices]
# 2. Filtra ulteriormente per quelle che sono aumentate di presenza
nearby_increases = []
for _, nearby_row in nearby_gdf.iterrows():
other_cell_tuple = (nearby_row['TileX'], nearby_row['TileY'])
if other_cell_tuple != cell_tuple:
other_current = current_presence.get(other_cell_tuple, 0)
other_next = next_presence.get(other_cell_tuple, 0)
if other_next > other_current:
distance = source_geom.distance(nearby_row.geometry)
increase = other_next - other_current
nearby_increases.append((other_cell_tuple, increase, distance))
# Il resto della logica di distribuzione rimane identico...
if nearby_increases:
total_increase = sum(inc[1] for inc in nearby_increases)
if total_increase > 0:
for dest_cell, dest_increase, distance in nearby_increases:
flow_proportion = (dest_increase / total_increase) * (1 / (1 + distance * DISTANCE_PENALTY_FACTOR))
flow_count = int(outflow * flow_proportion)
if flow_count > 0:
flows.append({
'origin_tile': f"{cell_x}_{cell_y}",
'destination_tile': f"{dest_cell[0]}_{dest_cell[1]}",
'flow_count': flow_count,
'hour': hour,
'distance': distance
})
flows_df = pd.DataFrame(flows)
print(f"Flussi inferiti (con indice spaziale): {len(flows_df)} movimenti")
return flows_df
def create_tessellation_accurate(flows_df, tiles_df):
"""
Crea una tessellazione ACCURATA usando le geometrie WKT dal file delle tile.
"""
print("Creazione tessellazione accurata da dati reali...")
# Estraiamo tutti i tile unici dai flussi
origin_tiles = flows_df['origin_tile'].unique()
destination_tiles = flows_df['destination_tile'].unique()
all_tile_ids = set(np.concatenate([origin_tiles, destination_tiles]))
# Mappiamo i tile_id (es. "140021_94871") ai TileX e TileY numerici
tile_coords = pd.DataFrame([tid.split('_') for tid in all_tile_ids], columns=['TileX', 'TileY'])
tile_coords['TileX'] = tile_coords['TileX'].astype(int)
tile_coords['TileY'] = tile_coords['TileY'].astype(int)
# Uniamo con il dataframe originale delle tile per ottenere le geometrie WKT
# Ci assicuriamo di avere solo le tile effettivamente usate nei flussi
tessellation_data = pd.merge(tile_coords, tiles_df, on=['TileX', 'TileY'], how='inner')
if 'Tile' not in tessellation_data.columns:
raise KeyError("La colonna 'Tile' contenente le geometrie WKT non è stata trovata nel file tiles_d.csv")
# Convertiamo le stringhe WKT in oggetti geometria di Shapely
tessellation_data['geometry'] = tessellation_data['Tile'].apply(wkt.loads)
# Creiamo il GeoDataFrame finale per la tessellazione
# Il tile_ID deve essere una stringa per il matching con FlowDataFrame
tessellation_data['tile_ID'] = tessellation_data['TileX'].astype(str) + '_' + tessellation_data['TileY'].astype(str)
tessellation_gdf = gpd.GeoDataFrame(tessellation_data[['tile_ID', 'geometry']], crs='EPSG:4326')
print(f"Tessellazione accurata creata: {len(tessellation_gdf)} tile")
return tessellation_gdf
def create_flowdataframe_optimized(flows_df, tiles_df):
"""
Crea un FlowDataFrame ottimizzato per analisi con scikit-mobility
"""
print("Creazione FlowDataFrame ottimizzato...")
if len(flows_df) == 0:
print("Nessun flusso inferito. Creazione FlowDataFrame vuoto.")
# Creiamo un FlowDataFrame vuoto ma valido
empty_data = pd.DataFrame({
'origin_tile': ['dummy'],
'destination_tile': ['dummy'],
'flow_count': [0]
})
# Tessellazione minima
tessellation = gpd.GeoDataFrame({
'tile_ID': ['dummy'],
'geometry': [Polygon([(12.20, 44.42), (12.21, 44.42), (12.21, 44.43), (12.20, 44.43), (12.20, 44.42)])]
}, crs='EPSG:4326')
fdf = FlowDataFrame(empty_data,
origin='origin_tile',
destination='destination_tile',
flow='flow_count',
tessellation=tessellation)
return fdf
# Creiamo la tessellazione accurata
tessellation = create_tessellation_accurate(flows_df, tiles_df)
# Prepariamo i dati per FlowDataFrame
flow_data = flows_df.copy()
# Creiamo il FlowDataFrame con tessellazione
fdf = FlowDataFrame(flow_data,
origin='origin_tile',
destination='destination_tile',
flow='flow_count',
tessellation=tessellation)
print(f"FlowDataFrame ottimizzato creato: {len(fdf)} flussi")
print(f"Colonne disponibili: {list(fdf.columns)}")
print(f"Origini uniche: {fdf[fdf.columns[0]].nunique()}")
print(f"Destinazioni uniche: {fdf[fdf.columns[1]].nunique()}")
return fdf
def analyze_flows_optimized(fdf):
"""
Analizza i flussi ottimizzati usando scikit-mobility
"""
print("\n=== ANALISI FLUSSI OTTIMIZZATA ===")
# 1. Analisi dei flussi origine-destinazione
print("Analisi dei flussi origine-destinazione...")
if origin_destination_flows is not None:
try:
od_flows = origin_destination_flows(fdf)
print("Analisi origine-destinazione completata")
except Exception as e:
print(f"Errore nell'analisi origine-destinazione: {e}")
od_flows = None
else:
print("Funzione origine-destinazione non disponibile, implementazione manuale")
if len(fdf) > 0:
od_flows = fdf.groupby(['origin', 'destination'])['flow'].sum().reset_index()
else:
od_flows = None
# 2. Modello gravitazionale
print("Applicazione modello gravitazionale...")
if Gravity is not None and len(fdf) > 0:
try:
gravity_model = Gravity()
gravity_flows = gravity_model.fit(fdf)
print("Modello gravitazionale applicato con successo")
except Exception as e:
print(f"Errore nel modello gravitazionale: {e}")
gravity_flows = None
else:
print("Modello gravitazionale non disponibile o dati insufficienti")
gravity_flows = None
# 3. Modello di radiazione
print("Applicazione modello di radiazione...")
if Radiation is not None and len(fdf) > 0:
try:
radiation_model = Radiation()
radiation_flows = radiation_model.fit(fdf)
print("Modello di radiazione applicato con successo")
except Exception as e:
print(f"Errore nel modello di radiazione: {e}")
radiation_flows = None
else:
print("Modello di radiazione non disponibile o dati insufficienti")
radiation_flows = None
return {
'od_flows': od_flows,
'gravity_flows': gravity_flows,
'radiation_flows': radiation_flows,
'flow_dataframe': fdf
}
def create_flow_visualizations_optimized(flow_analysis, output_dir="analysis_results"):
"""
Crea visualizzazioni ottimizzate per l'analisi dei flussi
"""
import os
os.makedirs(output_dir, exist_ok=True)
print("\n=== CREAZIONE VISUALIZZAZIONI FLUSSI OTTIMIZZATE ===")
# Creiamo una figura grande con multiple subplot
fig = plt.figure(figsize=(20, 15))
# 1. Distribuzione dei flussi
plt.subplot(2, 3, 1)
flow_counts = flow_analysis['flow_dataframe']['flow'].values
if len(flow_counts) > 0 and np.any(flow_counts > 0):
plt.hist(flow_counts, bins=50, alpha=0.7, color='blue', density=True)
plt.title('Distribuzione Intensità Flussi\n(Ottimizzata)', fontsize=12)
plt.xlabel('Intensità del flusso')
plt.ylabel('Densità di probabilità')
plt.grid(True, alpha=0.3)
else:
plt.text(0.5, 0.5, 'Nessun flusso inferito', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Distribuzione Intensità Flussi\n(Nessun dato)', fontsize=12)
# 2. Flussi per fascia oraria
plt.subplot(2, 3, 2)
if 'hour' in flow_analysis['flow_dataframe'].columns and len(flow_analysis['flow_dataframe']) > 0:
hourly_flows = flow_analysis['flow_dataframe'].groupby('hour')['flow'].sum()
if len(hourly_flows) > 0:
plt.bar(hourly_flows.index, hourly_flows.values, alpha=0.7, color='green')
plt.title('Flussi Totali per Fascia Oraria', fontsize=12)
plt.xlabel('Ora del giorno')
plt.ylabel('Flussi totali')
plt.grid(True, alpha=0.3)
else:
plt.text(0.5, 0.5, 'Nessun flusso per ora', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Flussi per Fascia Oraria\n(Nessun dato)', fontsize=12)
else:
plt.text(0.5, 0.5, 'Nessun dato temporale', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Flussi per Fascia Oraria\n(Nessun dato)', fontsize=12)
# 3. Top origini
plt.subplot(2, 3, 3)
if len(flow_analysis['flow_dataframe']) > 0:
top_origins = flow_analysis['flow_dataframe'].groupby('origin')['flow'].sum().nlargest(10)
if len(top_origins) > 0:
plt.barh(range(len(top_origins)), top_origins.values, alpha=0.7, color='red')
plt.yticks(range(len(top_origins)), [f"Tile {o}" for o in top_origins.index])
plt.title('Top 10 Origini per Flusso', fontsize=12)
plt.xlabel('Flusso totale')
plt.grid(True, alpha=0.3)
else:
plt.text(0.5, 0.5, 'Nessun flusso', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Top Origini\n(Nessun dato)', fontsize=12)
else:
plt.text(0.5, 0.5, 'Nessun dato', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Top Origini\n(Nessun dato)', fontsize=12)
# 4. Top destinazioni
plt.subplot(2, 3, 4)
if len(flow_analysis['flow_dataframe']) > 0:
top_destinations = flow_analysis['flow_dataframe'].groupby('destination')['flow'].sum().nlargest(10)
if len(top_destinations) > 0:
plt.barh(range(len(top_destinations)), top_destinations.values, alpha=0.7, color='orange')
plt.yticks(range(len(top_destinations)), [f"Tile {d}" for d in top_destinations.index])
plt.title('Top 10 Destinazioni per Flusso', fontsize=12)
plt.xlabel('Flusso totale')
plt.grid(True, alpha=0.3)
else:
plt.text(0.5, 0.5, 'Nessun flusso', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Top Destinazioni\n(Nessun dato)', fontsize=12)
else:
plt.text(0.5, 0.5, 'Nessun dato', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Top Destinazioni\n(Nessun dato)', fontsize=12)
# 5. Distanza vs flusso
plt.subplot(2, 3, 5)
if 'distance' in flow_analysis['flow_dataframe'].columns and len(flow_analysis['flow_dataframe']) > 0:
distances = flow_analysis['flow_dataframe']['distance'].values
flows = flow_analysis['flow_dataframe']['flow'].values
if len(distances) > 0 and np.any(flows > 0):
plt.scatter(distances, flows, alpha=0.6, s=20, color='purple')
plt.title('Distanza vs Intensità Flusso', fontsize=12)
plt.xlabel('Distanza (gradi)')
plt.ylabel('Intensità del flusso')
plt.grid(True, alpha=0.3)
else:
plt.text(0.5, 0.5, 'Nessun flusso valido', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Distanza vs Flusso\n(Nessun dato)', fontsize=12)
else:
plt.text(0.5, 0.5, 'Nessun dato di distanza', ha='center', va='center', transform=plt.gca().transAxes)
plt.title('Distanza vs Flusso\n(Nessun dato)', fontsize=12)
# 6. Statistiche aggregate
plt.subplot(2, 3, 6)
stats_data = {
'Metrica': ['Flussi totali', 'Origini uniche', 'Destinazioni uniche',
'Flusso medio', 'Flusso max'],
'Valore': [0, 0, 0, 0, 0]
}
if len(flow_analysis['flow_dataframe']) > 0:
stats_data['Valore'][0] = len(flow_analysis['flow_dataframe'])
stats_data['Valore'][1] = flow_analysis['flow_dataframe']['origin'].nunique()
stats_data['Valore'][2] = flow_analysis['flow_dataframe']['destination'].nunique()
stats_data['Valore'][3] = flow_analysis['flow_dataframe']['flow'].mean()
stats_data['Valore'][4] = flow_analysis['flow_dataframe']['flow'].max()
stats_df = pd.DataFrame(stats_data)
plt.bar(stats_df['Metrica'], stats_df['Valore'], color=['blue', 'green', 'red', 'purple', 'orange'])
plt.title('Statistiche Aggregate Flussi Ottimizzati', fontsize=12)
plt.xticks(rotation=45)
plt.ylabel('Valore')
plt.tight_layout()
plt.savefig(f'{output_dir}/flow_analysis_optimized.png', dpi=300, bbox_inches='tight')
plt.show()
print(f"Visualizzazioni flussi ottimizzate salvate in: {output_dir}/flow_analysis_optimized.png")
def generate_flow_report_optimized(flow_analysis, output_file="flow_analysis_optimized_report.txt"):
"""
Genera un report ottimizzato dell'analisi dei flussi
"""
print(f"\n=== GENERAZIONE REPORT FLUSSI OTTIMIZZATO ===")
with open(output_file, 'w', encoding='utf-8') as f:
f.write("REPORT ANALISI FLUSSI OTTIMIZZATA: INFERENZA DA DATI DI PRESENZA (31 MAGGIO 2025)\n")
f.write("=" * 80 + "\n\n")
f.write("METODOLOGIA OTTIMIZZATA:\n")
f.write("- Inferenza dei flussi da snapshot temporali consecutivi\n")
f.write("- Limitazione geografica: ricerca solo in celle vicine (max 0.01°)\n")
f.write("- Filtro di presenza: solo celle con > 5 persone\n")
f.write("- Filtro di flusso: solo movimenti > 10 persone\n")
f.write("- Distribuzione proporzionale dei flussi in uscita\n")
f.write("- Utilizzo di FlowDataFrame di scikit-mobility con tessellazione\n\n")
fdf = flow_analysis['flow_dataframe']
f.write("DATI ANALIZZATI:\n")
f.write(f"- Flussi inferiti: {len(fdf):,}\n")
f.write(f"- Origini uniche: {fdf['origin'].nunique():,}\n")
f.write(f"- Destinazioni uniche: {fdf['destination'].nunique():,}\n")
if len(fdf) > 0:
f.write(f"- Intensità media flusso: {fdf['flow'].mean():.2f}\n")
f.write(f"- Intensità massima flusso: {fdf['flow'].max():.0f}\n")
else:
f.write("- Intensità media flusso: Nessun dato\n")
f.write("- Intensità massima flusso: Nessun dato\n")
f.write("\n")
f.write("OTTIMIZZAZIONI IMPLEMENTATE:\n")
f.write("- Complessità ridotta da O(n²) a O(n × k) dove k << n\n")
f.write("- Limitazione geografica: max 0.01° di distanza\n")
f.write("- Filtri di presenza per ridurre il numero di celle\n")
f.write("- Calcolo efficiente delle distanze\n")
f.write("- Tessellazione automatica per FlowDataFrame\n")
f.write("- Gestione di casi edge e dati insufficienti\n\n")
f.write("ANALISI DEI FLUSSI:\n")
if len(fdf) > 0:
# Top origini
top_origins = fdf.groupby('origin')['flow'].sum().nlargest(5)
f.write("Top 5 Origini:\n")
for i, (origin, flow) in enumerate(top_origins.items(), 1):
f.write(f" {i}. Tile {origin}: {flow:.0f} flussi\n")
f.write("\nTop 5 Destinazioni:\n")
top_destinations = fdf.groupby('destination')['flow'].sum().nlargest(5)
for i, (dest, flow) in enumerate(top_destinations.items(), 1):
f.write(f" {i}. Tile {dest}: {flow:.0f} flussi\n")
else:
f.write("- Nessun flusso significativo inferito\n")
f.write("\nPATTERN TEMPORALI:\n")
if 'hour' in fdf.columns and len(fdf) > 0:
hourly_flows = fdf.groupby('hour')['flow'].sum()
f.write("Flussi per fascia oraria:\n")
for hour, flow in hourly_flows.items():
f.write(f" {hour:02d}:00 - {flow:.0f} flussi\n")
else:
f.write("- Nessun pattern temporale identificato\n")
f.write("\nIMPLICAZIONI PER LA PIANIFICAZIONE:\n")
f.write("- Identificazione di corridoi di mobilità principali\n")
f.write("- Ottimizzazione dei servizi di trasporto\n")
f.write("- Pianificazione di infrastrutture di mobilità\n")
f.write("- Analisi della domanda di trasporto per fascia oraria\n")
f.write("- Strategie di marketing per servizi di mobilità\n\n")
f.write("VALIDAZIONE DELLA METODOLOGIA:\n")
f.write("- Inferenza basata su principi di conservazione della massa\n")
f.write("- Limitazione geografica realistica per mobilità urbana\n")
f.write("- Distribuzione proporzionale ai cambiamenti di presenza\n")
f.write("- Considerazione della distanza per la distribuzione\n")
f.write("- Utilizzo di modelli consolidati (Gravity, Radiation)\n")
f.write("- Tessellazione automatica per analisi spaziale\n\n")
f.write("LIMITAZIONI:\n")
f.write("- L'inferenza non è deterministicamente accurata\n")
f.write("- Non possiamo tracciare individui specifici\n")
f.write("- I flussi sono stime basate su pattern aggregati\n")
f.write("- La distribuzione è basata su euristiche\n")
f.write("- Limitazione geografica può perdere flussi a lunga distanza\n\n")
f.write("RIFERIMENTI BIBLIOGRAFICI:\n")
f.write("- Pappalardo, L., et al. (2015). Returners and Explorers dichotomy in human mobility\n")
f.write("- Song, C., et al. (2010). Modelling the scaling properties of human mobility\n")
f.write("- Jiang, S., et al. (2016). The TimeGeo modeling framework for urban mobility\n")
f.write("- Alessandretti, L., et al. (2018). Evidence for a conserved quantity in human mobility\n")
print(f"Report flussi ottimizzato salvato in: {output_file}")
def main():
"""
Funzione principale per l'analisi ottimizzata dei flussi con scikit-mobility
"""
print("=== ANALISI FLUSSI OTTIMIZZATA CON SCIKIT-MOBILITY ===")
print("Inferenza di flussi da dati di presenza aggregata (versione ottimizzata)\n")
# Carichiamo i dati di mobilità
traffic_data = load_and_prepare_traffic_data()
if traffic_data is None:
print("Errore: Impossibile caricare i dati di mobilità.")
return
# Carichiamo anche i dati delle tile per la tessellazione accurata
tiles_df = pd.read_csv('superset_data/datasets/pixel/tiles_d.csv')
# Creiamo snapshot temporali ottimizzati
snapshots = create_temporal_snapshots_optimized(traffic_data)
# Inferiamo i flussi con ottimizzazione e indice spaziale
flows_df = infer_flows_optimized_indexed(snapshots)
# Creiamo il FlowDataFrame ottimizzato
fdf = create_flowdataframe_optimized(flows_df, tiles_df)
# Analizziamo i flussi con scikit-mobility
flow_analysis = analyze_flows_optimized(fdf)
# Creiamo visualizzazioni
create_flow_visualizations_optimized(flow_analysis)
# Generiamo il report
generate_flow_report_optimized(flow_analysis)
print("\n=== ANALISI FLUSSI OTTIMIZZATA COMPLETATA ===")
print("L'analisi ottimizzata dei flussi con scikit-mobility è stata completata.")
print("I flussi sono stati inferiti dai dati di presenza aggregata con complessità ridotta.")
# Stampiamo alcune statistiche
print(f"\nStatistiche di verifica:")
print(f"- Numero di flussi inferiti: {len(fdf)}")
print(f"- Origini uniche: {fdf['origin'].nunique()}")
print(f"- Destinazioni uniche: {fdf['destination'].nunique()}")
if len(fdf) > 0:
print(f"- Intensità media: {fdf['flow'].mean():.2f}")
print(f"- Intensità massima: {fdf['flow'].max():.0f}")
else:
print("- Nessun flusso significativo inferito")
if __name__ == "__main__":
main() Decomposizione del metodo e dell’algoritmo
Il codice sopra implementa un pipeline di analisi in più fasi. Comprendere ogni passaggio è essenziale per valutare la validità dei risultati.
1. Parametrizzazione del modello, cioè definire le ipotesi
All’inizio del codice, abbiamo definito una serie di parametri chiave. Questa non è solo una buona pratica informatica, ma il cuore della nostra modellazione. Rendono le nostre ipotesi esplicite, trasparenti e modificabili.
MIN_PRESENCE_THRESHOLDeMIN_FLOW_THRESHOLD: Sono filtri statistici. Ci permettono di escludere il “rumore” di fondo (celle quasi vuote o micro-spostamenti insignificanti) per concentrarci sugli eventi di mobilità di massa che hanno maggiore probabilità di essere segnali reali.SEARCH_RADIUS_DEGREES: È il parametro più importante, poiché definisce la scala spaziale e, di conseguenza, il tipo di mobilità che stiamo studiando. Un raggio di ricerca limitato è l’ipotesi chiave che afferma “le persone si muovono nelle vicinanze”. Ma cosa significa “vicino”? Dipende da come ci si muove.- Per la mobilità pedonale: Una camminata di 15 minuti copre circa 1.25 km. Il parametro andrebbe impostato a
SEARCH_RADIUS_DEGREES = 0.011. - Per la mobilità veicolare: Un tragitto in auto di 15 minuti in città può coprire 5 km o più. Il parametro andrebbe impostato a
SEARCH_RADIUS_DEGREES = 0.045.
- Per la mobilità pedonale: Una camminata di 15 minuti copre circa 1.25 km. Il parametro andrebbe impostato a
2. L’Inferenza dei flussi con indice spaziale
La funzione infer_flows_optimized_indexed è il motore del modello. La sua logica si basa sui principi discussi nell’introduzione, ma con un’ottimizzazione cruciale: l’indice spaziale. Anziché confrontare una cella di partenza con tutte le altre decine di migliaia di celle (approccio O(N²)), per ogni ora creiamo un GeoDataFrame che organizza le celle in una struttura dati geografica efficiente (un R-tree). Quando dobbiamo trovare le destinazioni vicine, interroghiamo questo indice, che restituisce in modo quasi istantaneo solo le celle candidate all’interno del nostro raggio di ricerca. Questo riduce la complessità a O(N*k), dove k (le celle vicine) è un numero molto piccolo rispetto a N (le celle totali).
3. La tessellazione: da pixel a territorio
Perché tassellare? Un “pixel” con coordinate (TileX, TileY) è un’entità astratta. Per un’analisi di mobilità, dobbiamo sapere dove si trova nel mondo reale e che forma ha. La funzione create_tessellation_accurate svolge questo compito critico. Legge la colonna Tile dal file CSV, che contiene la definizione geometrica di ogni cella nel formato standard WKT (Well-Known Text). Usando la libreria Shapely, converte queste stringhe in oggetti Polygon geograficamente consapevoli. Questo processo:
- Ci dà l’accuratezza geografica dei calcoli di distanza e area.
- Abilita analisi spaziali avanzate, come l’intersezione dei flussi con mappe di zone urbane (aree industriali, residenziali, parchi).
- È il presupposto per creare visualizzazioni cartografiche corrette.
4. Creazione di FlowDataFrame
A questo punto, i flussi inferiti e la tessellazione vengono uniti in un FlowDataFrame, l’oggetto specializzato di scikit-mobility. Questo non è solo un contenitore di dati: è un oggetto “intelligente” che espone metodi per l’analisi avanzata della mobilità, come l’applicazione di modelli standard (Gravity, Radiation) con cui confrontare i nostri risultati.
Limiti/prospettive
Questo pipeline trasforma dati di presenza, chiaramente grezzi, in una rete di flussi georeferenziata, cioè uno strumento per comprendere alcune dinamiche urbane. È però importante essere consapevoli dei limiti di questo modello:
- È un’inferenza, non una certezza: Il modello stima i flussi più probabili, ma non può tracciare i singoli individui (un limite noto come fallacia ecologica).
- Le ipotesi sono semplificazioni: Il nostro modello di distribuzione del flusso è una delle tante euristiche possibili.
La vera forza di questo pipeline non risiede tanto nei risultati che produce, ma nel metodo flessibile. E infatti sono aperti sviluppi come:
- Confronto con modelli standard: Il primo passo successivo è confrontare i flussi che abbiamo inferito con quelli generati dai modelli accademici standard come il Modello gravitazionale e il Modello a radiazione (già disponibili in
scikit-mobility). Questo ci permetterebbe di quantificare la “bontà” del nostro modello di inferenza. - Modelli a raggio dinamico: Il nostro modello attuale usa un raggio di ricerca fisso per tutta la città. Un’evoluzione molto potente sarebbe renderlo dinamico: il raggio potrebbe variare a seconda della cella di partenza. Ad esempio, potrebbe essere più grande per una cella situata su una superstrada e molto più piccolo per una cella in un centro storico pedonale. Questo richiederebbe di arricchire i dati delle celle con informazioni sulla tipologia di zona (residenziale, commerciale, stradale, etc.).
- Modelli Multi-Modali: Un approccio alternativo, ma altrettanto potente, è quello multi-modale. Anziché un unico modello complesso, potremmo eseguire il nostro pipeline più volte con impostazioni diverse per isolare i differenti modi di trasporto:
- Run 1 (Pedonale): Con
SEARCH_RADIUS_DEGREES = 0.011per generare una mappa dei flussi a piedi. - Run 2 (Veicolare): Con
SEARCH_RADIUS_DEGREES = 0.045per generare una mappa dei flussi in auto.
- Run 1 (Pedonale): Con
Resta inoltre essenziale l’integrazione con altre fonti di dati (es. orari del trasporto pubblico, localizzazione di scuole e uffici) per validare e arricchire l’interpretazione dei pattern di mobilità emersi.