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

24 June 2021

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
index_x stationCode station_id num_bikes_available num_docks_available time_stamp index_y name lat lon capacity
0 0 16107 213688169 1 34 2021-06-21 17:26:27.825403 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35
1 0 16107 213688169 1 34 2021-06-21 17:27:26.935454 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35
2 0 16107 213688169 1 34 2021-06-21 17:28:26.940625 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35
3 0 16107 213688169 1 34 2021-06-21 17:29:26.936032 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35
4 0 16107 213688169 1 34 2021-06-21 17:30:26.935873 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35
... ... ... ... ... ... ... ... ... ... ... ...
2065141 1403 9104 129026597 7 13 2021-06-22 18:16:55.748422 1401 Caumartin - Provence 48.874423 2.328469 22
2065142 1403 9104 129026597 5 15 2021-06-22 18:17:55.748422 1401 Caumartin - Provence 48.874423 2.328469 22
2065143 1403 9104 129026597 5 15 2021-06-22 18:18:55.748371 1401 Caumartin - Provence 48.874423 2.328469 22
2065144 1403 9104 129026597 6 14 2021-06-22 18:19:55.748149 1401 Caumartin - Provence 48.874423 2.328469 22
2065145 1403 9104 129026597 5 15 2021-06-22 18:20:55.748413 1401 Caumartin - Provence 48.874423 2.328469 22

2065146 rows × 11 columns

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
index_x stationCode station_id num_bikes_available num_docks_available time_stamp index_y name lat lon capacity geometry
0 0 16107 213688169 1 34 2021-06-21 17:26:27.825403 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35 POINT (2.27572 48.86598)
1 0 16107 213688169 1 34 2021-06-21 17:27:26.935454 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35 POINT (2.27572 48.86598)
2 0 16107 213688169 1 34 2021-06-21 17:28:26.940625 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35 POINT (2.27572 48.86598)
3 0 16107 213688169 1 34 2021-06-21 17:29:26.936032 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35 POINT (2.27572 48.86598)
4 0 16107 213688169 1 34 2021-06-21 17:30:26.935873 0 Benjamin Godard - Victor Hugo 48.865983 2.275725 35 POINT (2.27572 48.86598)
... ... ... ... ... ... ... ... ... ... ... ... ...
2065141 1403 9104 129026597 7 13 2021-06-22 18:16:55.748422 1401 Caumartin - Provence 48.874423 2.328469 22 POINT (2.32847 48.87442)
2065142 1403 9104 129026597 5 15 2021-06-22 18:17:55.748422 1401 Caumartin - Provence 48.874423 2.328469 22 POINT (2.32847 48.87442)
2065143 1403 9104 129026597 5 15 2021-06-22 18:18:55.748371 1401 Caumartin - Provence 48.874423 2.328469 22 POINT (2.32847 48.87442)
2065144 1403 9104 129026597 6 14 2021-06-22 18:19:55.748149 1401 Caumartin - Provence 48.874423 2.328469 22 POINT (2.32847 48.87442)
2065145 1403 9104 129026597 5 15 2021-06-22 18:20:55.748413 1401 Caumartin - Provence 48.874423 2.328469 22 POINT (2.32847 48.87442)

2065146 rows × 12 columns

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
n_sq_qu n_sq_ar c_qu surface l_qu perimetre c_quinsee c_ar geometry
0 750000006 750000002 6 2.435508e+05 Vivienne 2058.472959 7510202 2 POLYGON ((260624.854 6252121.461, 260618.773 6...
1 750000010 750000003 10 2.717503e+05 Enfants-Rouges 2139.625388 7510302 3 POLYGON ((263504.516 6251415.017, 263522.754 6...
2 750000024 750000006 24 2.822799e+05 Saint-Germain-des-Prés 2565.899893 7510604 6 POLYGON ((260149.190 6250092.466, 260152.388 6...
3 750000037 750000010 37 9.268652e+05 Saint-Vincent-de-Paul 4072.789633 7511001 10 POLYGON ((262771.159 6253845.649, 262737.400 6...
4 750000042 750000011 42 8.379929e+05 Saint-Ambroise 4052.567737 7511102 11 POLYGON ((263931.799 6250766.570, 263629.859 6...
... ... ... ... ... ... ... ... ... ...
75 750000020 750000005 20 4.331978e+05 Sorbonne 2892.944068 7510504 5 POLYGON ((261516.706 6248520.128, 261476.634 6...
76 750000033 750000009 33 7.170916e+05 Saint-Georges 3429.188334 7510901 9 POLYGON ((260401.731 6253981.490, 260282.394 6...
77 750000034 750000009 34 5.434412e+05 Chaussée-d'Antin 3133.580092 7510902 9 POLYGON ((259981.114 6253987.888, 260164.693 6...
78 750000003 750000001 3 2.736968e+05 Palais-Royal 2166.839239 7510103 1 POLYGON ((260428.149 6251500.935, 260390.032 6...
79 750000048 750000012 48 1.235916e+06 Quinze-Vingts 4509.486974 7511204 12 POLYGON ((264183.946 6247852.678, 264095.921 6...

80 rows × 9 columns

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.