03_geopandas_TP.qmd

---
title: "Pratique de geopandas avec les données vélib"
draft: false
weight: 50
code-annotations: true
slug: geopandasTP
tags:
  - geopandas
  - Velib
  - Exercice
  - Cartographie
  - Manipulation
categories:
  - Manipulation
  - Exercice
type: book
description: |
  Ce chapitre illustre les fonctionalités de `GeoPandas` à partir des
  décomptes de vélo fournis par la ville de Paris
  en [opendata](https://opendata.paris.fr/explore/dataset/comptage-velo-donnees-compteurs/map/?disjunctive.id_compteur&disjunctive.nom_compteur&disjunctive.id&disjunctive.name&basemap=jawg.dark&location=12,48.85855,2.33754).
  Il prolonge
  le chapitre précédent avec des données un petit peu plus complexes
  à manipuler.
echo: false
image: https://minio.lab.sspcloud.fr/lgaliana/generative-art/pythonds/velib.png
---

::: {.content-visible when-format="html"}
{{< include "../../build/_printBadges.qmd" >}}
:::

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[^noteQGIS].
Grâce à la librairie [`Geopandas`](https://geopandas.org/en/stable/), 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. 

::: {.cell .markdown}
```{=html}
<details>
<summary>

Illustration du principe des données spatiales (documentation de `sf`, l'équivalent de `Geopandas` en `R`)

</summary>

![](https://user-images.githubusercontent.com/520851/50280460-e35c1880-044c-11e9-9ed7-cc46754e49db.jpg){width="70%"}

</details>
```
:::
[^noteQGIS]: 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.

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`](https://rgeo.linogaliana.fr/exercises/geospatial-wrangling.html).

::: {.cell .markdown}
```{=html}
<div class="alert alert-info" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-comment"></i> Note</h3>
```

Le package `cartiflette` est expérimental 
et n'est disponible que sur
[`Github`](https://github.com/InseeFrLab/cartogether), 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):

```{python}
#| eval: false
!pip install py7zr geopandas openpyxl tqdm s3fs
!pip install PyYAML xlrd
!pip install git+https://github.com/inseefrlab/cartiflette
```

Ces commandes permettent de récupérer l'ensemble du code 
source depuis [`Github`](https://github.com/InseeFrLab/cartiflette)

```{=html}
</div>
```
:::

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

```{python}
#| echo: true
#| eval: false

!pip install fiona shapely
!pip install pyproj rtree
!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 :

```{python}
#| echo: true
#| output: false
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.

```{python}
#| echo: true
#| output: false
!pip install requests py7zr geopandas openpyxl tqdm s3fs PyYAML xlrd
!pip install git+https://github.com/inseefrlab/cartiflette
```

```{python}
#| echo: true
#| output: false
from cartiflette import carti_download
```




## 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`](https://github.com/InseeFrLab/cartiflette.git)[^cartiflette-python].

[^cartiflette-python]:
  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](https://linogaliana.github.io/cartiflette-website/index.html).

::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 1: découverte des objets géographiques</h3>
```

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

```{python}
#| echo: true
#| output: false
communes_borders = carti_download(
    crs = 4326,
    values = ["75", "92", "93", "94"],
    borders="COMMUNE",
    vectorfile_format="geojson",
    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. 

```{python}
# 2) Regarder les premières lignes
communes_borders.head()
```

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. 

```{=html}
</div>
```
:::


```{python}
#| output: false

# 1. Chargement des données de Cartiflette
communes_borders = carti_download(
    crs = 4326,
    values = ["75", "92", "93", "94"],
    borders="COMMUNE",
    vectorfile_format="geojson",
    filter_by="DEPARTEMENT",
    source="EXPRESS-COG-CARTO-TERRITOIRE",
    year=2022)
```


```{python}
#| output: false
# 2. Regarder les premières lignes
communes_borders.head()
# Il y a une colonne geometry qui contient les informations nécessaires pour connaître les contours communaux
```

```{python}
#| output: false

# 3) Afficher le crs
communes_borders.crs
# Les données sont en WGS84, on les reprojette en lambert 93
```

```{python}
#| echo: false

# 4) afficher les communes du département 92
ax = communes_borders[communes_borders['INSEE_DEP'] == "92"].boundary.plot()
ax.set_axis_off()
```

```{python}
#| echo: false

# 5) Représenter la carte de Paris. Quel est le problème ?
ax = communes_borders[communes_borders['INSEE_DEP'] == "75"].boundary.plot()
ax.set_axis_off()
#On remarque ainsi facilement le problème pour Paris
#intra-muros: il manque les limites des arrondissements. 
#Cela appauvrit grandement la carte de Paris. 
```


```{python}
#| output: false

# 6. Chargement des données de Cartiflette
petite_couronne = 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)
    
petite_couronne.crs
petite_couronne = petite_couronne.to_crs(2154)
petite_couronne.crs
```

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

::: {.content-visible when-format="html"}

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

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

{{< include "_play_with_projection.qmd" >}}

```{ojs}
width_projected_map = screen.width/2
```

:::

::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 2 : Les projections, représentations et approximations</h3>
```

Voici un code utilisant encore
`cartiflette` 
pour récupérer les frontières françaises (découpées par région):

```{python}
#| output: false
#| echo: true
france = carti_download(
      values = ["France"],
      crs = 4326,
      borders = "REGION",
      vectorfile_format="geojson",
      simplification=50,
      filter_by="FRANCE_ENTIERE",
      source="EXPRESS-COG-CARTO-TERRITOIRE",
      year=2022)
france = france.loc[france['INSEE_REG']>10]
```

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.


```{=html}
</div>
```
:::


```{python}
#| output: false
# Question 1 : Tester différentes projections
france_2154 = france.to_crs(2154)
france_healpix = france.to_crs("+proj=healpix +lon_0=0 +a=1")
france_5070 = france.to_crs(5070)
france_3304 = france.to_crs(3304)
```

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 :

```{python}
#| label: fig-effet-mercator
#| fig-cap: "Comparaison des projections"
#| fig-subcap: 
#|   - "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)"
#| layout-ncol: 2

ax1 = france_2154.boundary.plot(edgecolor = "k", linewidth=0.5)
ax2 = france_healpix.boundary.plot(edgecolor = "k", linewidth=0.5)
ax3 = france_5070.boundary.plot(edgecolor = "k", linewidth=0.5)
ax4 = france_3304.boundary.plot(edgecolor = "k", linewidth=0.5)

ax1.set_axis_off()
ax2.set_axis_off()
ax3.set_axis_off()
ax4.set_axis_off()
```

```{python}
# Question 2
france = france.to_crs(3395)
france["superficie_4326"] = france.area
france = france.to_crs(2154)
france["superficie_2154"] = france.area
france["mismatch"] = france['superficie_2154']-france['superficie_4326']
```


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[^notecarte] pour se faire
une idée des régions où l'erreur de mesure est la plus importante.

[^notecarte]: Cette carte n'est pas trop soignée, c'est normal nous verrons comment
faire de belles cartes ultérieurement.

```{python}
ax = france.plot(column = "mismatch")
ax.set_axis_off()
```

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. 

::: {#fig-mercator}

![Exemple de reprojection de pays depuis le site [thetruesize.com](https://www.thetruesize.com/)](https://pythonds.linogaliana.fr/content/manipulation/truesize.png){#fig-surus}

!["Don't trust the Mercator projection" sur `Reddit`](https://rgeo.linogaliana.fr/exercises/img/mercator.jpg){#fig-mercator-funny}

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



Pour aller plus loin, la carte interactive
suivante, construite par Nicolas Lambert, issue de
ce [_notebook_ `Observable`](https://observablehq.com/@neocartocnrs/impact-of-projections-on-areas), illustre l'effet
déformant de la projection Mercator, et de quelques-unes autres,
sur notre perception de la taille des pays.

::: {.content-visible when-format="html"}

<details>
<summary>
Voir la carte interactive
</summary>

```{ojs}
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>`
```

</details>

```{ojs}
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](https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/table/) et 
requêtables directement en utilisant un URL

```{python}
#| echo: true
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:

```{python}
#| echo: true
#| output: false
idf = carti_download(
      values = ["11"],
      crs = 4326,
      borders = "DEPARTEMENT",
      vectorfile_format="geojson",
      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)
)
```

::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 3: importer et explorer les données velib</h3>
```

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](https://geopandas.org/en/stable/docs/user_guide/mapping.html#maps-with-layers)

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:

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

```{=html}
</div>
```
:::


```{python}
#| output: false
# 1. Importer les données velib
stations = gpd.read_file(url)
```

```{python}
#| output: false
# 2. Reprojection
stations.crs
stations = stations.to_crs(petite_couronne.crs) 
```


```{python}
#| output: false
# 3. Principales stations
principales_stations = stations.sort_values("capacity", ascending = False).head(50)
```

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

```{python}
# 4. petite couronne
base = petite_couronne.boundary.plot(edgecolor = "black", linewidth = 0.5)
base
```


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

```{python}
# 5. Ajout de la couche des départements
base = petite_couronne.boundary.plot(edgecolor = "black", linewidth = 0.5)
petite_couronne_departements.boundary.plot(ax = base, edgecolor = "blue", linewidth = 0.7)
base
```

Puis les stations (question 6).

```{python}
# 6. Ajout des stations
base = petite_couronne.boundary.plot(edgecolor = "black", linewidth = 0.5)
petite_couronne_departements.boundary.plot(ax = base, edgecolor = "blue", linewidth = 0.7)
principales_stations.plot(ax= base, markersize = "capacity", color = "red", alpha = 0.4)
base
```

La carte finale, après mise en forme:


```{python}
#7. sans axe et avec titre
base = petite_couronne.boundary.plot(edgecolor = "black", linewidth = 0.5)
petite_couronne_departements.boundary.plot(ax = base, edgecolor = "blue", linewidth = 0.7)
principales_stations.plot(ax= base, markersize = "capacity", color = "red", alpha = 0.4)
base.set_axis_off()
base.set_title("Les 50 principales stations de Vélib")
base
```

## 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](https://geopandas.org/mergingdata.html#spatial-joins). 
Une version pédagogique pour `R` se trouve dans la documentation [`utilitR`](https://www.book.utilitr.org/03_fiches_thematiques/fiche_donnees_spatiales#joindre-des-donn%C3%A9es-g%C3%A9ographiques-et-attributaires). 


::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 4: Associer les stations aux communes et arrondissements auxquels elles appartiennent</h3>
```

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.
2. Représenter la carte des stations du 19e arrondissement avec le code suivant :

```python
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)


```{=html}
</div>
```
:::

```{python}
#1.  Jointure spatiale entre stations et data_paris
stations_info = gpd.sjoin(stations, petite_couronne, predicate = 'within')
```

```{python}
#2. 19e arrondissement
stations_19 = stations_info.loc[stations_info['NOM'].str.contains("19e")]
arrondissement_19e = petite_couronne.loc[petite_couronne['NOM'].str.contains("19e")]
```

```{python}
# 3. Carto du 19e 
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
```



Carte obtenue à la question 4 :

```{python}
#4. Calcul et carte des capacity
stations_agg = (
  stations_info
  .groupby("INSEE_COG")
  .agg({"stationcode": "nunique", "capacity": "sum"})
  .reset_index()
)

petite_couronne_count = petite_couronne.merge(
  stations_agg
).to_crs(2154)

petite_couronne_count = petite_couronne_count.loc[petite_couronne_count["INSEE_DEP"]== "75"]
aplat = petite_couronne_count.plot(
    column = "capacity", cmap="coolwarm", legend=True)
aplat.set_axis_off()
aplat
```

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.

```{python}
#5. Calcul et carte des area et density
petite_couronne_count['area'] = petite_couronne_count.area
petite_couronne_count['area'] = petite_couronne_count['area'].div(1e6)
petite_couronne_count['density'] = petite_couronne_count['capacity']/petite_couronne_count['area']

aplat = petite_couronne_count.plot(
    column = "density", cmap="coolwarm", legend=True)
aplat.set_axis_off()
aplat
```


## 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ï](https://fr.wikipedia.org/wiki/Diagramme_de_Vorono%C3%AF),
une opération classique pour transformer des points en polygones. L'exercice suivant
permet de se familiariser avec cette approche[^notevoronoi].


[^notevoronoi]: Dans [ce document de travail](https://www.insee.fr/en/statistiques/4925202) 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.


__Exercice à venir__