-
Notifications
You must be signed in to change notification settings - Fork 45
/
03_geopandas_tutorial.Rmd
696 lines (519 loc) · 29.1 KB
/
03_geopandas_tutorial.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
---
jupyter:
jupytext:
text_representation:
extension: .Rmd
format_name: rmarkdown
format_version: '1.2'
jupytext_version: 1.6.0
kernelspec:
display_name: Python 3
language: python
name: python3
title: "Données spatiales: découverte de geopandas"
date: 2020-07-09T13:00:00Z
draft: false
weight: 40
slug: geopandas
tags:
- geopandas
- cartes
- velib
categories:
- Tutoriel
type: book
summary: |
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 complexe de manière simple.
---
```{r setup, include=FALSE}
dir_path <- gsub(here::here(), "..", here::here("course","manipulation"))
library(knitr)
library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
knitr::opts_chunk$set(fig.path = "")
knitr::opts_chunk$set(eval = TRUE, warning = FALSE, message = FALSE)
# Hook from Maelle Salmon: https://ropensci.org/technotes/2020/04/23/rmd-learnings/
knitr::knit_hooks$set(
plot = function(x, options) {
hugoopts <- options$hugoopts
paste0(
"{", "{<figure src=", # the original code is simpler
# but here I need to escape the shortcode!
'"', paste0(dir_path,"/",x), '" ',
if (!is.null(hugoopts)) {
glue::glue_collapse(
glue::glue('{names(hugoopts)}="{hugoopts}"'),
sep = " "
)
},
">}}\n"
)
}
)
```
```{r, echo = FALSE, results = 'asis', include = TRUE, eval = TRUE}
print_badges()
```
Dans ce tutoriel, nous allons utiliser:
<!----
* [Réseau de pistes cyclables de la ville de Paris](https://opendata.paris.fr/explore/dataset/reseau-cyclable/table/?disjunctive.typologie_simple&disjunctive.bidirectionnel&disjunctive.statut&disjunctive.sens_velo&disjunctive.voie&disjunctive.arrdt&disjunctive.bois&disjunctive.position&disjunctive.circulation&disjunctive.piste&disjunctive.couloir_bus&disjunctive.type_continuite&disjunctive.reseau&basemap=jawg.streets&location=12,48.85943,2.3493)
------>
* [Localisations des stations velib](https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
* [Carte des limites administratives françaises](https://www.data.gouv.fr/fr/datasets/r/07b7c9a2-d1e2-4da6-9f20-01a7b72d4b12)
* [Arrondissements parisiens](https://opendata.paris.fr/explore/dataset/arrondissements/download/?format=geojson&timezone=Europe/Berlin&lang=fr)
La représentation des données, notamment la cartographie, est présentée plus
amplement dans la partie [visualiser](#visualiser). Quelques méthodes
pour faire rapidement des cartes seront présentées ici mais
ce n'est pas l'objet de celle-ci.
Ce tutoriel s'inspire beaucoup d'un autre tutoriel que j'ai fait pour
`R` disponible
[dans la documentation `utilitr`](https://www.book.utilitr.org/spatdata.html).
Il peut servir de pendant à celui-ci pour l'utilisateur de `R`.
Quelques installations préalables sont nécessaires :
```{python, eval=FALSE, echo=TRUE}
!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.
```{python}
import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt
```
# 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). `R` 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](https://epsg.io/) est un bon aide-mémoire. Les plus fréquents, pour les utilisateurs français, sont les suivants (plus d'infos [ici](https://geodesie.ign.fr/contenu/fichiers/documentation/SRCfrance.pdf)):
* `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*](https://en.wikipedia.org/wiki/Simple_Features) 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
{{% panel status="hint" title="Conseil" icon="fa fa-lightbulb" %}}
Les manipulations de données sur un objet geopandas sont nettement plus lentes que sur un DataFrame traditionnel (car python doit gérer les informations géographiques pendant la manipulation des données). Lorsque vous manipulez des données de grandes dimensions, il peut être préférable d’effectuer les opérations sur les données avant de joindre une géométrie à celles-ci.
{{% /panel %}}
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 ;
* _geojson_ (`.json`) : ce format, non préconisé par l'OGC est largement utilisé pour le développement web comme dans la librairie `Leaflet`. 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.
* _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.
Cette [page](https://si.ecrins-parcnational.com/blog/2020-02-geojson-shapefile-geopackage.html) compare plus en détail ces trois types de formats de données géographiques.
L'aide de [geopandas](https://geopandas.org/io.html) 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. Elles peuvent être trouvées sur le
[site de l'IGN](https://geoservices.ign.fr/adminexpress#telechargementCog) ou récupérées sur
[data.gouv](https://www.data.gouv.fr/fr/datasets/decoupage-administratif-communal-francais-issu-d-openstreetmap/) (construites par `openstreetmap`).
Il est également possible d'utiliser l'une des [API de l'IGN](https://api.gouv.fr/les-api/api_carto_cadastre)
mais ces dernières ne sont pas encore très documentées pour des utilisateurs
de `Python`.
Nous choisirons la deuxième option (import depuis datagouv)
car les premières données sont
compressés au format `.7z`, ce qui n'est pas pratique à dézipper avec
`python` (format propriétaire).
L'inconvénient est que les arrondissements
parisiens ne sont pas présents dans le fichier proposé sur `data.gouv`. Il
faut donc utiliser une source complémentaire, issue de l'*opendata* de la
Mairie de Paris.
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](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 :
```{python download_unzip def}
import requests
import tempfile
import zipfile
url = 'https://www.data.gouv.fr/fr/datasets/r/07b7c9a2-d1e2-4da6-9f20-01a7b72d4b12'
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)
```
```{python apply download_unzip}
download_unzip(url)
communes = gpd.read_file(temporary_location + "/borders/communes-20190101.json")
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. 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é :
```{python}
paris = communes[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
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('map_paris1.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("map_paris1.png")
```
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 :
```{python}
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)
communes = communes[~communes.insee.str.startswith("75")].append(arrondissements)
```
En refaisant la carte ci-dessus, on obtient bien, cette fois, le résultat
espéré
```{python}
paris = communes[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
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('map_paris2.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("map_paris2.png")
```
# 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](https://opendata.paris.fr/explore/dataset/velib-emplacement-des-stations/table/) 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>
```{python}
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)
```
```{python}
communes['dep'] = communes.insee.str[:2]
```
On peut se rassurer en affichant les données sur la carte des communes d'Île-de-France.
```{python, include = FALSE}
fig,ax = plt.subplots(figsize=(10, 10))
stations.sample(200).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2)
communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
color = None)
ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
ax.set_axis_off()
ax
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('map1.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("map1.png")
```
Découvrez ci-dessous par étape les différentes lignes de commandes permettant d'afficher cette carte complète,
étape par étape:
:one:
Afficher le nuage de points de 200 stations vélibs prises au hasard
```{python, include = FALSE}
fig,ax = plt.subplots(figsize=(10, 10))
stations.sample(200).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2)
#communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
# color = None)
#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
# ax.set_axis_off()
ax
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.show()
plt.savefig('map2aa.png', bbox_inches='tight')
```
```{r, echo = FALSE, include = FALSE}
knitr::include_graphics("map2aa.png")
```
:two:
Ajouter à cette couche, en dessous, les contours des communes
```{python, include = FALSE}
fig,ax = plt.subplots(figsize=(10, 10))
stations.sample(200).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2)
communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
color = None)
#ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
# ax.set_axis_off()
ax
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.show()
plt.savefig('map2a.png', bbox_inches='tight')
```
```{r, echo = FALSE, include = FALSE}
knitr::include_graphics("map2a.png")
```
:three:
Ajouter un fond de carte de type *open street map* grâce au package
`contextily`
```{python, include = FALSE}
fig,ax = plt.subplots(figsize=(10, 10))
stations.sample(200).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2)
communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
color = None)
ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Watercolor)
# ax.set_axis_off()
ax
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('map3a.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("map3a.png")
```
:four:
Ne reste plus qu'à retirer l'axe des coordonnées, qui n'est pas très
esthétique. Pour cela:
```{python, include = FALSE, eval = FALSE}
fig,ax = plt.subplots(figsize=(10, 10))
stations.sample(200).plot(ax = ax, color = 'red', alpha = 0.4, zorder=2)
communes[communes['dep'].isin(['75','92','93','94'])].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
color = None)
ctx.add_basemap(ax, crs = stations.crs.to_string(), 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:
```{python}
stations.describe()
```
Pour connaître les plus grands départements:
```{python}
communes.groupby('dep').sum().sort_values('surf_ha', ascending = False)
```
ou les plus grandes communes de France métropolitaine :
```{python}
communes[communes.dep != "97"].sort_values('surf_ha', ascending = False)
```
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`:
```{python}
fig,ax = plt.subplots(figsize=(10, 10))
communes[communes.dep != "97"].dissolve(by='dep', aggfunc='sum').plot(ax = ax, column = "surf_ha")
ax.set_axis_off()
ax
```
```{python, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('dissolve1.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("dissolve1.png")
```
## Opérations sur les géométries
{{% panel status="warning" title="Warning" icon="fa fa-exclamation-triangle" %}}
Les données sont en système de coordonnées WGS 84 ou pseudo-Mercator (epsg: `4326`) et ne sont pas projetées. C'est un format approprié lorsqu'il s'agit d'utiliser un fonds
de carte openstreetmap, stamen, google maps, etc.
Mais ce n'est pas un
format sur lequel on désire faire des calculs car les distances sont faussées sans utiliser de projection. D'ailleurs, geopandas refusera certaines opérations
sur des données dont le crs est `4326`. On reprojete ainsi les données
dans la projection officielle pour la métropole, le Lambert 93
(epsg: 2154).
{{% /panel %}}
```{python, eval = FALSE}
communes = communes.to_crs(2154)
stations = stations.to_crs(2154)
```
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](https://shapely.readthedocs.io/en/latest/manual.html#general-attributes-and-methods).
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$):
```{python, eval = FALSE}
communes['superficie'] = communes.area.div(10**6)
communes
```
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
```{python}
departement = communes[communes.dep == "31"].copy()
departement['geometry'] = departement['geometry'].centroid
ax = departement.plot(figsize = (10,10), color = 'red', alpha = 0.4, zorder=2)
communes[communes['dep'] == "31"].plot(ax = ax, zorder=1, edgecolor = "black", facecolor="none",
color = None)
ctx.add_basemap(ax, crs = stations.crs.to_string(), source = ctx.providers.Stamen.Toner)
ax.set_axis_off()
ax
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('map2.png', bbox_inches='tight')
```
```{r, echo=FALSE}
knitr::include_graphics("map2.png")
```
# 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:
```python
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](https://www.earthdatascience.org/courses/earth-analytics/spatial-data-r/geographic-vs-projected-coordinate-reference-systems-UTM/) 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](https://www.earthdatascience.org/images/courses/earth-analytics/spatial-data/spatial-projection-transformations-crs.png)
La [documentation officielle de geopandas](https://geopandas.org/projections.html) est également très bien
faite sur le sujet. 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.
![](https://imgs.xkcd.com/comics/bad_map_projection_south_america.png)
*Image empruntée à <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`:
```{python}
communes.crs
```
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](https://epsg.io/) 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 !):
```{python, eval = FALSE}
communes = communes.set_crs(2154)
```
Alors que la reprojection (projection Albers: 5070) s'obtient de la manière suivante:
```{python}
fig,ax = plt.subplots(figsize=(10, 10))
communes[communes.dep != "97"].dissolve(by='dep', aggfunc='sum').to_crs(5070).plot(ax = ax)
ax
```
```{python, include = FALSE}
ax.set_axis_off()
ax
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('reproj1.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("reproj1.png")
```
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 !).
```{python, eval = TRUE}
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
```
```{python, echo = FALSE, include = FALSE}
plt.tight_layout(pad=0, h_pad = 0)
plt.savefig('reproj2.png', bbox_inches='tight')
```
```{r, echo = FALSE}
knitr::include_graphics("reproj2.png")
```
# 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:
~~~~python
geopandas_object.merge(pandas_object)
~~~~
Si on utilise deux objets géographiques mais ne désire conserver qu'une seule
dimension géographique[^1], on fera
~~~~python
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*.
[^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.
## Prolongation possible : joindre données sur dimension géographique
{{% panel status="hint" title="Conseil" icon="fa fa-lightbulb" %}}
Les jointures spatiales peuvent être très gourmandes en ressources (car il peut être nécessaire de croiser toutes les géométries de `x` avec toutes les géométries de `y`). Voici deux conseils qui peuvent vous aider:
- Il est préférable de tester les jointures géographiques sur un petit échantillon de données, pour estimer le temps et les ressources nécessaires à la réalisation de la jointure.
- Il est parfois possible d'écrire une fonction qui réduit la taille du problème. Exemple: vous voulez déterminer dans quelle commune se situe un logement dont vous connaissez les coordonnées et le département; vous pouvez écrire une fonction qui réalise pour chaque département une jointure spatiale entre les logements situés dans ce département et les communes de ce département, puis empiler les 101 tables de sorties.
{{% /panel %}}