Pratique de geopandas avec les données vélib

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 SIG¹. 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 :

!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

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

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.

1. 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 package cartiflette:

communes_borders = download_vectorfile_url_all(
    crs = 4326,
    values = ["75", "92", "93", "94"],
    borders="COMMUNE",
    vectorfile_format="topojson",
    filter_by="DEPARTEMENT",
    source="EXPRESS-COG-CARTO-TERRITOIRE",
    year=2022)

2. 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, ...

3. Afficher le crs de communes_borders. Ce dernier contrôle la transformation de l’espace tridimensionnel terrestre en une surface plane. Utiliser to_crs pour transformer les données en Lambert 93, le système officiel (code EPSG 2154).

4. Afficher les communes des Hauts de Seine (département 92) et utiliser la méthode plot

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

6. 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:

html`<div>${container_projection}</div>`

container_projection = html`<div class="container">
  <div class="row">
    <div class="projection">
      <div class="projection-label">Choisir une projection</div>
      ${viewof projection}
    </div>
  </div>
  <div class="row">
    <div class="projectedMap">
      ${projectedMap}
    </div>
  </div>
</div>`

viewof projection = projectionInput({
  name: "",
  value: "Mercator"
})

import {projectionInput} from "@fil/d3-projections"
import {map} from "@linogaliana/base-map"

projectedMap = map(projection,
                   {
                     //svg: true,
                     value: projection.options,
                     width: width_projected_map,
                     //height: 300,
                     //rotate: [0, -90],
                     //inertia: true,
                     show_equator: true,
                     background: "#f1f0eb"
                     
                     //show_structure: true
                   })

width_projected_map = screen.width/2

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

france = download_vectorfile_url_all(
      values = "metropole",
      crs = 4326,
      borders = "REGION",
      vectorfile_format="topojson",
      filter_by="FRANCE_ENTIERE",
      source="EXPRESS-COG-CARTO-TERRITOIRE",
      year=2022)

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

2. 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 :

(a) Mercator WGS84 (EPSG: 4326)

(b) Projection healpix (+proj=healpix +lon_0=0 +a=1)

(c) Projection prévue pour Tahiti (EPSG: 3304)

(d) Projection Albers prévue pour Etats-Unis (EPSG: 5070)

Figure 1: Comparaison des projections

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 carte³ 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.

(a) Exemple de reprojection de pays depuis le site thetruesize.com

(b) “Don’t trust the Mercator projection” sur Reddit

Figure 2: La projection Mercator, une vision déformante

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

html`<div class="grid-container">
  <div class="viewof-projection">${viewof projectionBertin}</div>
  <div class="viewof-mycountry">${viewof mycountry}</div>
  <div class="map-bertin">${mapBertin}</div>
</div>`

import {map as mapBertin, viewof projection as projectionBertin, viewof mycountry} from "@neocartocnrs/impact-of-projections-on-areas"

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

url = "https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr"

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:

idf = download_vectorfile_url_all(
      values = "11",
      crs = 4326,
      borders = "DEPARTEMENT",
      vectorfile_format="topojson",
      filter_by="REGION",
      source="EXPRESS-COG-CARTO-TERRITOIRE",
      year=2022)

petite_couronne_departements = idf.loc[idf['INSEE_DEP'].isin(["75","92","93","94"])].to_crs(2154)

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.

1. En utilisant l’URL précédent, importer les données velib sous le nom station
2. Vérifier la projection géographique de station (attribut crs). 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
3. 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

4. En premier lieu, grâce à boundary.plot, représenter la couche de base des limites des communes et arrondissements:

   • Utiliser les options edgecolor = "black" et linewidth = 0.5
   • Nommer cet objet base

5. Ajouter la couche des départements avec les options edgecolor = "blue" et linewidth = 0.7

6. 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 options color = "red" et alpha = 0.4.

7. Retirer les axes et ajouter un titre avec les options ci-dessous:

base.set_axis_off()
base.set_title("Les 50 principales stations de Vélib")

La couche de base obtenue à l’issue de la question 4.

<Axes: >

Puis en y ajoutant les limites départementales (question 5).

<Axes: >

Puis les stations (question 6).

<Axes: >

La carte finale, après mise en forme:

<Axes: title={'center': 'Les 50 principales stations de Vélib'}>

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.

1. Faire une jointure spatiale pour enrichir les données de stations en y ajoutant des informations de petite_couronne. Appeler cet objet stations_info.
2. Créer les objets stations_19e et arrondissement_19e pour stocker, respectivement, les stations appartenant au 19e et les limites de l’arrondissement.
3. Représenter la carte des stations du 19e arrondissement avec le code suivant :

base = petite_couronne.loc[petite_couronne['INSEE_DEP']=="75"].boundary.plot(edgecolor = "k", linewidth=0.5)
arrondissement_19e.boundary.plot(ax = base, edgecolor = "red", linewidth=0.9)
stations_19.plot(ax = base, color = "red", alpha = 0.4)
base.set_axis_off()
base.set_title("Les stations Vélib du 19e arrondissement")
base

3. Compter le nombre de stations velib et le nombre de places velib par arrondissement ou commune. Représenter sur une carte chacune des informations
4. Représenter les mêmes informations mais en densité (diviser par la surface de l’arrondissement ou commune en km2)

<Axes: title={'center': 'Les stations Vélib du 19e arrondissement'}>

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.

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

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 approche⁴.

Exercice à venir

Footnotes

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

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

3. Cette carte n’est pas trop soignée, c’est normal nous verrons comment faire de belles cartes ultérieurement.

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

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.