La respiration de Paris II : GeoPandas pour l'analyse de données spatiales

Dans un précédent billet, je vous ai montré comment lancer une acquisition des données Vélib en utilisant Python et la librairie Pandas et comment effectuer des manipulations simples. Commençons par charger ces données et par faire une jointure entre nos deux tables. Si vous avez besoin d'un rappel, vous pouvez consulter le précédent billet.

import pandas as pd

data_stations = pd.read_sql_table("stations", "sqlite:///database.db")
data_localisation = pd.read_sql_table("localisation", "sqlite:///database.db")
data_stations = data_stations.merge(data_localisation, on="station_id")

data_stations

GeoDataFrame : un dataframe, mais géographique

Aujourd'hui, nous allons analyser ces données spatiales en utilisant la librairie GeoPandas. Cette librairie présente beaucoup de similarité avec Pandas. L'objet de base s'appelle un GeoDataFrame, il s'agit d'un DataFrame avec une colonne spéciale nommée geometry qui contient des renseignements géographiques. Cette geometry peut être un point, comme dans le cas présent, mais aussi une ligne pour désigner une frontière ou bien un polygone pour matérialiser un territoire. Nous reviendrons sur ce dernier point dans quelques instants.

La première étape consiste à convertir notre couple (longitude, latitude) en un object Point reconnu par GeoPandas, puis à créer notre GeoDataFrame.

import geopandas as gpd

# Convert the longitude and latitude to a format recognized by geoPandas
geometry = gpd.points_from_xy(data_stations["lon"], data_stations["lat"])

# Create a DataFrame with a geometry containing the Points
geo_stations = gpd.GeoDataFrame(
    data_stations, crs="EPSG:4326", geometry=geometry
)

geo_stations

L'acronyme crssignifie Coordinate Reference System, c'est une indication du système de projection utilisé. En regardant la documentation Vélib, on voit que le référentiel de projection utilisé est WGS84. C'est le système de projection le plus commun aujourd'hui et il est notamment utilisé par les systèmes de positionnement par satellite GPS. Ce système est référencé 4326 en deux dimensions (X,Y) et 4979 en trois dimensions (X,Y,Z) selon la liste des codes EPSG, et c'est ce que nous donnons comme indication à notre GeoDataFrame.

Nos premières cartes

Avant de faire une carte, nous allons sélectionner des données à un temps fixé, et les enregistrer dans un nouveau GeoDataFrame

some_time = geo_stations["time_stamp"][0]
some_data = geo_stations[geo_stations["time_stamp"] == some_time]

some_time

Timestamp('2021-06-21 17:26:27.825403')

À ce stade, GeoPandas connait la geometry de nos données, c'est-à-dire comment elles s'agencent dans l'espace. Il ne nous reste plus qu'à afficher les données

some_data.plot()

Notre première carte n'a pas beaucoup d'allure, mais je vous assure que le plus gros du travail est derrière nous ! Pour l'instant, GeoPandas ne fait qu'afficher un marqueur bleu à l'emplacement de chaque station Vélib. Demandons à GeoPandas de changer la taille de nos marqueurs en fonction du nombre de vélibs dans la station

some_data.plot(markersize="num_bikes_available")

On se représente déjà mieux la distribution du nombre de vélibs dans Paris ! Mais au lieu de jouer avec la taille des marqueurs, on peut aussi varier leur couleur en spécifiant une carte de couleur cmap. Pour montrer une concentration, j'apprécie beaucoup la ColorMap OrRd présente par défaut dans Matplotlib et Pandas.

some_data.plot("num_bikes_available", cmap="OrRd")

On peut même combiner ces deux effets pour être encore plus explicite !

some_data.plot("num_bikes_available", markersize="num_bikes_available", cmap="OrRd")

Une belle mise en forme

Jusque-là nous avons laissé faire GeoPandas tout seul, mais il va falloir l'aider pour obtenir de belles cartes. Dans un premier temps, on va utiliser matplotlib pour reprendre la main sur la mise en forme. On va pouvoir ajouter un titre, déterminer la taille de la carte et enlever ces horribles axes

import matplotlib.pyplot as plt

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

ax.set_title(some_time.strftime("%A %B %d %H:%M"))
some_data.plot("num_bikes_available", markersize="num_bikes_available", cmap="OrRd", ax=ax)
ax.set_axis_off()


plt.show()

Il y a du mieux ! Maintenant nous allons pouvoir ajouter un fond de carte provenant, par exemple, de OpenStreetMap. Il existe un module python pour faire ça en une ligne et cette librairie magique s'appelle contextily. Le seul pré-requis est de convertir les coordonnées GPS au format EPSG:3857, mais GeoPandas fait ça très bien

import contextily as ctx

# Conversion des coordonnées
some_data = some_data.to_crs(epsg=3857)

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

ax.set_title(some_time.strftime("%A %B %d %H:%M"))
some_data.plot("num_bikes_available", markersize="num_bikes_available", cmap="OrRd", ax=ax)
ax.set_axis_off()

# Ajout du fond de carte
ctx.add_basemap(ax)

plt.show()

J'espère vous avoir convaincu de la puissance de GeoPandas !

Aller plus loin : rassembler les données par quartiers

L'inconvénient majeur de nos cartes est qu'il est difficile de distinguer et de hiérarchiser les zones à forte concentration de vélibs. La raison en est que les marqueurs ont tendance à se chevaucher et que l'on perd en lisibilité. Pour remédier à ce problème, nous allons rassembler les stations appartenant à un même quartier. Ainsi, nous allons tracer des surfaces sur la carte au lieu d'un nuage de points. La liste des quartiers administratifs de Paris est mise à disposition sur le site de la ville de Paris. GeoPandas sait parfaitement importer les fichiers GeoJSON

districts = gpd.read_file("quartier_paris.geojson").to_crs(epsg=3857)
districts

Chaque quartier est donc défini comme un polygone, avec un certain périmètre et une surface. On peut directement tracer les quartiers de Paris sur la carte.

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

districts.plot(ax=ax, alpha=0.5, edgecolor="white")
ax.set_axis_off()
ctx.add_basemap(ax)

plt.show()

On aimerait disposer du nombre de vélibs présent dans chacun de ces quartiers, et c'est là que se déploie toute la puissance de GeoPandas, voyez plutôt

districts["velib_number"] = districts.apply(
    lambda district: (
        some_data.within(district.geometry) * some_data["num_bikes_available"]
    ).sum(),
    axis=1
)

Ce bloc mérite quelques explications. On crée une nouvelle colonne nommée velib_number, c'est à dire le nombre de vélibs présents dans le quartier. Pour ce faire, on somme tous les vélos qui sont inclus dans le quartier en question. Cette opération est rendue possible par la méthode .within() qui va déterminer si la station est bien localisée dans le polygone représentant le quartier.

On peut faire la même opération pour calculer le nombre maximum de vélibs dans le quartier, et donc le pourcentage d'occupation des stations dans le quartier

districts["max_velib_number"] = districts.apply(
        lambda district: (
            some_data.within(district.geometry) * some_data["capacity"]
        ).sum(),
        axis=1,
)

districts["occupation"] = districts["velib_number"] / districts["max_velib_number"]

Enfin, on va tracer le taux d'occupation des stations Vélib, quartier par quartier

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

districts.plot("occupation", cmap="OrRd", ax=ax, alpha=0.5, edgecolor="white")
ax.set_axis_off()
ctx.add_basemap(ax)

plt.show()

Moteur et action !

Maintenant que vous savez produire une carte, il n'y a plus d'obstacle sur votre route pour en obtenir quelques milliers et produire une vidéo ! Je vous épargne ces dernières étapes et vous laisse contempler le résultat de notre analyse ci-dessous :

On peut distinguer à l'oeil quelques tendances suivant les quartiers, certains sont actifs la nuit et d'autres le jour. Mais est-il possible de faire mieux qu'une simple analyse à l'oeil ? Ce sera l'objet du dernier épisode de cette série sur les données Vélibs.