La respiration de Paris III : classification automatique non supervisée

Ce billet est le dernier de notre série sur les données publiques Vélib. Dans un premier temps, nous avons vu comment acquérir et enregistrer ces données avec Python et Pandas. Et dans un second temps, comment représenter ces données graphiquement avec l'aide de GeoPandas. Ce faisant, nous avons pu "voir" Paris respirer pendant $24$ heures consécutives et identifier à l'oeil nu les différents types de quartier dans la capitale.

Aujourd'hui, je vous propose une approche complètement automatique pour classifier les différents types de quartier dans Paris, en utilisant l'algorithme des $k$-moyennes.

Mais ne brûlons pas les étapes, commençons par importer nos précieuses données. Je vous renvoie aux deux précédents billets pour comprendre ce qui se passe ici.

import pandas as pd
import geopandas as gpd

# Importation des données
districts = gpd.read_file("quartier_paris.geojson")
data_localisation = pd.read_sql_table("localisation", "sqlite:///database.db")
data_observations = pd.read_sql_table("stations", "sqlite:///database.db")

# Fusions des tables localisation et observations
data_merged = data_observations.merge(data_localisation, on="station_id")

# On interprète ces données comme des données géographiques
geo_stations = gpd.GeoDataFrame(
    data_merged,
    crs="EPSG:4326",
    geometry=gpd.points_from_xy(data_merged["lon"], data_merged["lat"])
)

# On enlève les colonnes inutiles
geo_stations.drop(columns=["stationCode", "name", "lon", "lat", "index_x", "index_y"], inplace=True)
districts.drop(columns=["n_sq_qu", "n_sq_ar", "c_qu", "surface", "perimetre", "c_quinsee", "c_ar"], inplace=True)

# Notre table contenant les stations
geo_stations.head()

Séries temporelles par quartier

Le réseau Vélib compte environ $1400$ stations et il nous serait difficile de les étudier toutes en même temps. Nous commençons donc par regrouper les données par quartier et par heure.

# Jointure spatiale : on identifie le quartier de chaque observation
geo_districts = gpd.sjoin(districts, geo_stations, how="left")

# On aggrège les données par quartier et par heure
occupation_data = geo_districts.groupby(["l_qu", "time_stamp"], as_index=False).agg({"num_bikes_available":"sum", "geometry": "first"})

occupation_data.head()

La syntaxe est ici plus compacte que dans le précédent épisode donc nous allons la décortiquer ensemble. Dans un premier, nous groupons les données en fonction de l_qu, le nom du quartier et de time_stamp, c'est-à-dire l'heure d'observation. À l'issue de cette opération, on obtient un objet DataFrameGroupBy. Notez l'option as_index=False qui est nécessaire ici.

Puis nous aggrégeons les données avec la méthode .agg(). En argument, on lui dit qu'on veut sommer le nombre de vélos disponibles dans chaque groupe, et garder la géométrie associée au quartier. On peut ne garder que la première occurence, puisque toutes ces géométries sont identiques.

Nous avons donc l'historique, minute par minute, du nombre de vélibs stationnés dans chaque quartier, que nous pouvons visualiser de la façon suivante

import matplotlib.pyplot as plt
import matplotlib.dates as mdates

fig, ax = plt.subplots(figsize=(8,6))

occupation_data.groupby('l_qu').plot(x="time_stamp", y="num_bikes_available", kind="line", ax=ax, legend=False)
ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))

plt.show()

Chacune de ces courbes est une série temporelle qui caractérise l'activité du réseau Vélib dans son quartier. Ce diagramme est un peu brouillon, mais on peut distinguer à l'oeil quelques tendances. Certaines stations semblent se vider en journée quand d'autres semblent se remplir dans le même interval. Mais peut-on faire mieux et de manière automatique ?

Algorithme des $k$-moyennes

La question que l'on se pose est la suivante : "Est-il possible de classer ces séries temporelles de manières automatique ?". La réponse est oui, et pour cela nous allons utiliser l'algorithme des $k$-moyennes. Cet algorithme permet de diviser des données en $k$ partitions. Sa particularité est de pouvoir être mis en place sans supervisition. C'est à dire qu'il va apprendre tout seul à faire la distinction entre les différentes partitions sans que nous ayons besoin d'intervenir, et c'est précisément ce qu'on lui demande.

D'un point de vue formel, nous commençons avec un ensemble de points $(x_1, x_2, ..., x_n)$ que l'on cherche à partionner en $k$ ensembles notés $(S_1, ... S_k)$. L'algorithme des $k$-moyennes est une heuristique ayant pour but de minimiser $$ \sum_{i=1}^k \sum_{x_j \in S_i} || x_j - \mu_i ||^2, $$ où $\mu_i$ est le barycentre de $S_i$. Dit autrement, on veut construire $k$ ensembles disjoints au sein desquels les points sont les plus groupés. C'est un problème d'optimisation compliqué en théorie et l'algorithme des $k$-moyennes est une méthode efficace pour trouver une décomposition, mais il ne garantit pas que ce soit la solution optimale. Cette décomposition sera tout de même suffisante pour notre cas d'étude.

L'algorithme fonctionne de la manière suivante :

• On choisit au hasard les barycentres $\mu_i^{(1)}$ des différentes partitions $S_i$
• On répète le processus suivant jusqu'à convergence :
  • On affecte chaque point $x_j$ à la partition la plus proche, ce qui permet de définir les partitions $$ S_i^{(t)} = \left\{ x_j: \forall \ell,~ ||x_j - \mu_i(t)|| \leq ||x_j - \mu_\ell ^{(t)} || \right\} $$
  • On met à jour le barycentre de chacune des partitions $$ \mu_i^{(t+1)} = \frac{1}{|S_i^{(t)}|} \sum_{x_j \in S_i^{(t)}} x_j. $$

Nous utiliserons la librairie tslearn, littéralement "apprentissage pour les séries temporelles", qui repose notamment sur sklearn, une librairie très connue de machine learning. Elle va nous permettre d'appliquer l'algorithme des $k$-moyennes où chaque $x_i$ correspond à une de nos séries temporelles et où la norme $||x_j - \mu_i(t)||$ est la norme euclidienne. La norme euclidienne veut simplement dire que l'on compare nos deux séries terme à terme, soit à même temps. Il existe des normes plus élaborées comme la norme DTW pour "Dynamic Time Warping" qui sont fréquemment utilisées dans l'analyse des séries temporelles, et sont surtout adaptées à des séries temporelles qui ne seraient pas parfaitement synchrones. Dans le cas présent, il serait contre-productif de choisir une telle norme, d'une part parce qu'elle est beaucoup plus lente à déterminer, et d'autre part parce que nous souhaitons avoir des classes parfaitement synchronisées.

Applications aux données temporelles

Commençons l'analyse en mettant nos données sous la forme de séries temporelles pour la librairie tslearn

from tslearn.utils import to_time_series_dataset
from tslearn.preprocessing import TimeSeriesScalerMeanVariance

labels, time_series = [] , []
for l_qu, group in occupation_data.groupby("l_qu"):
    labels.append(l_qu)
    time_series.append(group["num_bikes_available"].array)

# On formatte les séries temporelles
time_series = to_time_series_dataset(time_series)
# On normalise ces mêmes séries
time_series = TimeSeriesScalerMeanVariance().fit_transform(time_series)

time_series.shape

(80, 1473, 1)

Rien de particulier ici, si ce n'est que la fonction to_time_series_dataset crée un jeu de données pour tslearn à partir d'un tableau, et que nous normalisons ces données avec l'aide de TimeSeriesScalerMeanVariance qui fait en sorte que chaque série ait pour moyenne $0$ et variance $1$. Cette étape est très importante puisqu'elle permet de s'affranchir de nombreux biais comme le nombre absolu de vélibs par quartier, ou la proportion d'habitants utilisant le réseau.

On peut maintenant lancer l'apprentissage non-supervisé de notre modèle avec l'algorithme des $k$-moyennes

from tslearn.clustering import TimeSeriesKMeans
n_classes = 3

model = TimeSeriesKMeans(n_clusters=n_classes, metric="euclidean")
model.fit(time_series)

TimeSeriesKMeans()

Visualisons les barycentres (en rouge) et les différents éléments (gris) de chaque classe que notre modèle a trouvé

import matplotlib.dates as mdates

fig = plt.figure(figsize=(8,n_classes * 3))

# On récupère les dates des différentes observations
time_labels = occupation_data["time_stamp"].unique()

# Pour chaque classe
for yi in range(n_classes):
    ax = fig.add_subplot(n_classes, 1, 1 + yi)

    # On sélectionne les séries qui correspondent à cette classe
    for xx in time_series[model.labels_ == yi]:
        ax.plot(time_labels, xx, "k-", alpha=.2)

    # Le barycentre de la classe
    ax.plot(time_labels, model.cluster_centers_[yi].ravel(), "r-")

    # Pour formatter l'heure des observations
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))

plt.show()

Interprétation et visualisation

Notre algortihme a classé automatiquement tous les quartiers en trois catégories bien distinctes que l'on peut interpréter comme étant :

1. Les quartiers d'affaires qui se remplissent à 9h du matin et se vident aux alentours de 18h
2. Les quatiers résidentiels qui se remplissent durant la nuit et se vident à 9h
3. Les quartiers d'activités nocturnes qui sont actifs de 18h à minuit environ

Maintenant que nous avons déterminé la classe de chaque quartier, nous pouvons l'ajouter à notre liste des quartiers de Paris

interpretation = ["Affaire", "Résidentiel", "Nocturne"]

classification = pd.DataFrame({"l_qu":labels, "type": [interpretation[classe] for classe in model.labels_] })
districts_classified = districts.merge(classification, on="l_qu")

districts_classified.head()

Et produire la carte des types de quartiers dans Paris à partir des données Vélibs

import contextily as ctx

fig, ax = plt.subplots(figsize=(8, 6))

districts_classified = districts_classified.to_crs(epsg=3857)

districts_classified.plot("type", ax=ax, alpha=0.5, edgecolor="white", legend=True)
ax.set_axis_off()
ctx.add_basemap(ax)

plt.show()

Je vous laisse le soin de me dire ce que vous en pensez !

Cette carte marque la fin de cette série sur les données publiques du réseau Vélib. Cette série a été pour moi un chouette exercice pour mettre en application plein d'outils différents

• Pandas pour la gestion des données
• Request pour aller interroger différentes API
• APScheduler pour programmer des tâches à intervalles réguliers
• GeoPandas pour la manipulation de données géographiques et la production de cartes
• Matplotlib pour la mise en forme de graphiques et de cartes
• Tslearn qui est une porte d'entrée pour s'intéresser au machine learning.

J'espère que vous avez apprécié, je vous souhaite une bonne journée !