Données spatiales : découverte de geopandas

Les données géolocalisées se sont multipliées depuis quelques années, qu’il s’agisse de données open-data ou de traces numériques géolocalisées de type big-data. Pour les données spatiales, le package GeoPandas étend les fonctionalités de l’écosystème Pandas afin de permettre de manipuler des données géographiques complexes de manière simple.

Tutoriel
Manipulation
Author

Lino Galiana

Published

2024-04-27

Dans ce tutoriel, nous allons utiliser les données suivantes :

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 l’objet de ce chapitre porte davantage sur la manipulation des données géographiques.

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
!pip install topojson

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

import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt

1 Données spatiales

1.1 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 en 3 dimensions (longitude, latitude, altitude).

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 commune, les coordonnées d’un bâtiment;
  2. des données attributaires (ou attributs) : des mesures et des caractéristiques associées 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).

1.2 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…

1.3 Résumé

En résumé, un objet GeoPandas comporte les éléments suivantes :

  1. Les attributs. Ce sont les valeurs associées à chaque niveau géographique. Il s’agit de la dimension tabulaire usuelle, dont le traitement est similaire à celui d’un objet Pandas classique.
  2. Les géométries. Ce sont les valeurs numériques interprétées pour représenter la dimension géographique. Elles permettent de représenter dans un certain référentiel (le système de référence) la dimension géographique.
  3. Le système de référence. Il s’agit du système permettant de transformer les positions sur le globe (3 dimensions avec une boule asymétrique) en un plan en deux dimensions. Il en existe une multitude, identifiables à partir d’un code EPSG (4326, 2154…). Leur manipulation est facilitée par Geopandas qui s’appuie sur Shapely, de la même manière que Pandas s’appuie sur Numpy ou Arrow.

2 Le système de projection cartographique

2.1 Principe

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.).

L’analyse cartographique emprunte dès lors à la géométrie des concepts pour représenter des objets dans l’espace. Les projections sont au coeur de la gestion des données spatiales. Ces dernières consistent à transformer une position dans l’espace terrestre à une position sur un plan. Il s’agit donc d’une opération de projection d’un espace tri-dimensionnel dans un espace à deux dimensions. Ce post propose de riches éléments sur le sujet, notamment l’image suivante qui montre bien le principe d’une projection :

Les différents types de projection

Les différents types de projection

Cette opération n’est pas neutre. L’une des conséquences du théorème remarquable de Gauss est que la surface de la Terre ne peut être cartographiée sans distortion. Une projection ne peut simultanément conserver intactes les distances et les angles (i.e. les positions). Il n’existe ainsi pas de projection universellement meilleure, ce qui ouvre la porte à la coexistence de nombreuses projections différentes, pensées pour des tâches différentes. Un mauvais système de représentation fausse l’appréciation visuelle mais peut aussi entraîner des erreurs dans les calculs sur la dimension spatiale.

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 ou encore Web Mercator. 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. Il s’agit du système le plus usuel, notamment quand on travaille avec des fonds de carte web.

Comme évoqué plus haut, l’une des projections les plus connues est la projection Web Mercator dite WGS84 (code EPSG 4326). Il s’agit d’une projection conservant intacte les angles, ce qui implique qu’elle altère les distances. Celle-ci a en effet été pensée, à l’origine, pour représenter l’hémisphère Nord. Plus on s’éloigne de celui-ci, plus les distances sont distordues. Cela amène à des distorsions bien connues (le Groenland hypertrophié, l’Afrique de taille réduite, l’Antarctique démesuré…). En revanche, la projection Mercator conserve intacte les positions. C’est cette propriété qui explique son utilisation dans les systèmes GPS et ainsi dans les fonds de carte de navigation du type Google Maps.

Exemple de reprojection de pays depuis le site thetruesize.com

Exemple de reprojection de pays depuis le site thetruesize.com

Observez les variations significatives de proportions pour certains pays selon les projections choisies:

Pour aller plus loin, la carte interactive suivante, construite par Nicolas Lambert, issue de ce notebook Observable, illustre l’effet déformant de la projection Mercator, et de quelques-unes autres, sur notre perception de la taille des pays.

Voir la carte interactive

Il existe en fait de nombreuses représentations possibles du monde, plus ou moins alambiquées. Les projections sont très nombreuses et certaines peuvent avoir une forme suprenante. Par exemple, la projection de Spillhaus propose de centrer la vue sur les océans et non une terre. C’est pour cette raison qu’on parle parfois de monde tel que vu par les poissons à son propos.

3 Importer des données spatiales

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.
  • 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.
  • geojson (.json) : ce format, non préconisé par l’OGC, est largement utilisé pour le développement web comme dans la librairie leaflet.js. 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.
  • topojson (.json) : une variante du geojson qui se développe progressivement pour assister les visualisations web. Au lieu de stocker l’ensemble des points permettant de représenter une géométrie, seuls les arcs sont conservés. Cela allège substantiellement le poids du fichier et permet, avec une librairie adaptée, de reconstruire l’ensemble des contours géographiques.

Cette page compare plus en détail les principes 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.

3.1 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 des communes. Celui-ci peut être récupéré de plusieurs manières. En premier lieu, pour récupérer le fond de carte officiel, produit par l’IGN, sous le nom d’AdminExpress[^1], il est possible de se rendre sur le site de l’IGN et de le télécharger. 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. Le package cartiflette, issu d’un projet interministériel, propose une récupération facilitée de fonds de carte officiels de l’IGN. Ce projet vise à faciliter la récupération des sources officielles, notamment celles de l’IGN, et leur association à des jeux de données géographiques.

Ici, nous sommes intéressés par les contours des communes de la petite couronne. On pourrait désirer récupérer l’ensemble de la région Ile-de-France mais nous allons nous contenter de l’analyse de Paris intra-muros et des départements limitrophes. C’est l’un des avantage de cartiflette que de faciliter la récupération de fonds de carte sur un ensemble de département. Cela évite la récupération d’un fond de carte très volumineux (plus de 500Mo) pour une analyse restreinte (quelques départements). Un autre avantage de cartiflette est de faciliter la récupération de fonds de carte consolidés comme celui dont on a besoin ici : arrondissements dans Paris, communes ailleurs. Comme cela est expliqué dans un encadré à part, il s’agirait d’une opération pénible à mettre en oeuvre sans cartiflette.

Les contours de cet espace peuvent être récupérés de la manière suivante :

from cartiflette import carti_download

shp_communes = carti_download(
    crs=4326,
    values=["75", "92", "93", "94"],
    borders="COMMUNE_ARRONDISSEMENT",
    vectorfile_format="geojson",
    filter_by="DEPARTEMENT",
    source="EXPRESS-COG-CARTO-TERRITOIRE",
    year=2022,
)
shp_communes.head(3)
INSEE_DEP INSEE_REG ID NOM INSEE_COM STATUT POPULATION INSEE_COG ARR CV ... TAAV2017 TDAAV2017 CATEAAV2020 BV2012 LIBELLE_DEPARTEMENT LIBELLE_REGION PAYS SOURCE geometry AREA
0 75 11 ARR_MUNI0000000009736045 Paris 3e Arrondissement 75056 Arrondissement municipal 34025 75103 751 75ZZ ... 5 50 11 75056 Paris Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((2.35016 48.86199, 2.35350 48.86121, ... NaN
1 75 11 ARR_MUNI0000000009736046 Paris 2e Arrondissement 75056 Arrondissement municipal 21595 75102 751 75ZZ ... 5 50 11 75056 Paris Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((2.34792 48.87069, 2.34004 48.87197, ... NaN
2 75 11 ARR_MUNI0000000009736545 Paris 4e Arrondissement 75056 Arrondissement municipal 29131 75104 751 75ZZ ... 5 50 11 75056 Paris Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((2.36849 48.85581, 2.36828 48.85579, ... NaN

3 rows × 27 columns

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.

shp_communes.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Ici, les données sont dans le système WGS84 (code EPSG 4326). Ce n’est pas le Lambert-93 comme on pourrait s’y attendre, ce dernier étant le système légal de projection pour la France métropolitaine.

Pour s’assurer qu’on a bien récupéré les contours voulus, on peut représenter graphiquement les contours grâce à la méthode plot sur laquelle nous reviendrons :

ax = shp_communes.boundary.plot()
ax.set_axis_off()

4 Opérations sur les attributs et les géométries

4.1 Import des données vélib

Souvent, le découpage communal ne sert qu’en fond de cartes, pour donner des repères. En complément de celui-ci, on peut désirer exploiter un autre jeu de données. On va partir des données de localisation des stations velib, disponibles sur le site d’open data de la ville de Paris et requêtables directement par l’url https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr

velib_data = "https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr"
stations = gpd.read_file(velib_data)
stations.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Les données sont dans le système de projection WGS84 qui est celui du système GPS. Celui-ci s’intègre bien avec les fonds de carte OpenStreetMap ou Google Maps. En toute rigueur, si on désire effectuer certains calculs géométriques (mesurer des surfaces…), il est nécessaire de re-projeter les données dans un système qui préserve la géométrie (c’est le cas du Lambert 93).

Pour avoir une intuition de la localisation des stations, et notamment de la densité hétérogène de celles-ci, on peut afficher les données sur la carte des communes de la petite couronne. Il s’agit donc d’enrichir la carte précédente d’une couche supplémentaire, à savoir la localisation des stations. Au passage, on va utiliser un fond de carte plus esthétique:

fig, ax = plt.subplots(figsize=(10, 10))
stations.sample(200).to_crs(3857).plot(ax=ax, color="red", alpha=0.4, zorder=2)
shp_communes.to_crs(3857).plot(
    ax=ax, zorder=1, edgecolor="black", facecolor="none", color=None
)
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)
ax.set_axis_off()

Découvrez ci-dessous par étape les différentes lignes de commandes permettant d’afficher cette carte complète, étape par étape :

1️⃣ Afficher le nuage de points de 200 stations vélibs prises au hasard

fig, ax = plt.subplots(figsize=(10, 10))
stations.sample(200).to_crs(3857).plot(ax=ax, color="red", alpha=0.4, zorder=2)

2️⃣ Ajouter à cette couche, en-dessous, les contours des communes

fig, ax = plt.subplots(figsize=(10, 10))
stations.sample(200).to_crs(3857).plot(ax=ax, color="red", alpha=0.4, zorder=2)
shp_communes.to_crs(3857).plot(
    ax=ax, zorder=1, edgecolor="black", facecolor="none", color=None
)

3️⃣ Ajouter un fond de carte de type open street map grâce au package contextily

fig, ax = plt.subplots(figsize=(10, 10))
stations.sample(200).to_crs(3857).plot(ax=ax, color="red", alpha=0.4, zorder=2)
shp_communes.to_crs(3857).plot(
    ax=ax, zorder=1, edgecolor="black", facecolor="none", color=None
)
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)

4️⃣ Il ne reste plus qu’à retirer l’axe des coordonnées, qui n’est pas très esthétique.

fig, ax = plt.subplots(figsize=(10, 10))
stations.sample(200).to_crs(3857).plot(ax=ax, color="red", alpha=0.4, zorder=2)
shp_communes.to_crs(3857).plot(
    ax=ax, zorder=1, edgecolor="black", facecolor="none", color=None
)
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)
ax.set_axis_off()
ax

In fine, on obtient la carte désirée.

4.2 Opérations sur les attributs

Toutes les opérations possibles sur un objet Pandas le sont également sur un objet GeoPandas. Pour manipuler les données, et non la géométrie, on parlera d’opérations sur les attributs.

Par exemple, si on désire connaître quelques statistiques sur la taille des stations, l’approche est identique à si on avait un objet Pandas classique :

stations.describe()
capacity
count 1468.000000
mean 31.168256
std 12.291324
min 0.000000
25% 23.000000
50% 29.000000
75% 37.000000
max 76.000000

Pour classer les départements de la petite couronne, du plus grand au plus petit, procédons en deux étapes:

  1. Récupérons le contour des communes grâce à cartiflette. Notons qu’on pourrait récupérer directement les contours départementaux mais pour l’exercice, nous allons le créer nous-mêmes comme agrégation des contours communaux (voir plus bas ainsi que ce notebook Observable pour la méthode plus légère qui utilise pleinement les fonctionnalités de cartiflette).

  2. Calculons la surface totale de ce territoire (méthode area sur un objet GeoPandas.GeoDataFrame ramenée en km², attention néamoins au système de projection comme cela est expliqué plus bas)

shp_communes["surface"] = shp_communes.area.div(10**6)

Les plus grands départements s’obtiennent par une agrégation des surfaces communales :

shp_communes.groupby("INSEE_DEP").sum(numeric_only=True).sort_values(
    "surface", ascending=False
)
INSEE_REG POPULATION ZE2020 TUU2017 TDUU2017 TAAV2017 TDAAV2017 CATEAAV2020 surface
INSEE_DEP
94 517 1407124 52123 376 3760 235 2350 564 244.816689
93 440 1644903 44384 320 3200 200 2000 480 236.867414
92 396 1624357 39924 288 2880 180 1800 432 175.571633
75 220 2165423 22180 160 1600 100 1000 220 105.424231

Si on veut directement les plus grandes communes de la petite couronne parisienne :

shp_communes.sort_values("surface", ascending=False).head(10)
INSEE_DEP INSEE_REG ID NOM INSEE_COM STATUT POPULATION AREA ARR CV ... TDAAV2017 CATEAAV2020 BV2012 LIBELLE_DEPARTEMENT LIBELLE_REGION PAYS SOURCE geometry INSEE_COG surface
60 93 11 COMMUNE_0000000009735015 Tremblay-en-France 93073 Commune simple 36461 metropole 932 9320 ... 50 12 75056 Seine-Saint-Denis Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((663432.700 6875170.500, 663529.700 6... NaN 22.680549
12 75 11 ARR_MUNI0000000009736553 Paris 16e Arrondissement 75056 Arrondissement municipal 165523 NaN 751 75ZZ ... 50 11 75056 Paris Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((647190.200 6864524.400, 647155.000 6... 75116 16.412426
19 75 11 ARR_MUNI0000000009736532 Paris 12e Arrondissement 75056 Arrondissement municipal 139297 NaN 751 75ZZ ... 50 11 75056 Paris Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((655221.200 6858576.500, 655536.100 6... 75112 16.379110
44 93 11 COMMUNE_0000000009735500 Aulnay-sous-Bois 93005 Commune simple 86969 metropole 932 9302 ... 50 12 75056 Seine-Saint-Denis Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((660415.900 6872923.300, 660459.000 6... NaN 16.166721
32 92 11 COMMUNE_0000000009736056 Rueil-Malmaison 92063 Commune simple 78317 metropole 922 9222 ... 50 12 75056 Hauts-de-Seine Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((637679.300 6863761.400, 637603.700 6... NaN 14.541191
74 93 11 COMMUNE_0000000009736517 Noisy-le-Grand 93051 Commune simple 67871 metropole 932 9314 ... 50 12 75056 Seine-Saint-Denis Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((664453.800 6861358.900, 664703.700 6... NaN 13.144471
37 93 11 COMMUNE_0000000009735515 Saint-Denis 93066 Sous-préfecture 112852 metropole 933 9399 ... 50 12 75056 Seine-Saint-Denis Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((652098.200 6872022.600, 652002.600 6... NaN 12.371854
5 92 11 COMMUNE_0000000009736052 Nanterre 92050 Préfecture 96277 metropole 922 9299 ... 50 12 75056 Hauts-de-Seine Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((643490.700 6867612.300, 643408.400 6... NaN 12.227927
114 94 11 COMMUNE_0000000009737009 Vitry-sur-Seine 94081 Commune simple 95510 metropole 943 9499 ... 50 12 75056 Val-de-Marne Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((653537.000 6852608.400, 653581.100 6... NaN 11.664610
21 92 11 COMMUNE_0000000009735517 Gennevilliers 92036 Commune simple 48530 metropole 922 9214 ... 50 12 75056 Hauts-de-Seine Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((648070.700 6872566.900, 647847.400 6... NaN 11.633935

10 rows × 28 columns

Lors des étapes d’agrégation, groupby ne conserve pas les géométries. Autrement dit, si on effectue, par exemple, une somme en fonction d’une variable de groupe avec le combo groupby(...).sum(...) , on perd la dimension géographique.

Il est néanmoins possible d’aggréger à la fois les géométries et les attribus avec la méthode dissolve:

fig, ax = plt.subplots(figsize=(10, 10))
shp_communes.dissolve(by="INSEE_DEP", aggfunc="sum").plot(ax=ax, column="surface")
ax.set_axis_off()
ax

Pour produire l’équivalent de cette carte à un niveau France entière, il est néanmoins plus simple de directement récupérer les fonds officiels des départements plutôt que d’agréger les contours des communes:

dep = carti_download(
    values=["France"],
    crs=4326,
    borders="DEPARTEMENT",
    vectorfile_format="geojson",
    filter_by="FRANCE_ENTIERE",
    source="EXPRESS-COG-CARTO-TERRITOIRE",
    year=2022,
)
dep = dep.loc[dep["INSEE_DEP"].str.len() == 2]

dep["area"] = dep.to_crs(2154).area

Avant de calculer les surfaces des départements, pour éviter les déformations liées au système Mercator, nous faisons une reprojection des données à la volée. Plus de détails par la suite.

dep.sort_values("area", ascending=False).head(3)
INSEE_DEP PAYS LIBELLE_DEPARTEMENT POPULATION SOURCE geometry area
55 33 France Gironde 1623749 IGN:EXPRESS-COG-CARTO-TERRITOIRE MULTIPOLYGON (((-0.71884 45.32753, -0.71845 45... 1.007271e+10
30 40 France Landes 413690 IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((-0.24252 43.58483, -0.24275 43.58706... 9.354177e+09
68 24 France Dordogne 413223 IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((1.44824 45.01941, 1.44400 45.02011, ... 9.210808e+09 
ax = dep.plot(column="area")
ax.set_axis_off()

4.3 Opérations sur les géométries

Outre la représentation graphique simplifiée, sur laquelle nous reviendrons ultérieurement, l’intérêt principal d’utiliser GeoPandas est l’existence de méthodes efficaces pour manipuler la dimension spatiale. Un certain nombre proviennent du package Shapely.

Comme indiqué ci-dessus, nous reprojetons les données dans le système Lambert 93 qui ne fausse pas les calculs de distance et d’aires.

communes = shp_communes.to_crs(2154)
stations = stations.to_crs(2154)

Par exemple, on peut recalculer la taille d’une commune ou d’arrondissement avec la méthode area (et diviser par \(10^6\) pour avoir des \(km^2\) au lieu des \(m^2\)):

communes["superficie"] = communes.area.div(10**6)
communes.head(3)
INSEE_DEP INSEE_REG ID NOM INSEE_COM STATUT POPULATION AREA ARR CV ... CATEAAV2020 BV2012 LIBELLE_DEPARTEMENT LIBELLE_REGION PAYS SOURCE geometry INSEE_COG surface superficie
1 92 11 COMMUNE_0000000009736037 Levallois-Perret 92044 Commune simple 66082 metropole 922 9216 ... 12 75056 Hauts-de-Seine Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((647761.400 6867306.900, 647531.300 6... NaN 2.412607 2.412607
2 92 11 COMMUNE_0000000009736055 Bois-Colombes 92009 Commune simple 28841 metropole 922 9211 ... 12 75056 Hauts-de-Seine Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((646224.700 6867615.800, 646122.500 6... NaN 1.931211 1.931211
3 92 11 COMMUNE_0000000009736538 Malakoff 92046 Commune simple 30950 metropole 921 9218 ... 12 75056 Hauts-de-Seine Île-de-France France IGN:EXPRESS-COG-CARTO-TERRITOIRE POLYGON ((646995.300 6857373.400, 646723.300 6... NaN 2.073708 2.073708

3 rows × 29 columns

Une méthode qu’on utilise régulièrement est centroid qui, comme son nom l’indique, recherche le centroïde de chaque polygone et transforme ainsi des données surfaciques en données ponctuelles. Par exemple, pour représenter approximativement les centres des villages de la Haute-Garonne (31), après avoir téléchargé le fonds de carte adapté, fera

communes_31 = carti_download(
    values=["31"],
    crs=4326,
    borders="COMMUNE",
    vectorfile_format="geojson",
    filter_by="DEPARTEMENT",
    source="EXPRESS-COG-CARTO-TERRITOIRE",
    year=2022,
)

# on reprojete en 3857 pour le fond de carte
communes_31 = communes_31.to_crs(3857)

# on calcule le centroide
dep_31 = communes_31.copy()
communes_31["geometry"] = communes_31["geometry"].centroid

ax = communes_31.plot(figsize=(10, 10), color="red", alpha=0.4, zorder=2)
dep_31.to_crs(3857).plot(
    ax=ax, zorder=1, edgecolor="black", facecolor="none", color=None
)
# ctx.add_basemap(ax, source = ctx.providers.Stamen.Toner)
ax.set_axis_off()
ax

5 Gérer le système de projection

Précédemment, nous avons appliqué une méthode to_crs pour reprojeter les données dans un système de projection différent de celui du fichier d’origine :

communes = communes.to_crs(2154)
stations = stations.to_crs(2154)

Concernant la gestion des projections avec GeoPandas, la documentation officielle est très bien faite. Elle fournit notamment l’avertissement suivant qu’il est bon d’avoir en tête :

Be aware that most of the time you don’t have to set a projection. Data loaded from a reputable source (using the geopandas.read_file() command) should always include projection information. You can see an objects current CRS through the GeoSeries.crs attribute.

From time to time, however, you may get data that does not include a projection. In this situation, you have to set the CRS so geopandas knows how to interpret the coordinates.

Image empruntée à XKCD https://xkcd.com/2256/ qu’on peut également trouver sur https://blog.chrislansdown.com/2020/01/17/a-great-map-projection-joke/

Image empruntée à XKCD https://xkcd.com/2256/ qu’on peut également trouver sur https://blog.chrislansdown.com/2020/01/17/a-great-map-projection-joke/

Pour déterminer le système de projection d’une base de données, on peut vérifier l’attribut crs :

communes.crs
<Projected CRS: EPSG:2154>
Name: RGF93 v1 / Lambert-93
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: France - onshore and offshore, mainland and Corsica (France métropolitaine including Corsica).
- bounds: (-9.86, 41.15, 10.38, 51.56)
Coordinate Operation:
- name: Lambert-93
- method: Lambert Conic Conformal (2SP)
Datum: Reseau Geodesique Francais 1993 v1
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Les deux principales méthodes pour définir le système de projection utilisé sont :

  • df.set_crs : cette commande sert à préciser quel est le système de projection utilisé, c’est-à-dire comment les coordonnées (x,y) sont reliées à la surface terrestre. Cette commande ne doit pas être utilisée pour transformer le système de coordonnées, seulement pour le définir.
  • df.to_crs : cette commande sert à projeter les points d’une géométrie dans une autre, c’est-à-dire à recalculer les coordonnées selon un autre système de projection.

Dans le cas particulier de production de carte avec un fond OpenStreetMaps ou une carte dynamique leaflet, il est nécessaire de dé-projeter les données (par exemple à partir du Lambert-93) pour atterrir dans le système non-projeté WGS 84 (code EPSG 4326). Ce site dédié aux projections géographiques peut être utile pour retrouver le système de projection d’un fichier où il n’est pas indiqué.

La définition du système de projection se fait de la manière suivante (:warning: avant de le faire, se souvenir de l’avertissement !) :

communes = communes.set_crs(2154)

Alors que la reprojection (projection Albers : 5070) s’obtient de la manière suivante :

shp_region = carti_download(
    values=["France"],
    crs=4326,
    borders="REGION",
    vectorfile_format="geojson",
    filter_by="FRANCE_ENTIERE",
    source="EXPRESS-COG-CARTO-TERRITOIRE",
    year=2022,
)
shp_region = shp_region.loc[shp_region["INSEE_REG"] > 10]
fig, ax = plt.subplots(figsize=(10, 10))
shp_region.to_crs(5070).plot(ax=ax)
ax

<Figure size 672x480 with 0 Axes>

On le voit, cela modifie totalement la représentation de l’objet dans l’espace. Clairement, cette projection n’est pas adaptée aux longitudes et latitudes françaises. C’est normal, il s’agit d’une projection adaptée au continent nord-américain (et encore, pas dans son ensemble !).

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))

fig, ax = plt.subplots(figsize=(10, 10))
world[world.continent == "North America"].to_crs(5070).plot(
    alpha=0.2, edgecolor="k", ax=ax
)
ax

<Figure size 672x480 with 0 Axes>

6 Joindre des données

6.1 Joindre des données sur des attributs

Ce type de jointure se fait entre un objet géographique et un deuxième objet, géographique ou non. A l’exception de la question des géométries, il n’y a pas de différence par rapport à Pandas.

La seule différence avec Pandas est dans la dimension géographique. Si on désire conserver la dimension géographique, il faut faire attention à faire :

geopandas_object.merge(pandas_object)

Si on utilise deux objets géographiques mais ne désire conserver qu’une seule dimension géographique1, on fera

geopandas_object1.merge(geopandas_object2)

Seule la géométrie de l’objet de gauche sera conservée, même si on fait un right join.

6.2 Prolongation possible : joindre des données sur dimension géographique

Le chapitre suivant permettra de mettre en oeuvre des jointures géographiques.

Informations additionnelles

environment files have been tested on.

Latest built version: 2024-04-27

Python version used:

'3.11.6 | packaged by conda-forge | (main, Oct  3 2023, 10:40:35) [GCC 12.3.0]'
Package Version
affine 2.4.0
aiobotocore 2.12.2
aiohttp 3.9.3
aioitertools 0.11.0
aiosignal 1.3.1
alembic 1.13.1
aniso8601 9.0.1
annotated-types 0.6.0
appdirs 1.4.4
archspec 0.2.3
astroid 3.1.0
asttokens 2.4.1
attrs 23.2.0
Babel 2.14.0
bcrypt 4.1.2
beautifulsoup4 4.12.3
black 24.4.2
blinker 1.7.0
blis 0.7.11
bokeh 3.4.0
boltons 23.1.1
boto3 1.34.51
botocore 1.34.51
branca 0.7.1
Brotli 1.1.0
cachetools 5.3.3
cartiflette 0.0.2
Cartopy 0.23.0
catalogue 2.0.10
cattrs 23.2.3
certifi 2024.2.2
cffi 1.16.0
charset-normalizer 3.3.2
click 8.1.7
click-plugins 1.1.1
cligj 0.7.2
cloudpathlib 0.16.0
cloudpickle 3.0.0
colorama 0.4.6
comm 0.2.2
commonmark 0.9.1
conda 24.3.0
conda-libmamba-solver 24.1.0
conda-package-handling 2.2.0
conda_package_streaming 0.9.0
confection 0.1.4
contextily 1.6.0
contourpy 1.2.1
cryptography 42.0.5
cycler 0.12.1
cymem 2.0.8
cytoolz 0.12.3
dask 2024.4.1
dask-expr 1.0.10
debugpy 1.8.1
decorator 5.1.1
dill 0.3.8
distributed 2024.4.1
distro 1.9.0
docker 7.0.0
duckdb 0.10.1
en-core-web-sm 3.7.1
entrypoints 0.4
et-xmlfile 1.1.0
exceptiongroup 1.2.0
executing 2.0.1
fastjsonschema 2.19.1
fiona 1.9.6
flake8 7.0.0
Flask 3.0.2
folium 0.16.0
fontawesomefree 6.5.1
fonttools 4.51.0
frozenlist 1.4.1
fsspec 2023.12.2
GDAL 3.8.4
gensim 4.3.2
geographiclib 2.0
geopandas 0.12.2
geoplot 0.5.1
geopy 2.4.1
gitdb 4.0.11
GitPython 3.1.43
google-auth 2.29.0
graphene 3.3
graphql-core 3.2.3
graphql-relay 3.2.0
graphviz 0.20.3
great-tables 0.5.0
greenlet 3.0.3
gunicorn 21.2.0
htmltools 0.5.1
hvac 2.1.0
idna 3.6
imageio 2.34.1
importlib_metadata 7.1.0
importlib_resources 6.4.0
inflate64 1.0.0
ipykernel 6.29.3
ipython 8.22.2
ipywidgets 8.1.2
isort 5.13.2
itsdangerous 2.1.2
jedi 0.19.1
Jinja2 3.1.3
jmespath 1.0.1
joblib 1.3.2
jsonpatch 1.33
jsonpointer 2.4
jsonschema 4.21.1
jsonschema-specifications 2023.12.1
jupyter-cache 1.0.0
jupyter_client 8.6.1
jupyter_core 5.7.2
jupyterlab_widgets 3.0.10
kaleido 0.2.1
kiwisolver 1.4.5
kubernetes 29.0.0
langcodes 3.4.0
language_data 1.2.0
lazy_loader 0.4
libmambapy 1.5.7
llvmlite 0.42.0
locket 1.0.0
lxml 5.2.1
lz4 4.3.3
Mako 1.3.2
mamba 1.5.7
mapclassify 2.6.1
marisa-trie 1.1.0
Markdown 3.6
MarkupSafe 2.1.5
matplotlib 3.8.3
matplotlib-inline 0.1.6
mccabe 0.7.0
menuinst 2.0.2
mercantile 1.2.1
mizani 0.11.2
mlflow 2.11.3
mlflow-skinny 2.11.3
msgpack 1.0.7
multidict 6.0.5
multivolumefile 0.2.3
munkres 1.1.4
murmurhash 1.0.10
mypy 1.9.0
mypy-extensions 1.0.0
nbclient 0.10.0
nbformat 5.10.4
nest_asyncio 1.6.0
networkx 3.3
nltk 3.8.1
numba 0.59.1
numpy 1.26.4
oauthlib 3.2.2
opencv-python-headless 4.9.0.80
openpyxl 3.1.2
OWSLib 0.28.1
packaging 23.2
pandas 2.2.1
paramiko 3.4.0
parso 0.8.4
partd 1.4.1
pathspec 0.12.1
patsy 0.5.6
Pebble 5.0.7
pexpect 4.9.0
pickleshare 0.7.5
pillow 10.3.0
pip 24.0
pkgutil_resolve_name 1.3.10
platformdirs 4.2.0
plotly 5.19.0
plotnine 0.13.5
pluggy 1.4.0
polars 0.20.18
preshed 3.0.9
prometheus_client 0.20.0
prometheus-flask-exporter 0.23.0
prompt-toolkit 3.0.42
protobuf 4.25.3
psutil 5.9.8
ptyprocess 0.7.0
pure-eval 0.2.2
py7zr 0.20.8
pyarrow 15.0.0
pyarrow-hotfix 0.6
pyasn1 0.5.1
pyasn1-modules 0.3.0
pybcj 1.0.2
pycodestyle 2.11.1
pycosat 0.6.6
pycparser 2.21
pycryptodomex 3.20.0
pydantic 2.7.1
pydantic_core 2.18.2
pyflakes 3.2.0
Pygments 2.17.2
PyJWT 2.8.0
pylint 3.1.0
PyNaCl 1.5.0
pynsee 0.1.7
pyOpenSSL 24.0.0
pyparsing 3.1.2
pyppmd 1.1.0
pyproj 3.6.1
pyshp 2.3.1
PySocks 1.7.1
python-dateutil 2.9.0
python-dotenv 1.0.1
python-magic 0.4.27
pytz 2024.1
pyu2f 0.1.5
pywaffle 1.1.0
PyYAML 6.0.1
pyzmq 25.1.2
pyzstd 0.15.10
QtPy 2.4.1
querystring-parser 1.2.4
rasterio 1.3.10
referencing 0.34.0
regex 2023.12.25
requests 2.31.0
requests-cache 1.2.0
requests-oauthlib 2.0.0
rpds-py 0.18.0
rsa 4.9
Rtree 1.2.0
ruamel.yaml 0.18.6
ruamel.yaml.clib 0.2.8
s3fs 2023.12.2
s3transfer 0.10.1
scikit-image 0.23.2
scikit-learn 1.4.1.post1
scipy 1.13.0
seaborn 0.13.2
setuptools 69.2.0
shapely 2.0.3
six 1.16.0
smart-open 6.4.0
smmap 5.0.0
snuggs 1.4.7
sortedcontainers 2.4.0
soupsieve 2.5
spacy 3.7.4
spacy-legacy 3.0.12
spacy-loggers 1.0.5
SQLAlchemy 2.0.29
sqlparse 0.4.4
srsly 2.4.8
stack-data 0.6.2
statsmodels 0.14.1
tabulate 0.9.0
tblib 3.0.0
tenacity 8.2.3
texttable 1.7.0
thinc 8.2.3
threadpoolctl 3.4.0
tifffile 2024.4.24
tomli 2.0.1
tomlkit 0.12.4
toolz 0.12.1
topojson 1.8
tornado 6.4
tqdm 4.66.2
traitlets 5.14.2
truststore 0.8.0
typer 0.9.4
typing_extensions 4.11.0
tzdata 2024.1
Unidecode 1.3.8
url-normalize 1.4.3
urllib3 1.26.18
wasabi 1.1.2
wcwidth 0.2.13
weasel 0.3.4
webcolors 1.13
webdriver-manager 4.0.1
websocket-client 1.7.0
Werkzeug 3.0.2
wheel 0.43.0
widgetsnbextension 4.0.10
wordcloud 1.9.3
wrapt 1.16.0
xgboost 2.0.3
xlrd 2.0.1
xyzservices 2024.4.0
yarl 1.9.4
yellowbrick 1.5
zict 3.0.0
zipp 3.17.0
zstandard 0.22.0

View file history

SHA Date Author Description
06d003a 2024-04-23 10:09:22 Lino Galiana Continue la restructuration des sous-parties (#492)
8c316d0 2024-04-05 19:00:59 Lino Galiana Fix cartiflette deprecated snippets (#487)
ce33d5d 2024-01-16 15:47:22 Lino Galiana Adapte les exemples de code de cartiflette (#482)
005d89b 2023-12-20 17:23:04 Lino Galiana Finalise l’affichage des statistiques Git (#478)
3fba612 2023-12-17 18:16:42 Lino Galiana Remove some badges from python (#476)
1f23de2 2023-12-01 17:25:36 Lino Galiana Stockage des images sur S3 (#466)
69cf52b 2023-11-21 16:12:37 Antoine Palazzolo [On-going] Suggestions chapitres modélisation (#452)
09654c7 2023-11-14 15:16:44 Antoine Palazzolo Suggestions Git & Visualisation (#449)
8728352 2023-10-24 11:50:08 Lino Galiana Correction coquilles geopandas (#444)
102ce9f 2023-10-22 11:39:37 Thomas Faria Relecture Thomas, première partie (#438)
a771183 2023-10-09 11:27:45 Antoine Palazzolo Relecture TD2 par Antoine (#418)
f8831e7 2023-10-09 10:53:34 Lino Galiana Relecture antuki geopandas (#429)
154f09e 2023-09-26 14:59:11 Antoine Palazzolo Des typos corrigées par Antoine (#411)
9a4e226 2023-08-28 17:11:52 Lino Galiana Action to check URL still exist (#399)
a8f90c2 2023-08-28 09:26:12 Lino Galiana Update featured paths (#396)
a9198de 2023-08-25 18:33:00 Lino Galiana Geopandas tutorial
3bdf3b0 2023-08-25 11:23:02 Lino Galiana Simplification de la structure 🤓 (#393)
130ed71 2023-07-18 19:37:11 Lino Galiana Restructure les titres (#374)
f0c583c 2023-07-07 14:12:22 Lino Galiana Images viz (#371)
f21a24d 2023-07-02 10:58:15 Lino Galiana Pipeline Quarto & Pages 🚀 (#365)
8d81b5f 2023-02-18 18:21:59 Lino Galiana Change source get_vectorfile (#355)
d2eb6c2 2023-02-18 17:15:35 Lino Galiana Update index.qmd
3912a7e 2023-02-07 17:18:25 Lino Galiana Back to IGN provider (#350)
0312041 2022-12-11 13:43:49 Lino Galiana reprend box de geopandas (#332)
6662800 2022-10-28 11:14:27 Lino Galiana Change IGN dataset provider (#308)
f394b23 2022-10-13 14:32:05 Lino Galiana Dernieres modifs geopandas (#298)
8df6bbc 2022-10-12 11:50:57 Lino Galiana Corrige les tags du tuto geopandas (#295)
af763cc 2022-10-12 10:17:56 Lino Galiana Reprise exercice geopandas (#294)
1ef97df 2022-10-11 12:14:03 Lino Galiana Relecture chapitre geopandas (#289)
f10815b 2022-08-25 16:00:03 Lino Galiana Notebooks should now look more beautiful (#260)
494a85a 2022-08-05 14:49:56 Lino Galiana Images featured ✨ (#252)
d201e3c 2022-08-03 15:50:34 Lino Galiana Pimp la homepage ✨ (#249)
12965ba 2022-05-25 15:53:27 Lino Galiana :launch: Bascule vers quarto (#226)
9c71d6e 2022-03-08 10:34:26 Lino Galiana Plus d’éléments sur S3 (#218)
5cac236 2021-12-16 19:46:43 Lino Galiana un petit mot sur mercator (#201)
77120b8 2021-11-01 20:28:28 Lino Galiana Ajoute une section sur le geocodage (#173)
6777f03 2021-10-29 09:38:09 Lino Galiana Notebooks corrections (#171)
2a8809f 2021-10-27 12:05:34 Lino Galiana Simplification des hooks pour gagner en flexibilité et clarté (#166)
735e677 2021-10-19 09:46:12 Lino Galiana Règle problème des cartes qui s’affichent pas (#165)
5ad057f 2021-10-10 15:13:16 Lino Galiana Relectures pandas & geopandas (#159)
2e4d586 2021-09-02 12:03:39 Lino Galiana Simplify badges generation (#130)
4cdb759 2021-05-12 10:37:23 Lino Galiana :sparkles: :star2: Nouveau thème hugo :snake: :fire: (#105)
7f9f97b 2021-04-30 21:44:04 Lino Galiana 🐳 + 🐍 New workflow (docker 🐳) and new dataset for modelization (2020 🇺🇸 elections) (#99)
6d010fa 2020-09-29 18:45:34 Lino Galiana Simplifie l’arborescence du site, partie 1 (#57)
66f9f87 2020-09-24 19:23:04 Lino Galiana Introduction des figures générées par python dans le site (#52)
badc492 2020-09-22 18:36:33 Lino Galiana Finalize geopandas section (#48)
ffb05cf 2020-09-10 17:18:15 Lino Galiana Partie sur les données spatiales (#20) :warning: pas fini
Back to top

Footnotes

  1. Il est techniquement possible d’avoir un DataFrame comportant plusieurs géographies. Par exemple une géométrie polygone et une géométrie point (le centroid). C’est néanmoins souvent compliqué à gérer et donc peu recommandable.↩︎

Citation

BibTeX citation:
@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}
}
For attribution, please cite this work as:
Galiana, Lino. 2023. Python Pour La Data Science. https://doi.org/10.5281/zenodo.8229676.