#! pip install kneed

Tutorial per il clustering geografico

Teoria

import os
import geopandas as gpd
import numpy as np
from sklearn.cluster import KMeans, DBSCAN, HDBSCAN
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import silhouette_score, davies_bouldin_score
from kneed import KneeLocator
from scipy.spatial.distance import cdist
import folium
import plotly.graph_objects as go
import matplotlib.pyplot as plt

ilab_path = os.path.join(os.path.expanduser('~'), 'ILAB_DATA')
work_path = os.path.join(ilab_path, 'SARDEGNA')
cu_path = os.path.join(work_path, 'DBGT_CU')

def visualize_clusters(coordinates, labels, centers=None):
    """
    Crea una mappa interattiva con i cluster evidenziati
    """
    # Calcola il centro medio delle coordinate
    center_lat = np.mean(coordinates[:, 0])
    center_lon = np.mean(coordinates[:, 1])
    
    # Crea la mappa
    m = folium.Map(location=[center_lat, center_lon], zoom_start=13)
    
    # Colori per i cluster
    colors = ['red', 'blue', 'green', 'purple', 'orange', 'darkred',
                'lightred', 'beige', 'darkblue', 'darkgreen', 'cadetblue',
                'darkpurple', 'pink', 'lightblue', 'lightgreen']
    
    # Aggiungi punti alla mappa
    for idx, (lat, lon) in enumerate(coordinates):
        if labels[idx] == -1:
            color = 'gray'  # Punti di rumore
        else:
            color = colors[labels[idx] % len(colors)]
        
        folium.CircleMarker(
            location=[lat, lon],
            radius=2,
            color=color,
            fill=True,
            popup=f'Cluster {labels[idx]}'
        ).add_to(m)
    
    # Aggiungi i centri dei cluster se disponibili
    if centers is not None:
        for idx, center in enumerate(centers):
            folium.CircleMarker(
                location=[center[0], center[1]],
                radius=5,
                color='black',
                fill=True,
                popup=f'Centro Cluster {idx}'
            ).add_to(m)
    
    return m

Load data

Selezioniamo solo le attrezzature sportive di Cagliari

# comuni = gpd.read_parquet(os.path.join(ilab_path, 'ISTAT', 'DATA', 'comuni_pop.parquet'))
# cagliari = comuni[comuni['comune'] == 'Cagliari']
# poly = gpd.read_parquet(os.path.join(cu_path, 'DBGT_CU_Polygon.parquet'))
# poly_cagliari = poly.sjoin(cagliari.to_crs(4326), how='inner', predicate='intersects')

# attrezzatura_sportiva = poly_cagliari[poly_cagliari.Classe == 'ATTREZZATURA_SPORTIVA'].copy()
# attrezzatura_sportiva.to_parquet(os.path.join(cu_path,'attrezzatura_sportiva.parquet'))

attrezzatura_sportiva = gpd.read_parquet(os.path.join(cu_path, 'attrezzatura_sportiva.parquet'))

Calcolo i centroidi dei poligoni

attrezzatura_sportiva['geometry_c'] = attrezzatura_sportiva['geometry'].to_crs(32632).centroid.to_crs(4326)
attrezzatura_sportiva['x'] = attrezzatura_sportiva['geometry_c'].x
attrezzatura_sportiva['y'] = attrezzatura_sportiva['geometry_c'].y

X = attrezzatura_sportiva[['y','x']].values

results = {}

k-means

L'algoritmo k-means richiede il parametro k per funzionare. Selezioniamo dunque il miglior numero di cluster grazie a varie metriche

max_k = 800
random_state = 42

distortions = []
inertia = []
silhouette_scores = []
davies_bouldin_scores = []
K = range(2, max_k + 1)
for k in K:
    kmeans = KMeans(n_clusters=k, random_state=random_state)
    kmeans_labels = kmeans.fit_predict(X)
    distortions.append(sum(np.min(cdist(X, kmeans.cluster_centers_, 'euclidean'), axis=1)) / X.shape[0])
    inertia.append(kmeans.inertia_)
    silhouette_scores.append(silhouette_score(X, kmeans_labels))
    davies_bouldin_scores.append(davies_bouldin_score(X, kmeans_labels))

Plot dei risultati

plt.figure(figsize=(16, 10))

# Metodo del gomito - distortions
plt.subplot(2, 2, 1)
plt.plot(K, distortions, marker='o')
plt.title('Metodo del gomito - distortions')
plt.xlabel('Numero di cluster (k)')
plt.ylabel('Inertia')

# Metodo del gomito - inertia
plt.subplot(2, 2, 2)
plt.plot(K, inertia, marker='o', color='red')
plt.title('Metodo del gomito - inertia')
plt.xlabel('Numero di cluster (k)')
plt.ylabel('Inertia')

# Silhouette Score
plt.subplot(2, 2, 3)
plt.plot(K, silhouette_scores, marker='o', color='orange')
plt.title('Silhouette Score')
plt.xlabel('Numero di cluster (k)')
plt.ylabel('Silhouette Score')

# Davies-Bouldin Index
plt.subplot(2, 2, 4)
plt.plot(K, davies_bouldin_scores, marker='o', color='green')
plt.title('Davies-Bouldin Index')
plt.xlabel('Numero di cluster (k)')
plt.ylabel('Davies-Bouldin Index')

plt.tight_layout()
plt.show()

Selezioniamo il modello migliore

# Trova il punto di gomito
best_k_distorsions = KneeLocator(K, distortions, S=1.0, curve='convex', direction='decreasing').elbow
best_k_inertia = KneeLocator(K, inertia, S=1.0, curve='convex', direction='decreasing').elbow

print(f"Miglior k basato su distortions: {best_k_distorsions}")
print(f"Miglior k basato su inertia: {best_k_inertia}")

# Identifica il miglior k dai risultati
best_k_silhouette = K[np.argmax(silhouette_scores)]
best_k_davies_bouldin = K[np.argmin(davies_bouldin_scores)]

print(f"Miglior k basato sul Silhouette Score: {best_k_silhouette}")
print(f"Miglior k basato sul Davies-Bouldin Index: {best_k_davies_bouldin}")

Miglior k basato su distortions: 67
Miglior k basato su inertia: 34
Miglior k basato sul Silhouette Score: 177
Miglior k basato sul Davies-Bouldin Index: 800

Decidiamo che il migior modello sia quello che produce il maggior numero di cluster (escludendo 800) --> ovvero il miglior silhouette score

kmeans_optimal = KMeans(n_clusters=best_k_silhouette, random_state=random_state)
kmeans_labels = kmeans_optimal.fit_predict(X)
silhouette_avg = silhouette_score(X, kmeans_labels)

results['kmeans'] = {
    'labels': kmeans_labels,
    'n_clusters': best_k_silhouette,
    'centers': kmeans_optimal.cluster_centers_,
    'silhouette_score': silhouette_avg
}

visualize_clusters(X, results['kmeans']['labels'])#, results['kmeans']['centers'])

Make this Notebook Trusted to load map: File -> Trust Notebook

DB-scan

Svolgiamo una grid search manuale sui parametri principali, selezionando come modello migliore quello che produce un maggior numero di cluster

max_clusters = 0
best_configuration = None

metrics = ['euclidean', 'manhattan', 'haversine']
min_samples_values = [3, 5, 10]
eps_values = [0.00001, 0.0001, 0.001, 0.01]

for metric in metrics:
    for min_samples in min_samples_values:
        for eps in eps_values:
            dbscan = DBSCAN(
                metric=metric,
                min_samples=min_samples, 
                eps=eps, 
            )
            clusters = dbscan.fit_predict(X)
            n_clusters = len(set(clusters)) - (1 if -1 in clusters else 0)
            print(f"metric={metric}, min_samples={min_samples}, eps={eps}, n_clusters={n_clusters}")

            # Aggiorna la configurazione con il numero massimo di cluster
            if n_clusters > max_clusters:
                max_clusters = n_clusters
                best_configuration = {
                    'metric': metric,
                    'min_samples': min_samples,
                    'eps': eps,
                    'n_clusters': n_clusters
                }

metric=euclidean, min_samples=3, eps=1e-05, n_clusters=0
metric=euclidean, min_samples=3, eps=0.0001, n_clusters=16
metric=euclidean, min_samples=3, eps=0.001, n_clusters=87
metric=euclidean, min_samples=3, eps=0.01, n_clusters=1
metric=euclidean, min_samples=5, eps=1e-05, n_clusters=0
metric=euclidean, min_samples=5, eps=0.0001, n_clusters=0
metric=euclidean, min_samples=5, eps=0.001, n_clusters=51
metric=euclidean, min_samples=5, eps=0.01, n_clusters=1
metric=euclidean, min_samples=10, eps=1e-05, n_clusters=0
metric=euclidean, min_samples=10, eps=0.0001, n_clusters=0
metric=euclidean, min_samples=10, eps=0.001, n_clusters=10
metric=euclidean, min_samples=10, eps=0.01, n_clusters=1
metric=manhattan, min_samples=3, eps=1e-05, n_clusters=0
metric=manhattan, min_samples=3, eps=0.0001, n_clusters=11
metric=manhattan, min_samples=3, eps=0.001, n_clusters=87
metric=manhattan, min_samples=3, eps=0.01, n_clusters=2
metric=manhattan, min_samples=5, eps=1e-05, n_clusters=0
metric=manhattan, min_samples=5, eps=0.0001, n_clusters=0
metric=manhattan, min_samples=5, eps=0.001, n_clusters=42
metric=manhattan, min_samples=5, eps=0.01, n_clusters=1
metric=manhattan, min_samples=10, eps=1e-05, n_clusters=0
metric=manhattan, min_samples=10, eps=0.0001, n_clusters=0
metric=manhattan, min_samples=10, eps=0.001, n_clusters=9
metric=manhattan, min_samples=10, eps=0.01, n_clusters=2
metric=haversine, min_samples=3, eps=1e-05, n_clusters=4
metric=haversine, min_samples=3, eps=0.0001, n_clusters=69

metric=haversine, min_samples=3, eps=0.001, n_clusters=7
metric=haversine, min_samples=3, eps=0.01, n_clusters=1
metric=haversine, min_samples=5, eps=1e-05, n_clusters=0
metric=haversine, min_samples=5, eps=0.0001, n_clusters=17
metric=haversine, min_samples=5, eps=0.001, n_clusters=8
metric=haversine, min_samples=5, eps=0.01, n_clusters=1
metric=haversine, min_samples=10, eps=1e-05, n_clusters=0
metric=haversine, min_samples=10, eps=0.0001, n_clusters=0
metric=haversine, min_samples=10, eps=0.001, n_clusters=3

metric=haversine, min_samples=10, eps=0.01, n_clusters=1

print(f"Miglior configurazione: {best_configuration}")

Miglior configurazione: {'metric': 'euclidean', 'min_samples': 3, 'eps': 0.001, 'n_clusters': 87}

dbscan = DBSCAN(
    metric=best_configuration['metric'],
    min_samples=best_configuration['min_samples'], 
    eps=best_configuration['eps'],
)
dbscan_labels = dbscan.fit_predict(X)
results['dbscan'] = {
    'labels': dbscan_labels,
    'n_clusters': len(set(dbscan_labels)) - (1 if -1 in dbscan_labels else 0)
}

visualize_clusters(X, results['dbscan']['labels'])

Make this Notebook Trusted to load map: File -> Trust Notebook

HDB-scan

Svolgiamo una grid search manuale sui parametri principali, selezionando come modello migliore quello che produce un maggior numero di cluster

max_clusters = 0
best_configuration = None

metrics = ['euclidean', 'manhattan', 'haversine']
min_cluster_sizes = [3, 5, 10]
eps_values = [0.00001, 0.0001, 0.001, 0.01]

for metric in metrics:
    for min_cluster_size in min_cluster_sizes:
        for eps in eps_values:
            hdbscan = HDBSCAN(
                metric=metric,
                min_cluster_size=min_cluster_size, 
                cluster_selection_epsilon=eps, 
            )
            clusters = hdbscan.fit_predict(X)
            n_clusters = len(set(clusters)) - (1 if -1 in clusters else 0)
            print(f"metric={metric}, min_cluster_size={min_cluster_size}, eps={eps}, n_clusters={n_clusters}")

            # Aggiorna la configurazione con il numero massimo di cluster
            if n_clusters > max_clusters:
                max_clusters = n_clusters
                best_configuration = {
                    'metric': metric,
                    'min_cluster_size': min_cluster_size,
                    'eps': eps,
                    'n_clusters': n_clusters
                }

metric=euclidean, min_cluster_size=3, eps=1e-05, n_clusters=137
metric=euclidean, min_cluster_size=3, eps=0.0001, n_clusters=137
metric=euclidean, min_cluster_size=3, eps=0.001, n_clusters=108
metric=euclidean, min_cluster_size=3, eps=0.01, n_clusters=2
metric=euclidean, min_cluster_size=5, eps=1e-05, n_clusters=65
metric=euclidean, min_cluster_size=5, eps=0.0001, n_clusters=65
metric=euclidean, min_cluster_size=5, eps=0.001, n_clusters=58

metric=euclidean, min_cluster_size=5, eps=0.01, n_clusters=2
metric=euclidean, min_cluster_size=10, eps=1e-05, n_clusters=23
metric=euclidean, min_cluster_size=10, eps=0.0001, n_clusters=23
metric=euclidean, min_cluster_size=10, eps=0.001, n_clusters=23
metric=euclidean, min_cluster_size=10, eps=0.01, n_clusters=2
metric=manhattan, min_cluster_size=3, eps=1e-05, n_clusters=136
metric=manhattan, min_cluster_size=3, eps=0.0001, n_clusters=136

metric=manhattan, min_cluster_size=3, eps=0.001, n_clusters=113
metric=manhattan, min_cluster_size=3, eps=0.01, n_clusters=2
metric=manhattan, min_cluster_size=5, eps=1e-05, n_clusters=66
metric=manhattan, min_cluster_size=5, eps=0.0001, n_clusters=66

metric=manhattan, min_cluster_size=5, eps=0.001, n_clusters=60
metric=manhattan, min_cluster_size=5, eps=0.01, n_clusters=2
metric=manhattan, min_cluster_size=10, eps=1e-05, n_clusters=23
metric=manhattan, min_cluster_size=10, eps=0.0001, n_clusters=23
metric=manhattan, min_cluster_size=10, eps=0.001, n_clusters=23
metric=manhattan, min_cluster_size=10, eps=0.01, n_clusters=2
metric=haversine, min_cluster_size=3, eps=1e-05, n_clusters=124

metric=haversine, min_cluster_size=3, eps=0.0001, n_clusters=116
metric=haversine, min_cluster_size=3, eps=0.001, n_clusters=8
metric=haversine, min_cluster_size=3, eps=0.01, n_clusters=2

metric=haversine, min_cluster_size=5, eps=1e-05, n_clusters=56
metric=haversine, min_cluster_size=5, eps=0.0001, n_clusters=56
metric=haversine, min_cluster_size=5, eps=0.001, n_clusters=7

metric=haversine, min_cluster_size=5, eps=0.01, n_clusters=2
metric=haversine, min_cluster_size=10, eps=1e-05, n_clusters=16
metric=haversine, min_cluster_size=10, eps=0.0001, n_clusters=16

metric=haversine, min_cluster_size=10, eps=0.001, n_clusters=4
metric=haversine, min_cluster_size=10, eps=0.01, n_clusters=2

print(f"Miglior configurazione: {best_configuration}")

Miglior configurazione: {'metric': 'euclidean', 'min_cluster_size': 3, 'eps': 1e-05, 'n_clusters': 137}

hdbscan = HDBSCAN(
    metric=best_configuration['metric'],
    min_cluster_size=best_configuration['min_cluster_size'], 
    cluster_selection_epsilon=best_configuration['eps'], 
)
hdbscan_labels = hdbscan.fit_predict(X)
results['hdbscan'] = {
    'labels': hdbscan_labels,
    'n_clusters': len(set(hdbscan_labels)) - (1 if -1 in hdbscan_labels else 0)
}

visualize_clusters(X, results['hdbscan']['labels'])

Make this Notebook Trusted to load map: File -> Trust Notebook

Conclusioni

Da un'ispezione visiva dei cluster emerge che il miglior modello è l'hdb-scan, perché consente di non classificare tutti i punti (va bene che ci siano punti che non fanno parte di alcun cluster se sono isolati)

attrezzatura_sportiva.head(2)

	geometry	Nome	ClassID	Lotto	Strato	Classe	ClassREF	Tipo	index_right	cod_reg	...	provincia	regione	area_km2	TIME_PERIOD	females	males	total	geometry_c	x	y
381185	POLYGON Z ((9.12812 39.25162 0, 9.12809 39.251...	None	ATTR_SP_944	CU1	02	ATTREZZATURA_SPORTIVA	None	gradinata	7688	20	...	Cagliari	Sardegna	84.823594	2024	78689	68689	147378	POINT (9.12809 39.25163)	9.128088	39.251628
381186	POLYGON Z ((9.12687 39.25094 0, 9.1268 39.2508...	None	ATTR_SP_945	CU1	02	ATTREZZATURA_SPORTIVA	None	gradinata	7688	20	...	Cagliari	Sardegna	84.823594	2024	78689	68689	147378	POINT (9.12683 39.25093)	9.126833	39.250926

2 rows × 25 columns

attrezzatura_sportiva['clusters'] = hdbscan_labels

attrezzatura_sportiva.to_parquet(os.path.join(cu_path,'attrezzatura_sportiva_clustered.parquet'))