La cartographie est un excellent moyen de diffuser
une connaissance, y compris à des publics peu
familiers de la statistique. Ce chapitre permet
de découvrir la manière dont on peut
utiliser Python pour construire des
cartes standards (avec geopandas) ou
réactives (folium). Cela se fera
à travers un exercice permettant
de visualiser la fréquentation par les
vélos des routes parisiennes.
La pratique de la cartographie se fera, dans ce cours, en répliquant des cartes qu’on peut trouver sur
la page de l’open-data de la ville de Paris
ici.
Note
Produire de belles cartes demande du temps mais aussi du bon sens. En fonction de la structure des données, certaines représentations sont à éviter voire à exclure. L’excellent guide disponible ici propose quelques règles et évoque les erreurs à éviter lorsqu’on désire effectuer des
représentations spatiales.
Celui-ci reprend un guide de sémiologie cartographique
produit par l’Insee qui propose de nombreux conseils pratiques pour produire des représentations
cartographiques sensées.
Ce TP vise à initier :
Au module graphique de geopandas ainsi qu’aux packages geoplot et
contextily pour la construction de cartes figées. geoplot est construit sur seaborn et constitue ainsi une extension des graphiques de base.
Au package folium qui est un point d’accès vers la librairie JavaScript leaflet permettant de produire des cartes interactives
Les données utilisées sont :
Un sous-ensemble des données de paris open data a été mis à disposition
sur pour faciliter l’import (élimination des colonnes
qui ne nous servirons pas mais ralentissent l’import)
Certaines librairies géographiques dépendent de rtree qui est parfois difficile à installer car
ce package dépend de librairies compilées qui sont compliquées à installer sur Windows.
Pour installer rtree sur Windows, le mieux est d’utiliser Anaconda.
Avant de pouvoir commencer, il est nécessaire d’installer quelques
packages au préalable:
# Sur colab!pip install pandas fiona shapely pyproj rtree # à faire obligatoirement en premier pour utiliser rtree ou pygeos pour les jointures spatiales!pip install contextily!pip install geopandas!pip install geoplot
Dans la première partie, nous allons utiliser les packages suivants :
import pandas as pdimport geopandas as gpdimport contextily as ctximport geoplotimport matplotlib.pyplot as pltimport folium
/opt/mamba/lib/python3.9/site-packages/geopandas/_compat.py:123: UserWarning:
The Shapely GEOS version (3.12.0-CAPI-1.18.0) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
/tmp/ipykernel_4003/3556787358.py:2: UserWarning:
Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:
import os
os.environ['USE_PYGEOS'] = '0'
import geopandas
In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
Première carte avec l’API matplotlib de geopandas
Exercice 1: Importer les données
Importer les données de compteurs de vélos en deux temps.
Faire attention à deux valeurs aberrantes. Utiliser
la fonctionalité str.contains pour exclure les
observations contenant “Bike IN” ou “Bike OUT”
dans la variable
nom_compteur
Utiliser la méthode plot pour représenter les localisations des compteurs dans l’espace. C’est, on peut l’avouer, peu informatif sans apport extérieur. Il va donc falloir travailler un peu l’esthétique
/tmp/ipykernel_4003/3602001210.py:1: UserWarning:
This pattern is interpreted as a regular expression, and has match groups. To actually get the groups, use str.extract.
Warning
On serait tenté de faire un merge de la base compteurs et comptages.
En l’occurrence, il s’agirait d’un produit cartésien puisqu’il s’agit de faire exploser la base spatiale.
Avec des données spatiales, c’est souvent une très mauvaise idée. Cela duplique les points, créant des difficultés à représenter les données mais aussi ralentit les calculs.
Sauf à utiliser la méthode dissolve (qui va agréger k fois la même géométrie…), les géométries sont perdues lorsqu’on effectue des groupby.
Maintenant, tout est prêt pour une première carte. matplotlib fonctionne selon
le principe des couches. On va de la couche la plus lointaine à celle le plus
en surface. L’exception est lorsqu’on ajoute un fond de carte contextily via
ctx.add_basemap: on met cet appel en dernier.
Exercice 2: Première carte
Représenter une carte des compteurs avec le fonds de carte des arrondissements
Faire attention à avoir des arrondissements dont l’intérieur est transparent (argument à utiliser: facecolor).
Faire des bordures d’arrondissements noires et affichez les compteurs en rouge.
Pour obtenir un graphique plus grand, vous pouvez utiliser l’argument figsize = (10,10).
Pour les localisations, les points doivent être rouges en étant plus transparent au centre (argument à utiliser: alpha)
Vous devriez obtenir cette carte:
Exercice 3 : Ajouter un fonds de carte avec contextily
Repartir de la carte précédente.
Utiliser ctx.add_basemap pour ajouter un fonds de carte. Pour ne pas afficher les axes, vous pouvez utiliser ax.set_axis_off().
:warning: Par défaut, contextily désire un système de projection (crs) qui est le Web Mercator (epsg: 3857). Il faut changer la valeur de l’argument crs.
:warning: Avec les versions anciennes des packages, il faut utiliser .to_string sur un objet CRS pour qu’il soit reconnu par contextily. Sur des versions récentes, la valeur numérique du code EPSG est suffisante.
Trouver un fonds de carte plus esthétique, qui permette de visualiser les grands axes, parmi ceux possibles. Pour tester l’esthétique, vous pouvez utiliser cet url. La documentation de référence sur les tuiles disponibles est ici
ax.get_figure()
ax.get_figure()
ax.get_figure().savefig("featured_maps.png")
Le principe de la heatmap est de construire, à partir d’un nuage de point bidimensionnel, une distribution 2D lissée. La méthode repose sur les estimateurs à noyaux qui sont des méthodes de lissage local.
Exercice 3 : Ajouter un fonds de carte avec contextily
Pour le moment, la fonction geoplot.kdeplot n’incorpore pas toutes les fonctionalités de seaborn.kdeplot. Pour être en mesure de construire une heatmap avec des données pondérées (cf. cette issue dans le dépôt seaborn), il y a une astuce. Il faut simuler k points de valeur 1 autour de la localisation observée. La fonction ci-dessous, qui m’a été bien utile, est pratique
import numpy as npdef expand_points(shapefile, index_var ="grid_id", weight_var ='prop', radius_sd =100, crs =2154):""" Multiply number of points to be able to have a weighted heatmap :param shapefile: Shapefile to consider :param index_var: Variable name to set index :param weight_var: Variable that should be used :param radius_sd: Standard deviation for the radius of the jitter :param crs: Projection system that should be used. Recommended option is Lambert 93 because points will be jitterized using meters :return: A geopandas point object with as many points by index as weight """ shpcopy = shapefile shpcopy = shpcopy.set_index(index_var) shpcopy['npoints'] = np.ceil(shpcopy[weight_var]) shpcopy['geometry'] = shpcopy['geometry'].centroid shpcopy['x'] = shpcopy.geometry.x shpcopy['y'] = shpcopy.geometry.y shpcopy = shpcopy.to_crs(crs) shpcopy = shpcopy.loc[np.repeat(shpcopy.index.values, shpcopy.npoints)] shpcopy['x'] = shpcopy['x'] + np.random.normal(0, radius_sd, shpcopy.shape[0]) shpcopy['y'] = shpcopy['y'] + np.random.normal(0, radius_sd, shpcopy.shape[0]) gdf = gpd.GeoDataFrame( shpcopy, geometry = gpd.points_from_xy(shpcopy.x, shpcopy.y), crs = crs)return gdf
Exercice 4 : Data cleaning avant de pouvoir faire une heatmap
Calculer le trafic moyen, pour chaque station, entre 7 heures et 10 heures (bornes incluses) et nommer cet objet df1. Faire la même chose, en nommant df2, pour le trafic entre 17 et 20 heures (bornes incluses)
Nous allons désormais préparer les données de manière à faire une heatmap. Après avoir compris ce que permet de faire la fonction expand_points ci-dessus, créer une fonction explode_data qui suive les étapes suivantes.
Convertir un DataFrame dans le système de projection Lambert 93 (epsg: 2154)
Appliquer expand_points aux noms de variable adéquats. Vous pouvez fixer la valeur de radius_sd à 100.
Reconvertir l’output au format WGS84 (epsg: 4326)
Appliquer cette fonction à df1 et df2
Exercice 5 : Heatmap, enfin !
Représenter, pour ces deux moments de la journée, la heatmap du trafic de vélo avec geoplot.kdeplot. Pour cela :
Appliquer la fonction geoplot.kdeplot avec comme consignes :
d’utiliser les arguments shade=True et shade_lowest=True pour colorer l’intérieur des courbes de niveaux obtenues ;
d’utiliser une palette de couleur rouge avec une transparence modérée (alpha = 0.6)
d’utiliser l’argument clip pour ne pas déborder hors de Paris (en cas de doute, se référer à l’aide de geoplot.kdeplot)
L’argument bw (pour bandwidth) détermine le plus ou moins fort lissage spatial. Vous pouvez partir d’un bandwidth égal à 0.01 et le faire varier pour voir l’effet sur le résultat
Ne pas oublier d’ajouter les arrondissements. Avec geoplot, il faut utiliser geoplot.polyplot.
ax.get_figure()
Des cartes réactives grâce à folium
De plus en plus de données de visualisation reposent sur la cartographie réactive. Que ce soit dans l’exploration des données ou dans la représentation finale de résultats, la cartographie réactive est très appréciable.
folium offre une interface très flexible et très facile à prendre à main. Les cartes sont construites grâce à la librairie JavaScript Leaflet.js mais, sauf si on désire aller loin dans la customisation du résultat, il n’est pas nécessaire d’avoir des notions dans le domaine.
Un objet folium se construit par couche. La première est l’initialisation de la carte. Les couches suivantes sont les éléments à mettre en valeur. L’initialisation de la carte nécessite la définition d’un point central (paramètre location) et d’un zoom de départ (zoom_start). Plutôt que de fournir manuellement le point central et le zoom on peut :
Déterminer le point central en construisant des colonnes longitudes et latitudes et en prenant la moyenne de celles-ci ;
Utiliser la méthode fit_bounds qui cale la carte sur les coins sud-ouest et nord-est. En supposant que la carte s’appelle m, on fera m.fit_bounds([sw, ne])
Le bout de code suivant permet de calculer le centre de la carte
Si un fond gris s’affiche, c’est qu’il y a un problème de localisation ou d’accès à internet. Pour le premier cas, cela provient généralement d’un problème de projection ou d’une inversion des longitudes et latitudes.
Les longitudes représentent les x (axe ouest-est) et les latitudes y (axe sud-nord). De manière contrintuitive, folium attend qu’on lui fournisse les données sous la forme [latitude, longitude] donc [y,x]
Exercice 6 : Visualiser la localisation des stations
Calculer le centre centerde la carte des données compteurs. Il s’obtient en agrègeant l’ensemble des géométries, calculant le centroid et récupèrant la valeur sous forme de liste. Avec une logique similaire, calculez les bornes du sud-ouest sw et du nord-est ne de la carte.
Représenter la localisation des stations en utilisant un zoom optimal.
# Afficher la cartem
Exercice 7: Représenter les stations
Faire la même carte, avec des ronds proportionnels au nombre de comptages :
Pour le rayon de chaque cercle, vous pouvez appliquer la règle 500*x/max(x) (règle au doigt mouillé)
Vous pouvez réduire la taille des bordures de cercle avec l’option weight = 1 et fixer la couleur avec color = 'grey'
(Optionnel) Colorer en rouge les 10 plus grosses stations. L’opacité étant, par défaut, un peu faible, le paramètre fill_opacity = 0.4 améliore le rendu.
(Optionnel) Afficher, en supplément du nom du compteur lorsqu’on clique, la valeur du comptage en revenant à la ligne
La carte obtenue doit ressembler à la suivante :
# Afficher la cartem
Exercices supplémentaires
Densité de population dans la petite couronne parisienne
Pour cet exercice, le package cartiflette
va être pratique pour récupérer un fonds de carte mélangeant arrondissements
parisiens et communes dans les autres villes.
Nous allons privilégier une carte à ronds proportionnels (bubble map)
aux cartes chorolèpthes qui trompent
l’oeil. Les instructions d’installation du package topojson
sont disponibles dans la partie manipulation
Exercice: bubble map de densité des populations
Récupérer le fond de carte des départements 75, 92, 93 et 94
avec cartiflette. Pour cela, utiliser download_vectorfile_url_all
depuis cartiflette.s3 en fixant l’option level à COMMUNE_ARRONDISSEMENT.
Nommer cet objet df.
Afin que les calculs ultérieurs de surface ne soient pas faussés,
assurez-vous que les données sont en Lambert 93 en reprojetant
nos contours (code EPSG: 2154).
Créer un objet departements avec dissolve pour également disposer
d’un fond de carte des départements
Créer une variable surface et utilisant la méthode area. L’unité
doit être le km², il faut donc diviser par \(10^6\)
Créer une variable densite
Utiliser pd.cut avec les seuils 5000, 15000 et 30000 personnes
par km². Vous pouvez utiliser l’option label pour dénommer les tranches
Créer un GeoDataFrame de points en utilisant la méthode centroid. Celui-ci
nous servira à localiser le centre de nos ronds.
Représenter la densité communale sous forme de carte avec ronds proportionnels.
Vous pouvez utiliser la variable créée à la question 5 pour les couleurs.
La carte obtenue devrait ressembler à celle-ci:
/opt/mamba/lib/python3.9/site-packages/geopandas/plotting.py:730: FutureWarning:
is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
@book{galiana2023,
author = {Galiana, Lino},
title = {Python Pour La Data Science},
date = {2023},
url = {https://pythonds.linogaliana.fr/},
doi = {10.5281/zenodo.8229676},
langid = {en}
}