!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 pygeos
!pip install topojson
Dans ce TP,
nous allons apprendre à importer et
manipuler des données spatiales avec
Python
.
Ce langage propose
des fonctionnalités très intéressantes pour ce type de
données complexes qui le rendent capable de se comporter
comme un logiciel de SIG1.
Grâce à la librairie Geopandas
, une extension
de Pandas
aux données spatiales, les
données géographiques pourront être manipulées
comme n’importe quel type de données avec Python
.
La complexité induite par la dimension spatiale ne sera pas ressentie.
Illustration du principe des données spatiales (documentation de `sf`, l'équivalent de `Geopandas` en `R`)
{width="70%"}Ce chapitre illustre à partir d’exemples pratiques certains principes centraux de l’analyse de données :
- Manipulations sur les attributs des jeux de données ;
- Manipulations géométriques ;
- Gestion des projections cartographiques ;
- Création rapide de cartes (ce sera approfondi dans un prochain chapitre).
Si vous êtes intéressés par R
,
une version très proche de ce TP est disponible dans ce cours de R
.
Note
Le package cartiflette
est expérimental
et n’est disponible que sur
Github
, pas sur PyPi
.
Il est amené à évoluer rapidement et cette page sera mise à jour
quand de nouvelles fonctionalités (notamment l’utilisation d’API
)
seront disponibles pour encore simplifier la récupération de
contours géographiques.
Pour installer cartiflette
, il est nécessaire d’utiliser les commandes suivantes
depuis un Jupyter Notebook
(si vous utilisez la ligne de commande directement,
vous pouvez retirer les !
en début de ligne):
Ces commandes permettent de récupérer l’ensemble du code
source depuis Github
Préliminaires
Avant de se lancer dans le TD, il est nécessaire d’installer quelques
librairies qui ne sont pas disponibles par défaut, dans l’environnement Python
de base de la data science. Pour installer celles-ci depuis une
cellule de notebook Jupyter
, le code suivant est à exécuter :
Après installations, les packages à importer pour progresser dans ce chapitre sont les suivants :
import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt
Les instructions d’installation du package cartiflette
sont quant à elles détaillées dans le chapitre
précédent.
!pip install requests py7zr geopandas openpyxl tqdm s3fs PyYAML xlrd
!pip install git+https://github.com/inseefrlab/cartiflette@80b8a5a28371feb6df31d55bcc2617948a5f9b1a
from cartiflette.s3 import download_vectorfile_url_all
Lire et enrichir des données spatiales
Dans cette partie,
nous utiliserons
les fonds de carte de l’IGN dont
la mise à disposition est facilitée
par le projet cartiflette
2.
Exercice 1: découverte des objets géographiques
En premier lieu, on récupère des données géographiques grâce
au package cartiflette
.
- Utiliser
le code ci-dessous pour
télécharger les données communales (produit
Admin Express
de l’IGN) des départements de la petite couronne (75, 92, 93 et 94) de manière simplifiée grâce au packagecartiflette
:
= download_vectorfile_url_all(
communes_borders = 4326,
crs = ["75", "92", "93", "94"],
values ="COMMUNE",
borders="topojson",
vectorfile_format="DEPARTEMENT",
filter_by="EXPRESS-COG-CARTO-TERRITOIRE",
source=2022) year
- Regarder les premières lignes des données. Identifier la différence avec un dataframe standard.
id | ID | NOM | NOM_M | INSEE_COM | STATUT | POPULATION | INSEE_CAN | INSEE_ARR | INSEE_DEP | INSEE_REG | SIREN_EPCI | source | geometry | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | COMMUNE_0000000009736048 | NaN | Paris | PARIS | 75056 | Capitale d'état | 2165423 | NR | 1 | 75 | 11 | 200054781 | IGN:EXPRESS-COG-CARTO-TERRITOIRE | POLYGON ((2.36421 48.81640, 2.36333 48.81615, ... |
0 | COMMUNE_0000000009736037 | NaN | Levallois-Perret | LEVALLOIS-PERRET | 92044 | Commune simple | 66082 | 16 | 2 | 92 | 11 | 200054781/200057982 | IGN:EXPRESS-COG-CARTO-TERRITOIRE | POLYGON ((2.28739 48.90364, 2.28846 48.90302, ... |
1 | COMMUNE_0000000009736055 | NaN | Bois-Colombes | BOIS-COLOMBES | 92009 | Commune simple | 28841 | 11 | 2 | 92 | 11 | 200054781/200057990 | IGN:EXPRESS-COG-CARTO-TERRITOIRE | POLYGON ((2.26639 48.90629, 2.26645 48.90615, ... |
2 | COMMUNE_0000000009736538 | NaN | Malakoff | MALAKOFF | 92046 | Commune simple | 30950 | 18 | 1 | 92 | 11 | 200054781/200057966 | IGN:EXPRESS-COG-CARTO-TERRITOIRE | POLYGON ((2.27818 48.81425, 2.28066 48.81469, ... |
3 | COMMUNE_0000000009736038 | NaN | Clichy | CLICHY | 92024 | Commune simple | 63089 | 09 | 2 | 92 | 11 | 200054781/200057990 | IGN:EXPRESS-COG-CARTO-TERRITOIRE | POLYGON ((2.30377 48.89415, 2.30258 48.89487, ... |
Afficher le
crs
decommunes_borders
. Ce dernier contrôle la transformation de l’espace tridimensionnel terrestre en une surface plane. Utiliserto_crs
pour transformer les données en Lambert 93, le système officiel (code EPSG 2154).Afficher les communes des Hauts de Seine (département 92) et utiliser la méthode
plot
Ne conserver que Paris et réprésenter les frontières sur une carte : quel est le problème pour une analyse de Paris intramuros?
On remarque rapidement le problème. On ne dispose ainsi pas des limites des arrondissements parisiens, ce qui appauvrit grandement la carte de Paris.
- Cette fois, utiliser l’argument
borders="COMMUNE_ARRONDISSEMENT"
pour obtenir un fonds de carte consolidé des communes avec les arrondissements dans les grandes villes. Convertir en Lambert 93.
Le système de projection
Un concept central dans les logiciels de SIG est la notion de
projection. L’exercice précédent imposait parfois certaines projections
sans expliquer l’importance de ces choix. Python
, comme
tout SIG, permet une gestion cohérente des projections.
Observez les variations significatives de proportions pour certains pays selon les projections choisies:
Exercice 2 : Les projections, représentations et approximations
Voici un code utilisant encore
cartiflette
pour récupérer les frontières françaises (découpées par région):
= download_vectorfile_url_all(
france = "metropole",
values = 4326,
crs = "REGION",
borders ="topojson",
vectorfile_format="FRANCE_ENTIERE",
filter_by="EXPRESS-COG-CARTO-TERRITOIRE",
source=2022) year
- S’amuser à représenter les limites de la France avec plusieurs projections:
- Mercator WGS84 (EPSG: 4326)
- Projection healpix (
+proj=healpix +lon_0=0 +a=1
) - Projection prévue pour Tahiti (EPSG: 3304)
- Projection Albers prévue pour Etats-Unis (EPSG: 5070)
- Calculer la superficie en \(km^2\) des régions françaises dans les deux systèmes de projection suivants : World Mercator WGS84 (EPSG: 3395) et Lambert 93 (EPSG: 2154). Calculer la différence en \(km^2\) pour chaque région.
Avec la question 1 illustrant quelques cas pathologiques, on comprend que les projections ont un effet déformant qui se voit bien lorsqu’on les représente côte à côte sous forme de cartes :
Cependant le problème n’est pas que visuel, il est également numérique. Les calculs géométriques amènent à des différences assez notables selon le système de référence utilisé.
On peut représenter ces approximations sur une carte3 pour se faire une idée des régions où l’erreur de mesure est la plus importante.
/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
Ce type d’erreur de mesure est normal à l’échelle du territoire français. Les projections héritères du Mercator déforment les distances, surtout lorqu’on se rapproche de l’équateur ou des pôles.
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 n’est donc pas suprenant que nos déformations soient exacerbées aux extrèmes du territoire métropolitain. Si les approximations sont légères sur de petits territoires, les erreurs peuvent être non négligeables à l’échelle de la France.
Il faut donc systématiquement repasser les données dans le système de projection Lambert 93 (le système officiel pour la métropole) avant d’effectuer des calculs géométriques.
Utiliser des données géographiques comme des couches graphiques
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 en utilisant un URL
= "https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr" url
Dans le prochain exercice, nous proposons de créer rapidement une carte comprenant trois couches :
- Les localisations de stations sous forme de points ;
- Les bordures des communes et arrondissements pour contextualiser ;
- Les bordures des départements en traits plus larges pour contextualiser également.
Nous irons plus loin dans le travail cartographique dans le prochain chapitre. Mais être en mesure de positionner rapidement ses données sur une carte est toujours utile dans un travail exploratoire.
En amont de l’exercice,
utiliser la fonction suivante du package cartiflette
pour récupérer
le fonds de carte des départements de la petite couronne:
= download_vectorfile_url_all(
idf = "11",
values = 4326,
crs = "DEPARTEMENT",
borders ="topojson",
vectorfile_format="REGION",
filter_by="EXPRESS-COG-CARTO-TERRITOIRE",
source=2022)
year
= idf.loc[idf['INSEE_DEP'].isin(["75","92","93","94"])].to_crs(2154) petite_couronne_departements
Exercice 3: importer et explorer les données velib
On commence par récupérer les données nécessaires à la production de cette carte.
- En utilisant l’URL précédent, importer les données velib sous le nom
station
- Vérifier la projection géographique de
station
(attributcrs
). Si celle-ci est différente des données communales, reprojeter ces dernières dans le même système de projection que les stations de vélib - Ne conserver que les 50 principales stations (variable
capacity
)
On peut maintenant construire la carte de manière séquentielle avec la méthode plot
en s’aidant de cette documentation
En premier lieu, grâce à
boundary.plot
, représenter la couche de base des limites des communes et arrondissements:- Utiliser les options
edgecolor = "black"
etlinewidth = 0.5
- Nommer cet objet
base
- Utiliser les options
Ajouter la couche des départements avec les options
edgecolor = "blue"
etlinewidth = 0.7
Ajouter les positions des stations et ajuster la taille en fonction de la variable
capacity
. L’esthétique des points obtenus peut être contrôlé grâce aux optionscolor = "red"
etalpha = 0.4
.Retirer les axes et ajouter un titre avec les options ci-dessous:
base.set_axis_off()"Les 50 principales stations de Vélib") base.set_title(
La couche de base obtenue à l’issue de la question 4.
Puis en y ajoutant les limites départementales (question 5).
Puis les stations (question 6).
La carte finale, après mise en forme:
Jointures spatiales
Les jointures attributaires fonctionnent comme avec un Pandas
classique.
Pour conserver un objet spatial in fine, il faut faire attention à utiliser en premier (base de gauche) l’objet Geopandas
.
En revanche, l’un des intérêts des objets Geopandas
est qu’on peut également faire une jointure sur la dimension spatiale grâce à plusieurs fonctions.
La documentation à laquelle se référer est ici.
Une version pédagogique pour R
se trouve dans la documentation utilitR
.
Exercice 4: Associer les stations aux communes et arrondissements auxquels elles appartiennent
Dans cet exercice, on va supposer que :
- les localisations des stations
velib
sont stockées dans un dataframe nomméstations
- les données administratives
sont dans un dataframe nommé
petite_couronne
.
- Faire une jointure spatiale pour enrichir les données de stations en y ajoutant des informations de
petite_couronne
. Appeler cet objetstations_info
. - Créer les objets
stations_19e
etarrondissement_19e
pour stocker, respectivement, les stations appartenant au 19e et les limites de l’arrondissement. - Représenter la carte des stations du 19e arrondissement avec le code suivant :
= petite_couronne.loc[petite_couronne['INSEE_DEP']=="75"].boundary.plot(edgecolor = "k", linewidth=0.5)
base = base, edgecolor = "red", linewidth=0.9)
arrondissement_19e.boundary.plot(ax = base, color = "red", alpha = 0.4)
stations_19.plot(ax
base.set_axis_off()"Les stations Vélib du 19e arrondissement")
base.set_title( base
- Compter le nombre de stations velib et le nombre de places velib par arrondissement ou commune. Représenter sur une carte chacune des informations
- Représenter les mêmes informations mais en densité (diviser par la surface de l’arrondissement ou commune en km2)
Carte obtenue à la question 4 :
/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
<Axes: >
Avec la carte de la question 4, basée sur des aplats de couleurs (choropleth map), le lecteur est victime d’une illusion classique. Les arrondissements les plus visibles sur la carte sont les plus grands. D’ailleurs c’est assez logique qu’ils soient également mieux pourvus en velib. Même si l’offre de velib est probablement plus reliée à la densité de population et d’équipements, on peut penser que l’effet taille joue et qu’ainsi on est victime d’une illusion avec la carte précédente.
Si on représente plutôt la capacité sous forme de densité, pour tenir compte de la taille différente des arrondissements, les conclusions sont inversées et correspondent mieux aux attentes d’un modèle centre-périphérie. Les arrondissements centraux sont mieux pourvus, cela se voit encore mieux avec des ronds proportionnels plutôt qu’une carte chorolèpthe.
Exercice supplémentaire
Les exercices précédents ont permis de se familiariser au traitement de données spatiales. Néanmoins il arrive de devoir jongler plus avec la dimension géométrique par exemple pour changer d’échelle ou introduire des fusions/dissolutions de géométries.
Imaginons que chaque utilisateur de velib se déplace exclusivement vers la station la plus proche (à supposer qu’il n’y a jamais pénurie ou surcapacité). Quelle est la carte de la couverture des vélibs ? Pour répondre à ce type de question, on utilise fréquemment la la tesselation de Voronoï, une opération classique pour transformer des points en polygones. L’exercice suivant permet de se familiariser avec cette approche4.
Exercice à venir
Footnotes
D’ailleurs, le logiciel de cartographie spécialisé QGIS, s’appuie sur
Python
pour les manipulations de données nécessaires avant de réaliser une carte.↩︎La librairie
Python
est encore expérimentale mais les prochaines semaines devraient permettre de combler ce manque. Une documentation interactive illustrant le code nécessaire pour reproduire telle ou telle carte est disponible sur linogaliana.github.io/cartiflette-website.↩︎Cette carte n’est pas trop soignée, c’est normal nous verrons comment faire de belles cartes ultérieurement.↩︎
Dans ce document de travail sur données de téléphonie mobile, on montre néanmoins que cette approche n’est pas sans biais sur des phénomènes où l’hypothèse de proximité spatiale est trop simplificatrice.↩︎
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}
}