Données spatiales: découverte de geopandas

Download nbviewer Onyxia
Binder Open In Colab githubdev

</p>

Dans ce tutoriel, nous allons utiliser:

La représentation des données, notamment la cartographie, est présentée plus amplement dans la partie visualiser. Quelques méthodes pour faire rapidement des cartes seront présentées ici mais ce n’est pas l’objet de ce chapitre.

Ce tutoriel s’inspire beaucoup d’un autre tutoriel que j’ai fait pour R disponible dans la documentation utilitr. Il peut servir de pendant à celui-ci pour l’utilisateur de R.

Quelques installations préalables sont nécessaires :

!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

Pour être en mesure d’exécuter ce tutoriel, les imports suivants seront utiles.

In [2]:
import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt

Données spatiales: quelle différence avec des données traditionnelles ?

Le terme “données spatiales” désigne les données qui portent sur les caractéristiques géographiques des objets (localisation, contours, liens). Les caractéristiques géographiques des objets sont décrites à l’aide d’un système de coordonnées qui permettent une représentation dans un espace euclidien ($(x,y)$). Le passage de l’espace réel (la Terre, qui est une sphère) à l’espace plan se fait grâce à un système de projection. Voici quelques exemples de données spatiales:

  • Une table décrivant des bâtiments, avec les coordonnées géographiques de chaque bâtiment;
  • Le découpage communal du territoire, avec le contour du territoire de chaque commune;
  • Les routes terrestres, avec les coordonnées décrivant leur parcours.

Les données spatiales rassemblent classiquement deux types de données :

  1. des données géographiques (ou géométries): objets géométriques tels que des points, des vecteurs, des polygones, ou des maillages (raster). Exemple: la forme de chaque chaque commune, les coordonnées d’un bâtiment;
  2. des données attributaires (ou attributs): des mesures et des caractéristiques associés aux objets géométriques. Exemple: la population de chaque commune, le nombre de fenêtres et le nombre d’étages d’un bâtiment.

Les données spatiales sont fréquemment traitées à l’aide d’un système d’information géographique (SIG), c’est-à-dire un système d’information capable de stocker, d’organiser et de présenter des données alphanumériques spatialement référencées par des coordonnées dans un système de référence (CRS). Python dispose de fonctionnalités lui permettant de réaliser les mêmes tâches qu’un SIG (traitement de données spatiales, représentations cartographiques).

Les systèmes de projection font l’objet de standards internationaux et sont souvent désignés par des codes dits codes EPSG. Ce site est un bon aide-mémoire. Les plus fréquents, pour les utilisateurs français, sont les suivants (plus d’infos ici):

  • 2154: système de projection Lambert 93. Il s’agit du système de projection officiel. La plupart des données diffusées par l’administration pour la métropole sont disponibles dans ce système de projection.
  • 27572: Lambert II étendu. Il s’agit de l’ancien système de projection officiel. Les données spatiales anciennes peuvent être dans ce format.
  • 4326: WGS 84 ou système de pseudo-Mercator. Attention, ce n’est en réalité pas un système de projection mais un système de coordonnées (longitude / latitude) qui permet simplement un repérage angulaire sur l’ellipsoïde. Il est utilisé pour les données GPS.

De pandas à geopandas

Le package geopandas est une boîte à outils conçue pour faciliter la manipulation de données spatiales. La grande force de geopandas est qu’il permet de manipuler des données spatiales comme s’il s’agissait de données traditionnelles, car il repose sur le standard ISO 19125 simple feature access défini conjointement par l’Open Geospatial Consortium (OGC) et l’International Organization for Standardization (ISO).

Par rapport à un DataFrame standard, un objet geopandas comporte une colonne supplémentaire: geometry. Elle stocke les coordonnées des objets géographiques (ou ensemble de coordonnées s’agissant de contours). Un objet geopandas hérite des propriétés d’un DataFrame pandas mais propose des méthodes adaptées au traitement des données spatiales.

Ainsi, grâce à geopandas, on pourra effectuer des manipulations sur les attributs des données comme avec pandas mais on pourra également faire des manipulations sur la dimension spatiale des données. En particulier,

  • Calculer des distances et des surfaces;
  • Agréger rapidement des zonages (regrouper les communes en département par exemple);
  • Trouver dans quelle commune se trouve un bâtiment à partir de ses coordonnées géographiques;
  • Recalculer des coordonnées dans un autre système de projection.
  • Faire une carte, rapidement et simplement

Par rapport à un logiciel spécialisé comme QGIS, python permettra d’automatiser le traitement et la représentation des données. D’ailleurs, QGIS utilise lui-même python

Importer des données spatiales

Les données spatiales sont plus riches que les données traditionnelles car elles incluent, habituellement, des éléments supplémentaires pour placer dans un espace cartésien les objets. Cette dimension supplémentaire peut être simple (un point comporte deux informations supplémentaire: $x$ et $y$) ou assez complexe (polygones, lignes avec direction, etc.)

Les formats les plus communs de données spatiales sont les suivants:

  • shapefile (.shp): format (propriétaire) le plus commun de données géographiques. La table de données (attributs) est stockée dans un fichier séparé des données spatiales. En faisant geopandas.read_file("monfichier.shp"), le package fait lui-même le lien entre les observations et leur représentation spatiale ;
  • geojson (.json) : ce format, non préconisé par l’OGC est largement utilisé pour le développement web comme dans la librairie Leaflet. La dimension spatiale est stockée dans le même fichier que les attributs. Ces fichiers sont généralement beaucoup plus légers que les shapefiles mais possèdent des limites s’agissant de gros jeux de données.
  • geopackage (.gpkg) : ce (relativement) nouveau format libre en un seul fichier également (lui recommandé par l’OGC) vise progressivement à se substituer au shapefile. Il est par exemple le format par défaut dans QGIS.

Cette page compare plus en détail ces trois types de formats de données géographiques.

L’aide de geopandas propose des bouts de code en fonction des différentes situations dans lesquelles on se trouve.

Exemple: récupérer les découpages territoriaux

L’un des fonds de carte les plus fréquents qu’on utilise est celui des limites administratives. Elles peuvent être trouvées sur le site de l’IGN ou récupérées sur data.gouv (construites par openstreetmap). Il est également possible d’utiliser l’une des API de l’IGN mais ces dernières ne sont pas encore très documentées pour des utilisateurs de Python.

Nous allons utiliser le package pynsee qui propose un module dédié à la récupération de fonds de carte officiels pour valoriser des données d’open data. Une précédente version de ce cours utilisait les données de data.gouv, elle est reportée en annexe

L’inconvénient est que les arrondissements parisiens ne sont pas présents dans le fichier proposé sur data.gouv. Il faut donc utiliser une source complémentaire, issue de l’opendata de la Mairie de Paris.

In [3]:
from pynsee.geodata import get_geodata
shp_communes = gpd.GeoDataFrame(
  get_geodata('ADMINEXPRESS-COG-CARTO.LATEST:commune')
)
shp_communes = shp_communes.rename({"insee_com": 'insee'}, axis = 'columns')
shp_communes = shp_communes.set_crs(3857)
In [4]:
shp_communes.head()
id nom nom_m insee statut population insee_can insee_arr insee_dep insee_reg siren_epci geometry bbox crs
0 COMMUNE_0000001269976025 Le Diamant LE DIAMANT 97206 Commune simple 5511 NR 3 972 02 249720053 MULTIPOLYGON (((-6794714.852 1625025.923, -679... [-6796892.81639186, 1624969.98776564, -6788212... EPSG:3857
1 COMMUNE_0000001269976026 Les Anses-d'Arlet LES ANSES-D'ARLET 97202 Commune simple 3494 NR 3 972 02 249720053 MULTIPOLYGON (((-6796770.257 1632364.571, -679... [-6801724.22854182, 1626341.46982216, -6795643... EPSG:3857
2 COMMUNE_0000001269976019 Rivière-Salée RIVIERE-SALEE 97221 Commune simple 11874 NR 3 972 02 249720053 MULTIPOLYGON (((-6788283.345 1637546.648, -678... [-6790993.527111, 1631234.68464338, -6781597.8... EPSG:3857
3 COMMUNE_0000001269975998 Saint-Pierre SAINT-PIERRE 97225 Sous-préfecture 4121 NR 4 972 02 200041788 MULTIPOLYGON (((-6814293.278 1664074.082, -681... [-6814293.27835075, 1657968.05307498, -6805323... EPSG:3857
4 COMMUNE_0000001269975994 Sainte-Marie SAINTE-MARIE 97228 Commune simple 15487 NR 2 972 02 200041788 MULTIPOLYGON (((-6786330.097 1662972.643, -678... [-6796534.21176401, 1658613.16621139, -6786313... EPSG:3857

On reconnaît la structure d’un DataFrame pandas. A cette structure s’ajoute une colonne geometry qui enregistre la position des limites des polygones de chaque observation.

Comme vu précédemment, le système de projection est un élément important. Il permet à python d’interpréter les valeurs des points (deux dimensions) en position sur la terre, qui n’est pas un espace plan. Ici, les données ne sont pas projetées puisqu’elles sont dans le système WSG84 (epsg: 4326) ce qui permet de facilement ajouter un fonds de carte openstreetmap ou stamen pour rendre une représentation graphique plus esthétique. En toute rigueur, pour faire une carte statique d’un pays en particulier, il faudrait reprojeter les données dans un système de projection adapté à la zone géographique étudiée (par exemple le Lambert-93 pour la France métropolitaine).

On peut ainsi représenter Paris pour se donner une idée de la nature du shapefile utilisé :

In [5]:
paris = shp_communes.loc[shp_communes['insee'].str.startswith("75")]

fig,ax = plt.subplots(figsize=(10, 10))
paris.plot(ax = ax, alpha=0.5, edgecolor='blue')
ctx.add_basemap(ax, crs = paris.crs.to_string())
ax.set_axis_off()
ax
<AxesSubplot:>

On voit ainsi que les données pour Paris ne comportent pas d’arrondissement, ce qui est limitant pour une analyse focalisée sur Paris. On va donc les récupérer sur le site d’open data de la ville de Paris et les substituer à Paris :

En refaisant la carte ci-dessus, on obtient bien, cette fois, le résultat espéré

In [7]:
paris = communes.loc[communes.insee.str.startswith("75")]

fig,ax = plt.subplots(figsize=(10, 10))

paris.plot(ax = ax, alpha=0.5, edgecolor='k')
ctx.add_basemap(ax, crs = paris.crs.to_string())
ax.set_axis_off()
ax
<AxesSubplot:>