Importazione delle librerie¶

In [1]:
import os
import folium
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import geopandas as gpd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
In [2]:
ilab_path = os.path.join(os.path.expanduser('~'), 'ILAB_DATA')
work_path = os.path.join(ilab_path, 'PATRIMONIO', 'DATA')
out_path = os.path.join(ilab_path, 'PATRIMONIO', 'OUT')

Caricamento dei dati¶

In [3]:
# Dataset principale
grid_data = pd.read_csv(os.path.join(out_path,'grid_08_adv.csv'))
# Metadati
metadata = pd.read_csv(os.path.join(out_path,'metadati_dataset_08_cat.csv'))

# Crea un dizionario per accesso rapido ai metadati
metadata_dict = {row['nome_colonna']: 
                {'category': row['category'], 'description': row['descrizione']} 
                 for _, row in metadata.iterrows()}

print(f"Dataset originale: {grid_data.shape[0]} righe, {grid_data.shape[1]} colonne")

Dataset originale: 2286 righe, 336 colonne

Preparazione dei dati¶

In [4]:
# Seleziona solo colonne numeriche, escluse quelle non analitiche
numeric_cols = grid_data.select_dtypes(include=[np.number]).columns
selected_cols = [col for col in numeric_cols if 
                col != 'cell_a_m2' and 
                'geometry' not in col and 
                not col.startswith('h3_')]

# Rimuovi colonne con più del 30% di valori mancanti
threshold = len(grid_data) * 0.3
cols_to_keep = []

for col in selected_cols:
    missing_count = grid_data[col].isna().sum()
    if missing_count <= threshold:
        cols_to_keep.append(col)

# Seleziona dataset con colonne rilevanti
data_subset = grid_data[['h3_08'] + cols_to_keep]

# Rimuovi righe con valori mancanti
data_clean = data_subset.dropna()

# Separa identificatori e dati per l'analisi
cell_ids = data_clean['h3_08']
features = data_clean.drop('h3_08', axis=1)

print(f"Dataset pulito: {features.shape[0]} righe, {features.shape[1]} colonne")

Dataset pulito: 1565 righe, 323 colonne

Standardizzazione dei dati¶

In [5]:
# Standardizza le variabili (media=0, dev. std=1)
scaler = StandardScaler()
scaled_data = scaler.fit_transform(features)
feature_names = features.columns.tolist()

print("Dati standardizzati completati.")

Dati standardizzati completati.

Rimozione delle variabili correlate¶

In [6]:
# Identifica e rimuovi variabili altamente correlate
threshold = 0.85

# Calcola matrice di correlazione
corr_matrix = np.corrcoef(scaled_data, rowvar=False)

# Identifica coppie di variabili correlate
correlated_pairs = []
for i in range(len(feature_names)):
    for j in range(i+1, len(feature_names)):
        if abs(corr_matrix[i, j]) > threshold:
            correlated_pairs.append((feature_names[i], feature_names[j], corr_matrix[i, j]))

# Crea un grafo di correlazioni
correlation_graph = {}
for col1, col2, _ in correlated_pairs:
    if col1 not in correlation_graph:
        correlation_graph[col1] = []
    if col2 not in correlation_graph:
        correlation_graph[col2] = []
    
    correlation_graph[col1].append(col2)
    correlation_graph[col2].append(col1)

# Trova gruppi connessi (DFS)
visited = set()
groups = []

# Funzione DFS definita inline
def dfs(node, group):
    visited.add(node)
    group.append(node)
    
    for neighbor in correlation_graph.get(node, []):
        if neighbor not in visited:
            dfs(neighbor, group)

# Trova tutti i gruppi correlati
for col in correlation_graph:
    if col not in visited:
        group = []
        dfs(col, group)
        if group:
            groups.append(group)

# Seleziona un rappresentante per ogni gruppo (il primo elemento)
representatives = [group[0] for group in groups]

# Colonne da mantenere: non correlate + rappresentanti
all_correlated = set([col for group in groups for col in group])
reduced_features = [col for col in feature_names if 
                  col not in all_correlated or col in representatives]

# Crea matrice ridotta
indices_to_keep = [feature_names.index(col) for col in reduced_features]
reduced_data = scaled_data[:, indices_to_keep]

print(f"Gruppi di variabili correlate trovati: {len(groups)}")
print(f"Dataset dopo rimozione correlate: {reduced_data.shape[1]} colonne")

Gruppi di variabili correlate trovati: 5
Dataset dopo rimozione correlate: 115 colonne

Analisi PCA¶

In [7]:
# Esegui l'analisi delle componenti principali
pca_model = PCA()
pc_scores = pca_model.fit_transform(reduced_data)

# Calcola varianza spiegata e varianza cumulativa
explained_variance = pca_model.explained_variance_ratio_
cumulative_variance = np.cumsum(explained_variance)

# Determina il numero di componenti per spiegare l'80% della varianza
n_components_80 = np.argmax(cumulative_variance >= 0.8) + 1

print(f"Componenti necessarie per 80% varianza: {n_components_80}")

# Mostra varianza spiegata
print("\nVarianza spiegata:")
for i in range(min(10, len(explained_variance))):
    print(f"- PC{i+1}: {explained_variance[i]*100:.2f}% (cumulativa: {cumulative_variance[i]*100:.2f}%)")

Componenti necessarie per 80% varianza: 38

Varianza spiegata:
- PC1: 23.17% (cumulativa: 23.17%)
- PC2: 6.71% (cumulativa: 29.89%)
- PC3: 4.54% (cumulativa: 34.43%)
- PC4: 3.37% (cumulativa: 37.80%)
- PC5: 2.98% (cumulativa: 40.77%)
- PC6: 2.59% (cumulativa: 43.36%)
- PC7: 2.13% (cumulativa: 45.49%)
- PC8: 1.99% (cumulativa: 47.48%)
- PC9: 1.93% (cumulativa: 49.41%)
- PC10: 1.77% (cumulativa: 51.18%)

Interpretazione delle componenti principali¶

In [8]:
# Interpreta le componenti principali
n_components = 5  # Numero di componenti da interpretare
n_top_features = 5  # Numero di variabili più importanti per componente

interpretations = []

# Per ogni componente, trova le variabili con maggior peso
for i in range(min(n_components, pca_model.components_.shape[0])):
    # Calcola i loadings (pesi assoluti)
    loadings = [(reduced_features[j], pca_model.components_[i, j], 
                metadata_dict.get(reduced_features[j], {}).get('category', 'N/A')) 
                for j in range(len(reduced_features))]
    
    # Ordina per peso assoluto decrescente
    sorted_loadings = sorted(loadings, key=lambda x: abs(x[1]), reverse=True)
    
    # Seleziona le n_top_features più importanti
    top_features = [(name, loading, '+' if loading >= 0 else '-', category) 
                    for name, loading, category in sorted_loadings[:n_top_features]]
    
    interpretations.append({
        'component': i + 1,
        'variance_explained': pca_model.explained_variance_ratio_[i],
        'top_features': top_features
    })

print("\nInterpretazione delle componenti principali:")
for interp in interpretations:
    print(f"\nPC{interp['component']} ({interp['variance_explained']*100:.2f}% varianza):")
    for feature, loading, sign, category in interp['top_features']:
        print(f"  {sign} {feature} ({category}): {abs(loading):.4f}")

Interpretazione delle componenti principali:

PC1 (23.17% varianza):
  + sanita_ospedale_transit||transit__60__ospedale (Accessibilità sanità - Ospedali): 0.1696
  + ovm_infr_infrastrutture||walking__30__stazione (N/A): 0.1556
  + sanita_ospedale_transit||transit__60__AZIENDA OSPEDALIERA (Accessibilità sanità - Ospedali): 0.1523
  + ovm_places_sport||walking__15__formazione_sportiva (Accessibilità sport): 0.1484
  + ovm_places_cultura||walking__30__evento_culturale (N/A): 0.1451

PC2 (6.71% varianza):
  - sanita_ospedale||traffic__30__IRCCS PRIVATO (N/A): 0.2399
  - sanita_ospedale||traffic__30__IRCCS PUBBLICO (N/A): 0.2232
  - sanita_ospedale||traffic__15__OSPED. CLASSIF. O ASSIMIL. ART 1 L.132/1968 (Accessibilità sanità - Ospedali): 0.2080
  - sanita_ospedale||traffic__30__IRCCS FONDAZIONE (N/A): 0.1997
  - sanita_ospedale||traffic__30__POLICLINICO UNIVERSITARIO PRIVATO (N/A): 0.1994

PC3 (4.54% varianza):
  - ovm_places_sport||traffic__30__luoghi_sportivi (N/A): 0.3147
  - sanita_ospedale||traffic__30__A.O. INTEGRATA CON L'UNIVERSITA' (N/A): 0.3041
  - sanita_ospedale_transit||transit__60__A.O. INTEGRATA CON L'UNIVERSITA' (Accessibilità sanità - Ospedali): 0.2739
  - ovm_places_sport||traffic__15__luoghi_sportivi (Accessibilità sport): 0.2677
  - sanita_ospedale||traffic__15__A.O. INTEGRATA CON L'UNIVERSITA' (Accessibilità sanità - Ospedali): 0.2182

PC4 (3.37% varianza):
  + ovm_places_istruzione||walking__30__formazione_adulti (N/A): 0.2776
  - sanita_ospedale||walking__30__OSPED. A GESTIONE DIRETTA PRESIDIO A.S.L. (N/A): 0.2058
  + ovm_places_istruzione_transit||transit__15__formazione_adulti (Accessibilità istruzione): 0.1984
  + ovm_places_istruzione_transit||transit__30__formazione_adulti (Accessibilità istruzione): 0.1919
  + ovm_places_istruzione||walking__15__formazione_adulti (Accessibilità istruzione): 0.1870

PC5 (2.98% varianza):
  + ovm_places_sanita||walking__30__cliniche_mediche (N/A): 0.2393
  - ovm_places_istruzione_transit||transit__30__formazione_alternativa (Accessibilità istruzione): 0.2270
  - sanita_ospedale||traffic__15__AZIENDA OSPEDALIERA (Accessibilità sanità - Ospedali): 0.2104
  + sanita_ospedale||walking__30__A.O. INTEGRATA CON IL SSN (N/A): 0.2085
  + ovm_places_sanita||walking__15__cliniche_mediche (Accessibilità sanità): 0.2052

Preparazione per il clustering¶

In [9]:
# Seleziona le prime n_components_80 componenti per il clustering
pc_data_for_clustering = pc_scores[:, :n_components_80]

Clustering K-means¶

In [10]:
# Esegui il clustering K-means e valuta il numero ottimale di cluster
max_k = 8

# Valuta diversi valori di k
silhouette_scores = []
inertia_values = []

for k in range(2, max_k + 1):
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    clusters = kmeans.fit_predict(pc_data_for_clustering)
    
    # Calcola silhouette score
    sil_score = silhouette_score(pc_data_for_clustering, clusters)
    silhouette_scores.append(sil_score)
    inertia_values.append(kmeans.inertia_)
    
    print(f"k={k}, Silhouette Score={sil_score:.4f}, Inertia={kmeans.inertia_:.2f}")

# Trova k ottimale (massimo silhouette score)
best_k = np.argmax(silhouette_scores) + 2  # +2 perché partiamo da k=2

# Esegui clustering finale con k ottimale
final_kmeans = KMeans(n_clusters=best_k, random_state=42, n_init=10)
final_clusters = final_kmeans.fit_predict(pc_data_for_clustering)

print(f"Numero ottimale di cluster: {best_k}")

k=2, Silhouette Score=0.5368, Inertia=115304.95

k=3, Silhouette Score=0.2980, Inertia=104854.58

k=4, Silhouette Score=0.2449, Inertia=97843.14

k=5, Silhouette Score=0.2632, Inertia=92741.36

k=6, Silhouette Score=0.2660, Inertia=88398.05

k=7, Silhouette Score=0.2486, Inertia=85244.38

k=8, Silhouette Score=0.2641, Inertia=82036.88

Numero ottimale di cluster: 2

Analisi dei cluster¶

In [11]:
# Analizza le caratteristiche di ogni cluster
cluster_profiles = []

# Per ogni cluster, calcola media delle componenti principali
for i in range(best_k):
    cluster_points = pc_data_for_clustering[final_clusters == i]
    means = np.mean(cluster_points, axis=0)
    cluster_profiles.append(means)

# Interpreta i cluster
cluster_interpretations = []

for i, profile in enumerate(cluster_profiles):
    interpretation = {
        'cluster': i + 1,
        'size': np.sum(final_clusters == i),
        'percentage': np.sum(final_clusters == i) / len(final_clusters) * 100,
        'significant_components': []
    }
    
    # Trova componenti significative (valore assoluto > 0.5)
    for j, value in enumerate(profile):
        if abs(value) > 0.5:
            direction = 'alto' if value > 0 else 'basso'
            component_interp = interpretations[j] if j < len(interpretations) else None
            
            interpretation['significant_components'].append({
                'component': j + 1,
                'value': value,
                'direction': direction,
                'interpretation': component_interp
            })
    
    cluster_interpretations.append(interpretation)

print("\nAnalisi dei cluster:")
for interp in cluster_interpretations:
    print(f"\nCluster {interp['cluster']}: {interp['size']} celle ({interp['percentage']:.2f}%)")
    
    if interp['significant_components']:
        for comp in interp['significant_components']:
            print(f"  PC{comp['component']} ({comp['direction']}): {comp['value']:.2f}")
            
            if comp['interpretation']:
                for feature, loading, sign, category in comp['interpretation']['top_features'][:3]:
                    effective_sign = sign if comp['direction'] == 'alto' else ('+' if sign == '-' else '-')
                    print(f"    {effective_sign} {feature} ({category})")
    else:
        print("  Cluster con caratteristiche nella media")

Analisi dei cluster:

Cluster 1: 175 celle (11.18%)
  PC1 (alto): 12.12
    + sanita_ospedale_transit||transit__60__ospedale (Accessibilità sanità - Ospedali)
    + ovm_infr_infrastrutture||walking__30__stazione (N/A)
    + sanita_ospedale_transit||transit__60__AZIENDA OSPEDALIERA (Accessibilità sanità - Ospedali)
  PC3 (alto): 0.54
    - ovm_places_sport||traffic__30__luoghi_sportivi (N/A)
    - sanita_ospedale||traffic__30__A.O. INTEGRATA CON L'UNIVERSITA' (N/A)
    - sanita_ospedale_transit||transit__60__A.O. INTEGRATA CON L'UNIVERSITA' (Accessibilità sanità - Ospedali)
  PC4 (basso): -1.01
    - ovm_places_istruzione||walking__30__formazione_adulti (N/A)
    + sanita_ospedale||walking__30__OSPED. A GESTIONE DIRETTA PRESIDIO A.S.L. (N/A)
    - ovm_places_istruzione_transit||transit__15__formazione_adulti (Accessibilità istruzione)

Cluster 2: 1390 celle (88.82%)
  PC1 (basso): -1.53
    - sanita_ospedale_transit||transit__60__ospedale (Accessibilità sanità - Ospedali)
    - ovm_infr_infrastrutture||walking__30__stazione (N/A)
    - sanita_ospedale_transit||transit__60__AZIENDA OSPEDALIERA (Accessibilità sanità - Ospedali)

Creazione del dataset finale¶

In [12]:
# Crea dataset finale con risultati
result_df = pd.DataFrame({
    'h3_08': cell_ids.values,
    'cluster': final_clusters + 1  # Cluster da 1 a K
})

# Aggiungi le prime componenti principali
for i in range(min(5, pc_scores.shape[1])):
    result_df[f'PC{i+1}'] = pc_scores[:, i]

print("\nEsempio di risultati finali:")
print(result_df.head())

# Salva risultati
result_df.to_csv('cluster_results.csv', index=False)

print("\nAnalisi completata. Risultati salvati in 'cluster_results.csv'")

Esempio di risultati finali:
             h3_08  cluster       PC1       PC2       PC3       PC4       PC5
0  881e805327fffff        2  0.396462 -4.952539  1.144988  1.800543  2.700710
1  881e805ae7fffff        2 -0.477435 -0.058688  1.099270  0.391273  0.609472
2  881e805803fffff        2  1.177851  2.222827 -0.481049 -2.274696 -0.240375
3  881e80c9dbfffff        2 -2.555567 -0.254078  0.612190  0.232882  0.698368
4  881e804301fffff        2 -1.068426 -2.186474 -0.253502 -0.521432  1.571695

Analisi completata. Risultati salvati in 'cluster_results.csv'

Visualizzazione dei risultati¶

In [13]:
# Visualizza i risultati dell'analisi PCA e del clustering
# Crea una figura con due subplot
fig, axes = plt.subplots(1, 2, figsize=(15, 6))

# 1. Visualizza i primi due componenti principali con colori per cluster
scatter = axes[0].scatter(pc_scores[:, 0], pc_scores[:, 1], c=final_clusters, cmap='viridis', alpha=0.7)
axes[0].set_xlabel(f'PC1 ({explained_variance[0]*100:.2f}%)')
axes[0].set_ylabel(f'PC2 ({explained_variance[1]*100:.2f}%)')
axes[0].set_title('Proiezione PCA con Cluster')
axes[0].grid(True, linestyle='--', alpha=0.7)

# Aggiungi legenda
legend = axes[0].legend(*scatter.legend_elements(), title="Cluster")
axes[0].add_artist(legend)

# 2. Visualizza l'andamento dell'indice di silhouette per diversi k
k_values = range(2, 2 + len(silhouette_scores))
axes[1].plot(k_values, silhouette_scores, 'o-', color='blue')
axes[1].set_xlabel('Numero di cluster (k)')
axes[1].set_ylabel('Indice Silhouette')
axes[1].set_title('Determinazione del numero ottimale di cluster')
axes[1].grid(True, linestyle='--', alpha=0.7)

# Evidenzia il k con massimo silhouette
best_k_idx = np.argmax(silhouette_scores)
best_k = k_values[best_k_idx]
axes[1].axvline(x=best_k, color='red', linestyle='--')
axes[1].text(best_k + 0.1, silhouette_scores[best_k_idx] - 0.05, 
            f'k ottimale = {best_k}', color='red', fontweight='bold')

plt.tight_layout()
plt.savefig('pca_cluster_results.png', dpi=300, bbox_inches='tight')
plt.show()
No description has been provided for this image

Visualizzazione dei risultati su mappa¶

In [14]:
grid = gpd.read_parquet(os.path.join(work_path,'grid_roma_08_32632.geoparquet'))
cluster = pd.read_csv('cluster_results.csv')
grid_clust = grid.merge(cluster, how='left')
In [15]:
m = None
components = ['PC1','PC2','PC3','PC4','PC5']
components_lbl = [
    'Accessibilità ai servizi urbani principali',
    'Presenza di strutture sanitarie specializzate',
    'Accessibilità a impianti sportivi e ospedali universitari',
    'Formazione adulti e presidio sanitario territoriale',
    'Servizi sanitari di prossimità vs formazione alternativa'
    ]
for comp,lbl in zip(components,components_lbl):
    m = grid_clust.explore(m=m, column=comp, name=lbl, show=False)
folium.LayerControl().add_to(m)
m
Out[15]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [ ]: