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

• Localisations des stations velib ;
• fonds de carte AdminExpress à travers un package Python nommé cartiflette facilitant la récupération de cette source.

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

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

/miniconda/envs/python-ENSAE/lib/python3.9/site-packages/geopandas/_compat.py:123: UserWarning:

The Shapely GEOS version (3.11.0-CAPI-1.17.0) is incompatible with the GEOS version PyGEOS was compiled with (3.10.3-CAPI-1.16.1). Conversions between both will be slow.

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 ;
• 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.

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 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¹, il est possible de se rendre sur le site de l’IGN et de le télécharger depuis le serveur FTP. 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 pynsee propose notamment un module dédié à la récupération de fonds de carte officiels pour valoriser des données d’open data. L’API sur laquelle il repose étant parfois lente, nous présentons le code dédié uniquement en annexe.

Nous proposons ici une méthode nouvelle de récupération de ces données qui s’appuie sur le projet interministériel cartiflette.
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.

!pip install --upgrade botocore==1.23.26  #Sur colab, sinon bug
!pip install --upgrade urllib3==1.22.0 #Sur colab, sinon bug
!pip install py7zr #Sur colab, sinon bug
!pip install s3fs #Sur colab, sinon bug
!git clone https://github.com/InseeFrLab/cartogether.git
%cd ./cartogether
!pip install -r requirements.txt
!pip install .

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:

import cartiflette.s3

shp_communes = cartiflette.s3.download_vectorfile_url_all(
    values = ["75", "92", "93", "94"],
    level="COMMUNE_ARRONDISSEMENT",
    vectorfile_format="geojson",
    decoupage="departement",
    year=2022)

shp_communes = shp_communes
shp_communes.head()

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

Note

Si on ne désire pas utiliser le niveau COMMUNE_ARRONDISSEMENT, il est nécessaire de mettre en oeuvre une construction du fond de carte en plusieurs phases.

En premier lieu, il est nécessaire de récupérer le niveau des communes.

import cartiflette.s3

shp_communes = cartiflette.s3.download_vectorfile_url_all(
    values = ["75", "92", "93", "94"],
    level="COMMUNE",
    vectorfile_format="GPKG",
    decoupage="departement",
    year=2022)

shp_communes.head()

	ID	NOM	NOM_M	INSEE_COM	STATUT	POPULATION	INSEE_CAN	INSEE_ARR	INSEE_DEP	INSEE_REG	SIREN_EPCI	geometry
0	COMMUNE_0000000009736048	Paris	PARIS	75056	Capitale d'état	2165423	NR	1	75	11	200054781	POLYGON ((653313.800 6857558.100, 653249.300 6...
0	COMMUNE_0000000009736037	Levallois-Perret	LEVALLOIS-PERRET	92044	Commune simple	66082	16	2	92	11	200054781/200057982	POLYGON ((647761.400 6867306.900, 647839.200 6...
1	COMMUNE_0000000009736055	Bois-Colombes	BOIS-COLOMBES	92009	Commune simple	28841	11	2	92	11	200054781/200057990	POLYGON ((646224.700 6867615.800, 646229.000 6...
2	COMMUNE_0000000009736538	Malakoff	MALAKOFF	92046	Commune simple	30950	18	1	92	11	200054781/200057966	POLYGON ((646995.300 6857373.400, 647177.500 6...
3	COMMUNE_0000000009736038	Clichy	CLICHY	92024	Commune simple	63089	09	2	92	11	200054781/200057990	POLYGON ((648952.800 6866241.200, 648865.900 6...

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

On peut remarquer que la ville de Paris ne comporte pas d’arrondissements sur cette carte. Pour vous en convaincre, vous pouvez exécuter la commande :

ax = shp_communes.loc[shp_communes['INSEE_DEP']=="75"].boundary.plot()
ax.set_axis_off()

Il faut donc utiliser une source complémentaire. Le contour officiel des arrondissements est produit par l’IGN séparemment des contours de communes. Les contours d’arrondissements sont également disponibles grâce à cartiflette:

arrondissements = cartiflette.s3.download_vectorfile_url_all(
    values = "75",
    level="ARRONDISSEMENT_MUNICIPAL",
    vectorfile_format="geojson",
    decoupage="departement",
    year=2022)

ax = arrondissements.plot(alpha = 0.8, edgecolor = "k")
ax.set_axis_off()

Il ne reste plus qu’à remplacer Paris par ses arrondissements dans shp_communes. Pour cela, on peut utiliser les méthodes vues dans le chapitre Pandas relatives aux filtres et à la concaténation de plusieurs DataFrames:

import pandas as pd

shp_communes = pd.concat(
  [
    shp_communes[shp_communes['INSEE_DEP'] != "75"].to_crs(2154),
    arrondissements.to_crs(2154)
  ])

Cette approche fonctionne mais elle nécessite un certain nombre de gestes, qui sont autant de risques d’erreurs. Il est donc recommandé de privilégier le niveau COMMUNE_ARRONDISSEMENT qui fait exactement ceci mais de manière fiable.

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

Import des données velib

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 Stamen, 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.Stamen.Watercolor)
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)

<AxesSubplot: >

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

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.Stamen.Watercolor)
ax.set_axis_off()
ax

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

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:

stations.describe()

Pour connaître les plus grands départements de France métropolitaine, procédons en deux étapes:

1. Récupérons le contour des communes de France métropolitaine dans son ensemble 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.
2. Calculons la surface (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)

from cartiflette.download import get_vectorfile_ign
france = get_vectorfile_ign(
  level = "COMMUNE",
  field = "metropole",
  source = "COG",
  provider="opendatarchives"
  )

france['surface'] = france.area.div(10**6)

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

france.groupby('INSEE_DEP').sum(numeric_only = True).sort_values('surface', ascending = False)

Si on veut directement les plus grandes communes de France métropolitaine :

france.sort_values('surface', ascending = False).head(10)

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))
france.dissolve(by='INSEE_DEP', aggfunc='sum').plot(ax = ax, column = "surface")
ax.set_axis_off()
ax

/miniconda/envs/python-ENSAE/lib/python3.9/site-packages/geopandas/geodataframe.py:1676: FutureWarning:

The default value of numeric_only in DataFrameGroupBy.sum is deprecated. In a future version, numeric_only will default to False. Either specify numeric_only or select only columns which should be valid for the function.

<AxesSubplot: >

Pour produire cette carte, il serait 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 = get_vectorfile_ign(
  level = "DEPARTEMENT", field = "metropole",
  source = "COG", provider="opendatarchives")
dep["area"] = dep.area
ax = dep.plot(column = "area")
ax.set_axis_off()

opendatarchives
COG
Data have been previously downloaded and are still available in /tmp/COG-2022

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)

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), on fera

communes_31 = cartiflette.s3.download_vectorfile_url_all(
    values = ["31"],
    level="COMMUNE",
    vectorfile_format="geojson",
    decoupage="departement",
    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

<AxesSubplot: >

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)

Le système de projection est fondamental pour que la dimension spatiale soit bien interprétée par Python. 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. Ce post propose de riches éléments sur le sujet, notamment l’image suivante qui montre bien le principe d’une projection :

La Terre peut ainsi être représentée de multiples manières, ce qui n’est pas neutre sur la manière de se représenter certains continents. L’Afrique apparaît beaucoup moins vaste qu’elle ne l’est en réalité sur les cartes utilisant cette projection. L’une des déformations les mieux connue est celle provoquée par la projection Mercator. Le Groënland paraît avoir la même surface que l’Amérique du Sud. Pourtant, cette dernière est 8 fois plus grande.

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.

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.

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

communes.crs

<Derived 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 (⚠️ 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 = get_vectorfile_ign(
  level = "REGION", field = "metropole",
  source = "COG", provider="opendatarchives")

fig,ax = plt.subplots(figsize=(10, 10))
shp_region.to_crs(5070).plot(ax = ax)
ax

opendatarchives
COG
Data have been previously downloaded and are still available in /tmp/COG-2022

<AxesSubplot: >

<Figure size 768x480 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

<AxesSubplot: >

<Figure size 768x480 with 0 Axes>

Joindre des données

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éographique², 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.

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

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

Annexe

Récupération des données depuis data.gouv

Avec cette méthode, les données des limites administratives demandent donc un peu de travail pour être importées car elles sont zippées (mais c’est un bon exercice !).

Le code suivant, dont les détails apparaîtront plus clairs après la lecture de la partie webscraping permet de :

1. Télécharger les données avec requests dans un dossier temporaire
2. Les dézipper avec le module zipfile

La fonction suivante automatise un peu le processus :

import requests
import tempfile
import zipfile

url = 'https://www.data.gouv.fr/fr/datasets/r/0e117c06-248f-45e5-8945-0e79d9136165'
temporary_location = tempfile.gettempdir()

def download_unzip(url, dirname = tempfile.gettempdir(), destname = "borders"):
  myfile = requests.get(url)
  open("{}/{}.zip".format(dirname, destname), 'wb').write(myfile.content)
  with zipfile.ZipFile("{}/{}.zip".format(dirname, destname), 'r') as zip_ref:
      zip_ref.extractall(dirname + '/' + destname)

download_unzip(url)

shp_communes = gpd.read_file(temporary_location + "/borders/communes-20220101.shp")

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é :

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

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 :

arrondissements = gpd.read_file("https://opendata.paris.fr/explore/dataset/arrondissements/download/?format=geojson&timezone=Europe/Berlin&lang=fr")
arrondissements = arrondissements.rename(columns = {"c_arinsee": "insee"})
arrondissements['insee'] = arrondissements['insee'].astype(str)
shp_communes = shp_communes[~shp_communes.insee.str.startswith("75")].append(arrondissements)

Pour produire la carte, il faudrait faire:

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='k')
ctx.add_basemap(ax, crs = paris.crs.to_string())
ax.set_axis_off()
ax

Récupération des données depuis le package pynsee

Pour connaître les contraintes d’installation du package pynsee, se référer à la partie de cours dédiée à Pandas.

#le téléchargement des données prend plusieurs minutes
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)