-
Notifications
You must be signed in to change notification settings - Fork 0
/
01-ManipulationDonneesSpatiales.qmd
2641 lines (2114 loc) · 122 KB
/
01-ManipulationDonneesSpatiales.qmd
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
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Manipulation des données spatiales dans R {#sec-chap01}
Dans ce chapitre, nous décrivons comment importer, manipuler et cartographier des données spatiales dans R. Pour une description plus détaillée du langage de programmation R -- objets et expression, opérateurs, structures de données (vecteurs, matrices, *arrays*, *DataFrame*), importation et manipulation de données --, lisez le chapitre intitulé [*Prise en main avec R*](https://serieboldr.github.io/MethodesQuantitatives/01-priseenmainR.html) [@RBoldAir].
::: bloc_package
::: bloc_package-header
::: bloc_package-icon
:::
**Liste des *packages* utilisés dans ce chapitre**
:::
::: bloc_package-body
- Pour importer et manipuler des fichiers géographiques :
- `sf` pour importer et manipuler des données vectorielles.
- `rmapshaper` pour simplifier des géométries en conservant la topologie.
- `terra` pour importer et manipuler des données *raster*.
- `gpx` pour importer des coordonnées GPS au format *GPS eXchange Format*.
- `foot` pour créer des enveloppes orientées sur les géométries.
- `concaveman` pour créer des enveloppes concaves.
- Pour cartographier des données :
- `ggplot2` est un *package* pour construire des graphiques qui peut être aussi utilisé pour visualiser des données spatiales.
- `tmap` pour construire des cartes thématiques.
- `RColorBrewer` pour construire une palette de couleur.
- `ggpurb` pour combiner des graphiques et des cartes.
- Pour importer des tables attributaires :
- `foreign` pour importer des fichiers *dBase*.
- `xlsx` pour importer des fichiers Excel.
:::
:::
::: bloc_attention
::: bloc_attention-header
::: bloc_attention-icon
:::
**Pas de panique!**
:::
::: bloc_attention-body
Ce document comprend de nombreuses notions sur l'importation, la manipulation et la cartographie de données spatiales dans R, soit des opérations que vous avez l'habitude de réaliser dans ArcGIS Pro ou QGIS.
Prenez le temps de lire ce premier chapitre à tête reposée et assurez-vous de bien comprendre chaque notion avant de passer à la suivante. L'objectif n'est pas de maîtriser parfaitement la syntaxe R pour toutes les opérations dès la première semaine!
Vous allez manipuler de nombreuses données spatiales avec R au fil de la lecture du livre. Par conséquent, n'hésitez pas à revenir sur ce chapitre lorsque nécessaire; considérez-le comme un aide-mémoire.
:::
:::
## Importation de données géographiques {#sec-011}
::: bloc_notes
::: bloc_notes-header
::: bloc_notes-icon
:::
**Quels *packages* choisir pour importer et manipuler des données spatiales?**
:::
::: bloc_notes-body
Pour les données vectorielles, il existe deux principaux *packages* (équivalent d'une librairie dans Python) : `sp` [@PackageSP1; @PackageSP2] et `sf` [@PackageSF]. Puisque le *package* `sp` est progressivement délaissé par R, il est donc fortement conseillé d'utiliser `sf`.
Pour les données *raster*, il est possible d'utiliser les *packages* `raster` [@PackageRaster] et `terra` [@PackageTerra], dont le dernier, plus récent, semblerait plus rapide.
Cette transition de `sp` à `sf` et de `raster` à `terra` est assez récente et encore en cours durant l'écriture de ce livre. Il existe encore de nombreux *packages* basés sur `sp` et `raster`. Il est donc possible que vous ayez à les utiliser, car leur transition n'est peut-être pas encore effectuée. Notez que la façon dont ces anciens *packages* intègrent les données vectorielles et matricielles dans R est très différente de celle des nouveaux *packages*. À titre d'exemple, la fonction `sp::readOGR` lit un fichier *shapefile*, tout comme la fonction `sf::st_read`, mais la première produit un objet de type `SpatialDataFrame`, alors que la seconde produit un `tbl_df`. Dans le premier cas, les géométries et les données sont stockées dans deux éléments séparés, alors que dans le second cas, le `tbl_df` est un `data.frame` avec une colonne contenant les géométries.
Pour les personnes intéressées aux motivations ayant conduit à cette transition, consultez cet excellent [article de blog](https://r-spatial.org/r/2022/04/12/evolution.html). Il existe deux raisons principales : le mainteneur des *packages* `rgdal` et `rgeos` servant de fondation à `raster` et `sp` a pris sa retraite. À cela s'ajoutent le côté « vieille école » de ces *packages* (ayant plus de 20 ans!) et l'apparition de *packages* plus modernes. Il s'agit d'un bon exemple de ce qui peut arriver dans une communauté *open source* et des évolutions constantes de l'environnement R.
En résumé, privilégiez l'utilisation de `sf` et de `terra.`
Il convient d'installer les deux *packages*. Notez que l'installation d'un *package* requiert une connexion Internet, car R accède au répertoire de *packages* *CRAN* pour télécharger le *package* et l'installer sur votre ordinateur. Cette opération est réalisée avec la fonction `install.packages("nom du package").` Notez qu'une fois que le *package* est installé, il est enregistré localement sur votre ordinateur et y reste à moins de le désinstaller avec la fonction `remove.packages("nom du package").`
Autrement dit, il n'est pas nécessaire de les installer à chaque ouverture de R! Pour utiliser les fonctions d'un *package*, vous devez préalablement le charger avec la fonction `library("Nom du package")` (équivalent à la fonction `import` de Python).
Pour plus d'informations sur l'installation et le chargement de *packages*, consultez la [section suivante](https://serieboldr.github.io/MethodesQuantitatives/01-priseenmainR.html#sec-0123) [@RBoldAir].
:::
:::
### Importation de données vectorielles {#sec-0111}
La fonction `st_read` de `sf` permet d'importer une multitude de formats de données géographiques, comme des fichiers *shapefile* (`shp`), *GeoPackage* (`GPKG`), *GeoJSON* (`geojson`), *sqlite* (`sqlite`), *geodatabase d'ESRI* (`FileGDB`), *Geoconcept* (`gxt`), *Keyhole Markup Language* (`kml`), *Geography Markup Language* (`gml`), etc.
#### Importation d'un fichier *shapefile* {#sec-01111}
Le code R ci-dessous permet d'importer des couches géographiques au format *shapefile*. Notez que la fonction `list.files(pattern = ".shp")` renvoie préalablement la liste des couches *shapefile* présentes dans le dossier de travail.
```{r}
#| echo: true
#| message: false
#| warning: false
## Chargement des packages
library("sf")
library("terra")
library("tmap")
library("ggplot2")
library("ggpubr")
library("foreign")
library("xlsx")
library("rmapshaper")
library("RColorBrewer")
## Obtention d'une liste des shapefiles dans le dossier de travail
list.files(path = "data/chap01/shp", pattern = ".shp")
## Importation des shapefiles avec sf
Arrondissements <- st_read("data/chap01/shp/Arrondissements.shp", quiet=TRUE)
InstallationSport <- st_read("data/chap01/shp/Installations_sportives_et_recreatives.shp", quiet=TRUE)
PistesCyclables <- st_read("data/chap01/shp/Pistes_cyclables.shp", quiet=TRUE)
Rues <- st_read("data/chap01/shp/Segments_de_rue.shp", quiet=TRUE)
```
Regardons à présent la structure des couches importées. Pour ce faire, nous utilisons la fonction `head(nom du DataFrame, n=2)`; notez que le paramètre `n` permet de spécifier le nombre des premiers enregistrements à afficher. Les informations suivantes sont ainsi disponibles :
- `6 fields` : six champs attributaires (`TYPE`, `DETAIL`, `NOM`, `SURFACE`, `ECLAIRAGE`, `OBJECTID`).
- `Geometry type POINT` : le type de géométrie est **point**.
- `Bounding box: xmin: -8009681 ymin: 5686891 xmax: -8001939 ymax: 5696536` : les quatre coordonnées définissant l'enveloppe de la couche.
- `Projected CRS: WGS 84 / Pseudo-Mercator` : la projection cartographique. Ici, une projection cartographique utilisée par Google Maps et OpenStreetMap.
- La géométrie est enregistrée dans le champ `geometry`. Pour le premier enregistrement, nous avons la valeur `POINT (-8001939 5686891)`, soit un point avec les coordonnées géographiques (x,y) entre parenthèses.
```{r}
#| echo: true
#| message: false
#| warning: false
head(InstallationSport, n=2) # Visualisation des deux premiers enregistrements
names(InstallationSport) # Noms de champs (colonnes)
View(InstallationSport) # Afficher l'ensemble de la table attributaire
```
Explorons les types de géométries et la projection des autres couches avec le code ci-dessous. En résumé, les types de géométries sont :
- Des géométries simples
- `point` : un seul point.
- `linestring` : une séquence de deux points et plus formant une ligne.
- `polygon` : un seul polygone formé par une séquence de points pouvant contenir un ou plusieurs polygones intérieurs formant des trous.
- Des géométries multiples
- `multipoint` : plusieurs points pour une même observation.
- `multilinestring` : plusieurs lignes pour une même observation.
- `multipolygon` : plusieurs polygones pour une même observation.
- Une collection de géométries (`Geometrycollection`) qui peut contenir différents types de géométries décrites ci-dessus pour une même observation.
```{r}
#| echo: true
#| message: false
#| warning: false
head(PistesCyclables, n=2)
head(Rues, n=2)
head(Arrondissements, n=2)
```
Visualisons quelques couches importées avec `ggplot()`.
```{r}
#| echo: true
#| message: false
#| warning: false
## Arrondissements et rues
ggplot()+ geom_sf(data = Arrondissements, lwd = .8)+
geom_sf(data = Rues, aes(colour = TYPESEGMEN))
## Arrondissements, pistes cyclables et installations sportives
ggplot()+ geom_sf(data = Arrondissements, lwd = .8)+
geom_sf(data = PistesCyclables, aes(colour = NOM), lwd = .5)+
geom_sf(data = InstallationSport)
```
#### Importation d'une couche dans un *GeoPackage* {#sec-01112}
Pour importer une couche stockée dans un *GeoPackage* (GPKG), vous devez spécifier le fichier et la couche avec respectivement les paramètres `dsn` et `layer` de la fonction `st_read`. Le code ci-dessous permet d'importer les secteurs de recensement de la région métropolitaine de recensement de Sherbrooke pour l'année 2021. Notez que la fonction `st_layers(dsn)` permet d'obtenir la liste des couches contenues dans le fichier GPKG, avec le type de géométrie, les nombre d'entités spatiales et de champs, et la projection cartographique pour chacune d'elles.
```{r}
#| echo: true
#| message: false
#| warning: false
## Nom du fichier GPKG
fichierGPKG <- "data/chap01/gpkg/Recen2021Sherbrooke.gpkg"
## Liste des couches dans le GPKG
st_layers(dsn=fichierGPKG, do_count = TRUE)
## Importation d'une couche
SR.RMRSherb <- st_read(dsn = fichierGPKG,
layer = "SherbSR", quiet=TRUE)
## Affichage des deux premiers enregistrements
head(SR.RMRSherb, n=2)
## Visualisation rapide des secteurs avec ggplot
ggplot()+ geom_sf(data = SR.RMRSherb, lwd = .5)
```
#### Importation d'une couche dans une *geodatabase* d'ESRI {#sec-01113}
La logique est la même qu'avec un *GeoPackage*, nous spécifions le chemin de la *geodatabase* et la couche avec les paramètres `dsn` et `layer`.
```{r}
#| echo: true
#| message: false
#| warning: false
AffectDuTerritoire <- st_read(dsn = "data/chap01/geodatabase/Sherbrooke.gdb",
layer = "AffectationsDuTerritoire", quiet=TRUE)
## Visualisation des affectations du sol avec ggplot
ggplot()+ geom_sf(data = AffectDuTerritoire, aes(fill = TYPE), lwd = .2)
```
#### Importation de données GPS {#sec-01114}
En géomatique appliquée, il est fréquent de collecter des données sur le terrain avec un appareil GPS. Les données ainsi collectées peuvent être enregistrées dans différents formats de données dépendamment de l'appareil GPS utilisé : *GPS eXchange Format* (GPX), *Garmin's Flexible and Interoperable Data Transfer* (FIT), *Training Center XML* (TCX), etc.
##### Importation de coordonnées GPS longitude/latitude au format *csv* {#sec-011142}
Une personne ayant collecté des données sur le terrain pourrait aussi vous les transmettre dans un fichier *csv* (fichier texte délimité par des virgules). Il convient d'importer le fichier de coordonnées GPS dans R dans un *DataFrame* (avec la fonction `read.csv`). Une fois importé, nous constatons qu'il comprend trois champs :
- `id` : un champ identifiant avec des valeurs uniques.
- `lon` : longitude.
- `lat` : latitude.
Les points sont projetés en longitude/latitude (WGS84 long/lat, EPSG : 4326).
```{r}
#| echo: true
#| message: false
#| warning: false
## Importation du fichier csv
PointsGPS <- read.csv("data/chap01/gps/pointsGPS.csv")
head(PointsGPS)
```
Pour convertir le *DataFrame* en un objet `sf`, nous utilisons la fonction `st_as_sf` en spécifiant les champs pour les coordonnées et la projection cartographique.
```{r}
#| echo: true
#| message: false
#| warning: false
## Importation du fichier csv
PointsGPS <- st_as_sf(PointsGPS, coords = c("lon","lat"), crs = 4326)
```
Les points ainsi créés sont localisés dans la ville de Sherbrooke.
```{r}
#| echo: true
#| message: false
#| warning: false
## Affichage des points avec le package tmap
tmap_mode("view") ## Mode actif de tmap
tm_shape(PointsGPS)+
tm_dots(size = 0.05, shape = 21, col = "red")
```
##### Importation de coordonnées GPS au format *GPX* {#sec-011141}
Le format [GPX](https://fr.wikipedia.org/wiki/GPX_(format_de_fichier)) est certainement le format de stockage et d'échange de coordonnées GPS le plus utilisé. Les informations géographiques (x,y) et temporelles (date et heure) sont respectivement enregistrées en degrés longitude/latitude (projection WSG) (WGS84, EPSG : 4326) et en temps universel coordonné (UTC, format ISO 8601).
Pour importer un fichier GPX, nous utilisons le *package* `gpx`. S'il n'est pas installé sur votre ordinateur, lancez la commande `install.packages("gpx")` dans la console de R; n'oubliez pas de le charger avec `library("gpx")`! Ensuite, importez le fichier GPX avec la fonction `read_gpx`, enregistrez la trace GPS dans un *DataFrame* et convertissez-la en objet `sf`.
```{r}
#| echo: true
#| message: false
#| warning: false
library("gpx")
## Importation du fichier GPX
TraceGPS <- read_gpx("data/chap01/gps/TraceGPS.gpx")
## Cette trace GPS comprend trois listes : routes, tracks et waypoints
## Les points sont stockés dans tracks
names(TraceGPS)
## Pour visualiser les données, il suffit de lancer la ligne
## ci-dessous (mise en commentaire car le résultat est un peu long...)
# head(TraceGPS)
TraceGPS <- TraceGPS$tracks$`ID1_PA_2021-12-03_TRAJET01.gpx`
## Conversion du DataFrame en objet sf
TraceGPS <- st_as_sf(TraceGPS, coords = c("Longitude","Latitude"), crs = 4326)
## Visualisation des premiers enregistrements
head(TraceGPS, n=2)
```
La trace GPS correspond à un trajet réalisé à vélo à Abidjan (Côte d'Ivoire) le 3 décembre 2021. Cette trace a été obtenue avec une montre Garmin et comprend un point chaque seconde.
```{r}
#| echo: true
#| message: false
#| warning: false
tmap_mode("view")
tm_basemap(leaflet::providers$OpenStreetMap)+
tm_shape(TraceGPS)+
tm_dots(size = 0.001, shape = 21, col = "red")
```
::: bloc_aller_loin
::: bloc_aller_loin-header
::: bloc_aller_loin-icon
:::
**La structure de la classe `sf`**
:::
::: bloc_aller_loin-body
La classe `sf` est composée de trois éléments (@fig-ClasseSf) :
- L'objet **`simple feature geometry (sfg)`** est la géométrie d'une observation. Tel que vu plus haut, elle est une géométrie simple (`point`, `linestring`, `polygon`), multiple (`multipoint`, `multilinestring`, `multipolygon`) ou une collection de géométries différentes (`Geometrycollection`). Pour définir chacune de ces géométries, nous utilisons les méthodes `st_point()`, `st_linestring()`, `st_polygon()`, `st_multipoint()`, `st_multilinestring()`, `st_multipolygon()` et `Geometrycollection()`.
- L'objet **`simple feature column (sfc)`** est simplement une liste de `simple feature geometry (sfg)`. Elle représente la colonne `geometry` d'une couche vectorielle `sf`.
- L'objet **`data.frame`** correspond à la table attributaire.
Une **`simple feature`** correspond ainsi à une observation (ligne) d'un objet `sf`, soit une entité spatiale comprenant l'information sémantique (attributs) et l'information spatiale (géométrie).
{#fig-ClasseSf width="90%" fig-align="center"}
Voyons un exemple concret : créons une couche `sf` comprenant les trois entités spatiales décrites dans la @fig-ClasseSf.
```{r}
#| echo: true
#| message: false
#| warning: false
## Création des géométries : simple feature geometry (sfg)
point1 = st_point(c(-8001939, 5686891))
point2 = st_point(c(-8009681, 5696536))
point3 = st_point(c(-7998695, 5689743))
## Création d'une liste de géométries : simple feature geometry (sfc)
## avec la projection cartographique EPSG 3857
points_sfc = st_sfc(point1, point2, point3, crs = 3857)
## Création de la table attributaire : objet data.frame
table_attr = data.frame(TYPE = c("Aréna", "Aréna","Aréna"),
NOM = c("Aréna Eugène-Lalonde",
"Aréna Philippe-Bergeron",
"Centre Julien-Ducharme"),
OBJECTID = c(1, 2, 3))
## Création de l'objet sf
Arena_sf = st_sf(table_attr, geometry = points_sfc)
## Le résultat est bien identique à celui de la figure ci-dessus
head(Arena_sf)
```
:::
:::
### Importation de données *raster* {#sec-0112}
La fonction `terra::rast` permet d'importer des images de différents formats (GeoTiff, ESRI, ENVI, ERDAS, BIN, GRID, etc.). Nous importons ci-dessous cinq feuillets de [modèles numériques d'altitude (MNA) à l'échelle du 1/20000](https://www.donneesquebec.ca/recherche/dataset/modeles-numeriques-d-altitude-a-l-echelle-de-1-20-000) couvrant la ville de Sherbrooke. La @fig-MNA présente l'un d'entre eux.
```{r}
#| echo: true
#| message: false
#| warning: false
#| label: fig-MNA
#| fig-align: center
#| fig-cap: "Modèle numérique d'élévation au 1/20000 (feuillet f21e05_101)"
#| out-width: 65%
## Liste des fichiers GeoTIFF dans le dossier
list.files(path="data/chap01/raster", pattern = ".tif")
## Importation des fichiers
f21e05_101 <- terra::rast("data/chap01/raster/f21e05_101.tif")
f21e05_201 <- terra::rast("data/chap01/raster/f21e05_201.tif")
f31h08_102 <- terra::rast("data/chap01/raster/f31h08_102.tif")
f31h08_202 <- terra::rast("data/chap01/raster/f31h08_202.tif")
f21e12_101 <- terra::rast("data/chap01/raster/f21e12_101.tif")
## Visualisation des informations sur l'image f21e05_101
f21e05_101
# Visualisation de l'image
terra::plot(f21e05_101,
main="Modèle numérique d'altitude à l’échelle de 1/20 000 (f21e05_101)")
```
## Manipulation de données vectorielles {#sec-012}
::: bloc_objectif
::: bloc_objectif-header
::: bloc_objectif-icon
:::
***Package*** `sf` **et opérations géométriques**
:::
::: bloc_objectif-body
Le *package* `sf` est une librairie extrêmement complète permettant de réaliser une multitude d'opérations géométriques sur des couches vectorielles comme dans un système d'information géographique (SIG). Notre objectif n'est pas de toutes les décrire, mais d'aborder les principales. Au fil de vos projets avec `sf`, vous apprendrez d'autres fonctions. Pour ce faire, n'hésitez pas à consulter :
- Une belle [*Cheatsheet*](https://github.com/rstudio/cheatsheets/blob/main/sf.pdf) sur `sf`. Allez y jeter un œil, cela vaut la peine!
- Sur le [site CRAN](https://cran.r-project.org/web/packages/sf/index.html) de `sf`, vous trouverez plusieurs vignettes explicatives (exemples de code documentés).
- La [documentation complète en PDF](https://cran.r-project.org/web/packages/sf/sf.pdf).
La syntaxe `methods(class = 'sfc')` renvoie la liste des méthodes implémentées dans le *package* `sf`. Pour accéder à l'aide en ligne de l'une d'entre elles, écrivez simplement `?Nom de la fonction` (ex. : `?st_buffer`).
```{r}
#| echo: true
#| message: false
#| warning: false
methods(class = 'sfc')
```
:::
:::
### Fonctions relatives à la projection cartographique {#sec-0121}
Les trois principales fonctions relatives à la projection cartographique des couches vectorielles sont :
- `st_crs(x)` pour connaître la projection géographique d'un objet `sf`.
- `st_transform(x, cr)` pour modifier la projection cartographique.
- `st_is_longlat(x)` pour vérifier si les coordonnées sont en degrés longitude/latitude.
```{r}
#| echo: true
#| message: false
#| warning: false
## Importation d'un shapefile pour la province de Québec
ProvinceQc <- st_read("data/chap01/shp/Quebec.shp", quiet = TRUE)
## La projection est EPSG:3347 - NAD83 / Statistics Canada Lambert,
## soit la projection conique conforme de Lambert
st_crs(ProvinceQc)
## Reprojection de la couche en WGS84 long/lat (EPSG:4326)
ProvinceQc.4326 <- st_transform(ProvinceQc, crs = 4326)
## longitude/latitude?
st_is_longlat(ProvinceQc)
st_is_longlat(ProvinceQc.4326)
```
La @fig-Projections démontre bien que les deux couches sont projetées différemment.
```{r}
#| echo: true
#| message: false
#| warning: false
#| label: fig-Projections
#| fig-cap: "Deux projections cartographiques"
#| out-width: 75%
#| fig-align: center
Map1 <- ggplot()+geom_sf(data = ProvinceQc)+coord_sf(crs = st_crs(ProvinceQc))+
labs(subtitle = "Conique conforme de Lambert (EPSG : 3347)")
Map2 <- ggplot()+geom_sf(data = ProvinceQc.4326)+coord_sf(crs = st_crs(ProvinceQc.4326))+
labs(subtitle = "WGS84 long/lat (EPSG : 4326)")
comp_plot <- ggarrange(Map1, Map2, ncol = 2, nrow = 1)
annotate_figure(comp_plot,
top = text_grob("Province de Québec",
face = "bold", size = 12, just = "center")
)
```
### Fonctions d'opérations géométriques sur une couche {#sec-0122}
Il existe une quinzaine de fonctions d'opérations géométriques sur une couche dans le *package* `sf` dont le résultat renvoie de nouvelles géométries ([voir la documentation suivante](https://r-spatial.github.io/sf/reference/geos_unary.html)). Nous décrivons ici uniquement celles qui nous semblent les plus utilisées :
- `st_bbox(x)` renvoie les coordonnées minimales et maximales des géométries d'un objet `sf`. Pour créer l'enveloppe d'un objet `sf`, il suffit donc d'écrire `st_as_sfc(st_bbox(x))`.
- `st_boundary(x)` renvoie les limites (contours) des géométries d'un objet `sf`.
- `st_convex_hull(x)` crée l'enveloppe convexe des géométries d'un objet `sf`.
- `st_combine(x)` regroupe les géométries d'un objet `sf` en une seule géométrie, sans les réunir ni résoudre les limites internes.
- `st_union(x)` fusionne les géométries d'un objet `sf` en une seule géométrie.
- `st_buffer(x, dist, endCapStyle = c("ROUND", "FLAT", "SQUARE"), joinStyle = c("ROUND", "MITRE", "BEVEL"))` crée des zones tampons d'une distance définie avec le paramètre `dist`. Cette fonction s'applique à des points, à des lignes et à des polygones.
- `st_centroid(x)` crée des points au centre de chaque géométrie d'un objet `sf`. Elle s'applique donc à des lignes et à des polygones.
- `st_point_on_surface(x)` crée un point au centre de chaque polygone d'un objet `sf` .
- `st_simplify(x, dTolerance)` simplifie les contours de géométries (lignes ou polygones) avec une tolérance exprimée en mètres (paramètre `dTolerance`) d'un objet `sf` .
- `st_voronoi(x, bOnlyEdges = TRUE)` crée des polygones de Thiessen, appelés aussi polygones de Voronoï pour des points. Attention, le paramètre `bOnlyEdges = TRUE` renvoie des lignes tandis que `bOnlyEdges = FALSE` renvoie des polygones.
#### Enveloppe et union d'une couche {#sec-01221}
Le code ci-dessous crée une enveloppe (en bleu) et un polygone fusionné (en rouge) pour les arrondissements de la ville de Sherbrooke (@fig-EnveloppeUnion). La couche résultante de l'opération `st_as_sfc(st_bbox(x))` est ainsi l'équivalent des outils `Emprise` de QGIS et `Minimum Bounding Geometry (Geometry Type = Envelope)` d'ArcGIS Pro.
```{r}
#| echo: true
#| message: false
#| warning: false
## Enveloppe sur les arrondissements de la ville de Sherbrooke
Arrond.Enveloppe <- st_as_sfc(st_bbox(Arrondissements))
## Fusionne les géométries en une seule en résolvant les limites internes
Arrond.Union <- st_union(Arrondissements)
```
```{r}
#| echo: true
#| message: false
#| warning: false
#| label: fig-EnveloppeUnion
#| fig-cap: "Enveloppe sur une couche"
#| fig-align: center
#| out-width: 65%
tmap_mode("plot")
tm_shape(Arrond.Enveloppe) + tm_borders(col = "blue", lwd=2)+
tm_shape(Arrond.Union) + tm_borders(col = "red", lwd=2)+
tm_layout(frame = FALSE)+
tm_scale_bar(c(0,5,10))
```
#### Enveloppe orientée {#sec-01222}
La fonction `st_bbox` de `sf` produit des rectangles englobant des géométries qui sont orientées nord-sud. Il est possible de générer des rectangles orientés autour de géométries pour minimiser leur emprise et ainsi mieux représenter l'orientation de la géométrie initiale. Il n'existe pas de fonction dans `sf` pour le faire, mais le *package* `foot` offre une implémentation facile d'utilisation. Notez que `foot` n'est pas déposé sur CRAN et doit être téléchargé depuis *Github* avec la ligne de code ci-dessous.
```{r}
#| echo: true
#| message: false
#| warning: false
#| eval: false
devtools::install_github("wpgp/foot", build_vignettes = FALSE)
```
La couche résultante de l'opération `fs_mbr(x, returnShape = TRUE)` (@fig-EnveloppeOrientee, b) est ainsi l'équivalent des outils `Emprise orientée minimale (OMBB)` de QGIS et `Minimum Bounding Geometry (Geometry Type = Rectangle by area)` d'ArcGIS Pro.
```{r}
#| echo: true
#| message: false
#| warning: false
library(foot)
## Rectangles (enveloppes) orientés
rectangles_oriented <- fs_mbr(Arrondissements, returnShape = TRUE)
rectangles_oriented <- st_as_sf(rectangles_oriented,
crs = st_crs(Arrondissements))
rectangles_oriented$NOM <- Arrondissements$NOM
## Rectangles non orientés (nord-sud)
st_bbox_by_feature = function(x) {
x = st_geometry(x)
f <- function(y) st_as_sfc(st_bbox(y))
do.call("c", lapply(x, f))
}
rectangles <- st_as_sf(st_bbox_by_feature(Arrondissements),
crs = st_crs(Arrondissements))
rectangles$NOM <- Arrondissements$NOM
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-EnveloppeOrientee
#| fig-cap: "Enveloppes classiques et orientées"
#| fig-align: center
#| out-width: 65%
## Cartographie
Carte1 <- tm_shape(st_cast(rectangles,"LINESTRING")) +
tm_lines(col = "NOM", lwd = 4) +
tm_shape(Arrondissements) +
tm_borders('black') +
tm_layout(legend.show = FALSE, frame = FALSE,
title = "a. Enveloppes nord-sud",
title.size = 1.2)
Carte2 <- tm_shape(st_cast(rectangles_oriented,"LINESTRING")) +
tm_lines(col = "NOM", lwd = 4) +
tm_shape(Arrondissements) +
tm_borders('black') +
tm_layout(legend.show = FALSE, frame = FALSE,
title = "b. Enveloppes orientées",
title.size = 1.2)
tmap_arrange(Carte1, Carte2, nrow = 1, ncol = 2)
```
#### Centroïdes et centre de surface {#sec-01223}
Le code ci-dessous extrait les centres géométriques, c'est-à-dire les centroïdes (en bleu) et les points à l'intérieur des polygones (en rouge) pour les arrondissements de la ville de Sherbrooke (@fig-CentroidesPointsInterieur). Ces deux opérations correspondent aux outils `centroïdes` et `Point dans la surface` de QGIS et `Feature to Point (avec l'option Inside)` d'ArcGIS Pro.
```{r}
#| echo: true
#| message: false
#| warning: false
## Centroïdes et points dans les polygones sur les arrondissements
Arrond.centroide <- st_centroid(Arrondissements)
Arrond.pointpoly <- st_point_on_surface(Arrondissements)
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-CentroidesPointsInterieur
#| fig-cap: "Centroïdes et points à l'intérieur des polygones"
#| fig-align: center
#| out-width: 65%
tm_shape(Arrondissements) + tm_borders(col = "black", lwd=2)+
tm_shape(Arrond.centroide) + tm_dots(size=0.3,col="blue")+
tm_shape(Arrond.pointpoly) + tm_dots(size=0.3,col="red")+
tm_layout(frame = FALSE)
```
#### Zone tampon (*buffer*) {#sec-01224}
Une simple ligne de code permet de créer des zones tampons (équivalent des outils `Analyse vectorielle/Tampon` dans QGIS et `Buffer` dans ArcGIS Pro). Une fois les zones créées, utilisez la fonction `st_union` pour fusionner les tampons en un polygone (@fig-ZoneTampon).
```{r}
#| echo: true
#| message: false
#| warning: false
## Zones tampons de 1000 mètres autour des installations sportives et récréatives
InstSports.buffer <- st_buffer(InstallationSport, dist = 1000)
## Si vous souhaitez fusionner les zones tampons, utilisez la fonction st_union
InstSports.bufferUnion <- st_union(InstSports.buffer)
## Zones tampons de 500 mètres autour des lignes
PistesCyclables.buffer <- st_buffer(PistesCyclables, dist = 500)
PistesCyclables.bufferUnion <- st_union(PistesCyclables.buffer)
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-ZoneTampon
#| fig-cap: "Zones tampons"
#| fig-align: center
#| out-width: 65%
tmap_mode("plot")
M1 <- tm_shape(InstSports.buffer) + tm_borders(col = "black", lwd=1) + tm_fill(col="azure2")+
tm_shape(InstallationSport) + tm_dots(size=0.05,col="red")+
tm_layout(frame = FALSE)
M2 <- tm_shape(InstSports.bufferUnion) + tm_borders(col = "black", lwd=1)+tm_fill(col="azure2")+
tm_shape(InstallationSport) + tm_dots(size=0.05,col="red")+
tm_layout(frame = FALSE)
M3 <- tm_shape(PistesCyclables.buffer) + tm_borders(col = "black", lwd=1) + tm_fill(col="azure2")+
tm_shape(PistesCyclables) + tm_lines(lwd=1,col="red")+
tm_layout(frame = FALSE)
M4 <- tm_shape(PistesCyclables.bufferUnion) + tm_borders(col = "black", lwd=1)+tm_fill(col="azure2")+
tm_shape(PistesCyclables) + tm_lines(lwd=1,col="red")+
tm_layout(frame = FALSE)
tmap_arrange(M1, M2, M3, M4, ncol = 2, nrow =2)
```
Notez que pour des polygones, il est possible de créer des polygones intérieurs comme suit : `st_buffer(x, dist = - Valeur)`. Par exemple, le code ci-dessous crée des polygones de 200 mètres autour et à l'intérieur du parc du Mont-Bellevue de la ville de Sherbrooke (@fig-ZoneTampon2).
```{r}
#| echo: true
#| message: false
#| warning: false
## Importation de la couche des aires aménagées de la ville de Sherbrooke
AiresAmenag <- st_read(dsn = "data/chap01/geodatabase/Sherbrooke.gdb",
layer = "AiresAmenagees", quiet = TRUE)
## Sélection du parc du Mont-Bellevue
MontBellevue <- subset(AiresAmenag, NOM == "Parc du Mont-Bellevue")
## Création d'une zone tampon autour du parc
MontBellevue.ZTA500 <- st_buffer(MontBellevue, dist = 200)
## Création d'une zone tampon à l'intérieur du parc
MontBellevue.ZTI500 <- st_buffer(MontBellevue, dist = -200)
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-ZoneTampon2
#| fig-cap: "Zone tampon intérieure et zone tampon extérieure"
#| fig-align: center
#| out-width: 65%
## Visualisation avec tmap
tmap_mode("plot")
tm_shape(MontBellevue.ZTA500) +tm_borders(col = "red", lwd=1)+
tm_shape(MontBellevue) + tm_borders(col = "black", lwd=1) + tm_fill(col="azure2")+
tm_shape(MontBellevue.ZTI500) +tm_borders(col = "blue", lwd=1)+
tm_layout(frame = FALSE)
```
#### Simplification de géométries {#sec-01225}
La simplification ou généralisation d'une couche de lignes ou de polygones permet de supprimer des sommets tout en gardant le même nombre de géométries dans la couche résultante. Cette opération peut être réalisée dans QGIS avec l'outil `simplifier` et dans ArcGIS Pro avec l'outil `Generalize`. Deux raisons principales peuvent motiver le recours à cette opération :
- La réduction de la taille du fichier, surtout si la couche est utilisée pour de la cartographie interactive sur Internet avec des formats vectoriels comme le SVG (*Scalable Vector Graphics*), le KML ou le GeoJSON.
- L'utilisation de la couche à une plus petite échelle cartographique nécessitant la suppression de détails.
Le code suivant permet de simplifier les contours des arrondissements de la ville de Sherbrooke avec des tolérances de 250, 500, 1000 et 2000 mètres. Plus la valeur de la tolérance est élevée, plus les contours sont simplifiés (@fig-Simplification). Notez que l'algorithme de Douglas-Peucker [@douglas1973algorithms] a été implémenté dans la fonction `st_simplify`. Bien qu'intéressant, cet algorithme ne conserve pas les frontières entre les polygones.
```{r}
#| echo: true
#| message: false
#| warning: false
## Simplification des contours avec différentes distances de tolérance
Arrond.simplify250m <- st_simplify(Arrondissements,
preserveTopology = TRUE,
dTolerance = 250)
Arrond.simplify500m <- st_simplify(Arrondissements,
preserveTopology = TRUE,
dTolerance = 500)
Arrond.simplify1000m <- st_simplify(Arrondissements,
preserveTopology = TRUE,
dTolerance = 1000)
Arrond.simplify2000m <- st_simplify(Arrondissements,
preserveTopology = TRUE,
dTolerance = 2000)
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-Simplification
#| fig-cap: "Simplification des contours de géométries"
#| fig-align: center
#| out-width: 65%
tmap_mode("plot")
M1 <- tm_shape(Arrond.simplify250m)+tm_borders(col = "blue", lwd=2)+tm_layout(frame = FALSE)
M2 <- tm_shape(Arrond.simplify500m)+tm_borders(col = "red", lwd=2)+tm_layout(frame = FALSE)
M3 <- tm_shape(Arrond.simplify1000m)+tm_borders(col = "green", lwd=2)+tm_layout(frame = FALSE)
M4 <- tm_shape(Arrond.simplify2000m)+tm_borders(col = "purple1", lwd=2)+tm_layout(frame = FALSE)
tmap_arrange(M1, M2, M3, M4, ncol = 2, nrow =2)
```
Pour remédier au problème des frontières non conservées, utilisez l'algorithme de Visvalingam et Whyatt [-@visvalingam1993line] avec la fonction `ms_simplify` du *package* `rmapshaper` (@fig-SimplificationVisvalingam), tel qu'illustré dans le code ci-dessous. À titre de rappel, pour l'installer et le charger sur votre ordinateur, tapez dans la console : `install.packages("rmapshaper")` et `library("rmapshaper")`. Le paramètre `keep` permet de définir la proportion de points à retenir : plus sa valeur est faible, plus la simplification est importante.
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-SimplificationVisvalingam
#| fig-cap: "Simplification des contours avec l'algorithme de Visvalingam–Whyatt"
#| fig-align: center
#| out-width: 65%
## Algorithme de Visvalingam–Whyatt
Arrond.simplifyV.005 <- rmapshaper::ms_simplify(Arrondissements,
keep = .005, # proportion des points à retenir (entre 0 et 1)
method = "vis", # Algorithme Visvalingam
keep_shapes = TRUE)
tm_shape(Arrond.simplifyV.005)+tm_borders(col="red")+
tm_layout(frame = FALSE)
```
#### Enveloppe convexe (*convex hull*) {#sec-01226}
Le code ci-dessous permet de créer l'enveloppe convexe pour des points (@fig-ConvexHull). Notez que cette fonction peut également s'appliquer à des lignes et à des polygones. Elle correspond aux outils `Enveloppe convexe` de QGIS et `Feature to Point (avec l'option Convex hull)` d'ArcGIS Pro.
```{r}
#| echo: true
#| message: false
#| warning: false
## Enveloppe convexe autour des points GPS
PointsGPS.Convexhull <- st_convex_hull(st_union(PointsGPS))
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-ConvexHull
#| fig-cap: "Enveloppe convexe autour de points"
#| fig-align: center
#| out-width: 65%
ggplot()+geom_sf(data = PointsGPS.Convexhull)+geom_sf(data = PointsGPS)
```
#### Enveloppe concave (*concave hull*) {#sec-01227}
Une extension possible du polygone convexe est le polygone concave qui a une superficie plus réduite. Il n'existe pas de fonction dans `sf` qui l'implémente. Il faut donc installer un *package* supplémentaire, soit `concaveman`.
```{r}
#| echo: true
#| message: false
#| warning: false
library(concaveman)
## Convex hull autour des points GPS
PointsGPS.Concavhull <- concaveman(PointsGPS)
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-ConcaveHull
#| fig-cap: "Enveloppe concave autour de points"
#| fig-align: center
#| out-width: 65%
ggplot()+geom_sf(data = PointsGPS.Concavhull)+geom_sf(data = PointsGPS)
```
Notez que comparativement au polygone convexe (@fig-ConcaveHull), le polygone concave peut auvoir plus d'une solution possible ([lire l'article suivant](https://portailsig.org/content/sur-la-creation-des-enveloppes-concaves-concave-hull-et-les-divers-moyens-d-y-parvenir-forme.html)). Plus spécifiquement, il faut choisir son degré de concavité. Dans la fonction `concaveman::concaveman`, le paramètre `concavity` prend une valeur numérique, qui, si elle tend vers l'infini, produit un polygone convexe (@fig-ConcaveHull2).
```{r}
#| echo: true
#| message: false
#| warning: false
#| label: fig-ConcaveHull2
#| fig-cap: "Enveloppe concave autour de points"
#| fig-align: center
#| out-width: 80%
library(ggpubr)
# test avec plusieurs valeur de concavité
concav_values <- c(1,1.5,3,8)
plots <- lapply(concav_values, function(i){
Concavhull <- concaveman(PointsGPS, concavity = i)
this_plot <- ggplot()+
geom_sf(data = Concavhull)+
geom_sf(data = PointsGPS)+
labs(subtitle = paste0("Concavité : ",i))+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank())
return(this_plot)
})
ggarrange(plotlist = plots)
```
Dans QGIS, il existe plusieurs *plugins* permettant de générer des enveloppes concaves, ainsi qu'une fonction installée de base avec GRASS (`v.concave.hull`).
### Fonctions d'opérations géométriques entre deux couches {#sec-01213}
Les opérations entre deux couches sont bien connues et largement utilisées dans les SIG. Bien entendu, plusieurs fonctions de ce type sont disponibles dans `sf` et renvoient une nouvelle couche géographique `sf` :
- `st_intersection(x, y)` génère l'intersection entre les géométries de deux couches. À ne pas confondre avec la fonction `st_intersects(x, y)` qui permet de construire une requête spatiale.
- `st_union(x, y)` génère l'union entre les géométries de deux couches.
- `st_difference(x, y)` crée une géométrie à partir de `x` qui n'est pas en intersection avec `y.`
- `st_sym_difference(x, y)` crée une géométrie représentant les portions des géométries `x` et `y` qui ne s'intersectent pas.
- `st_crop(x, y, xmin, ymin, xmax, ymax)` extrait les géométries de `x` comprises dans un rectangle.
En guise de comparaison, toutes ces fonctions sont disponibles dans la boîte à outils de traitement de QGIS (dans le groupe `recouvrement de vecteur`) et les outils de la catégorie `Overlay du Geoprocessing` d'ArcGIS Pro. Le code ci-dessous illustre comment réaliser des intersections et des unions entre deux couches polygonales.
```{r}
#| echo: true
#| message: false
#| warning: false
## Importation des deux couches
polysX <- st_read("data/chap01/shp/PolyX.shp", quiet = TRUE)
polysY <- st_read("data/chap01/shp/PolyY.shp", quiet = TRUE)
## Intersection des deux couches
## Les géométries récupèrent les attributs des deux couches
Inter.XY <- st_intersection(polysX, polysY)
head(Inter.XY)
## Intersection entre deux couches préalablement fusionnées :
## Le résutat est une seule géométrie
Inter.XYUnion <- st_intersection(st_union(polysX), st_union(polysY))
## Union des deux couches
Union.XY <- st_union(st_union(polysX), st_union(polysY))
```
```{r}
#| echo: false
#| message: false
#| warning: false
#| eval: true
Zone <- st_buffer(st_union(st_union(polysX), st_union(polysY)), 20)
tmap_mode("plot")
MapZone <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(polysX)+tm_polygons(col = "red", alpha = .5,
border.col ="black", lwd = 1)+
tm_text(text = "X_id", size = 1)+
tm_shape(polysY)+tm_polygons(col = "yellow", alpha = .5,
border.col ="black", lwd = 1)+
tm_text(text = "Y_id", size = .8, col="black")+
tm_layout(main.title = "Polygones X (rouge) et Y (jaune)", main.title.size = .8)
Inter.XY$id <- 1:nrow(Inter.XY)
MapInter <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Inter.XY)+tm_polygons(col = "tomato3", alpha = .5,
border.col ="black", lwd = 1)+
tm_text(text = "id", size = .8)+
tm_layout(main.title = "st_intersection(polysX, polysY)",
main.title.size = .8)
Inter.XY$id <- NULL
MapInterU <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Inter.XYUnion)+tm_polygons(col = "tomato4", alpha = .5,
border.col ="black", lwd = 1)+
tm_layout(main.title = "st_intersection(st_union(polysX), st_union(polysY))",
main.title.size = .8)
MapUnion <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Union.XY)+tm_polygons(col = "tomato4", alpha = .5,
border.col ="black", lwd = 1)+
tm_layout(main.title = "st_union(st_union(polysX), st_union(polysY))",
main.title.size = .8)
tmap_arrange(MapZone, MapInter, MapInterU, MapUnion, ncol=2, nrow=2)
```
La fonction `st_intersection` peut aussi être utilisée comme la méthode `clip` dans un SIG (ArcGIS Pro ou QGIS). En guise d'exemple, dans le code ci-dessous, nous extrayons les points GPS localisés sur le territoire de la ville de Sherbrooke (@fig-Clip).
```{r}
#| echo: true
#| message: false
#| warning: false
#| label: fig-Clip
#| fig-cap: "Fonction `st_intersection()` équivalente à la méthode clip dans un SIG"
#| fig-align: center
#| out-width: 85%
# Nous nous assurons que les deux couches ont la même projection
PointsGPS <- st_transform(PointsGPS, st_crs(Arrond.Union))
# Extraction des points
PointsGPS.Sherb <- st_intersection(PointsGPS, Arrond.Union)
# Visualisation avant et après
Carte1 <- ggplot()+geom_sf(data = Arrond.Union)+geom_sf(data = PointsGPS)+
labs(subtitle = "Avant l'intersection")+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank())
Carte2 <- ggplot()+geom_sf(data = Arrond.Union)+geom_sf(data = PointsGPS.Sherb)+
labs(subtitle = "Après l'intersection")+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank())
ggarrange(Carte1, Carte2, ncol = 2, nrow = 1)
```
Quelques lignes de code suffisent pour générer les différences de superposition entre les géométries de couches géographiques (@fig-Clip2).
```{r}
#| echo: true
#| message: false
#| warning: false
## Différences entre deux couches
Diff.XY <- st_difference(st_union(polysX), st_union(polysY))
Diff.YX <- st_difference(st_union(polysY), st_union(polysX))
Diff.symXY <- st_sym_difference(st_union(polysY), st_union(polysX))
```
```{r}
#| echo: true
#| message: false
#| warning: false
#| label: fig-Clip2
#| fig-cap: "Différences de superposition entre des géométries de différentes couches"
#| fig-align: center
#| out-width: 75%
tmap_mode("plot")
MapInterU <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Inter.XYUnion)+tm_polygons(col = "tomato4", alpha = .5,
border.col ="black", lwd = 1)
MapDiffXY <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Diff.XY)+
tm_polygons(col = "red", alpha = .5,border.col ="black", lwd = 1)+
tm_layout(main.title = "st_difference(st_union(polysX), st_union(polysY))",
main.title.size = .8)
MapDiffYX <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Diff.YX)+
tm_polygons(col = "yellow", alpha = .5,border.col ="black", lwd = 1)+
tm_layout(main.title = "st_difference(st_union(polysY), st_union(polysX))",
main.title.size = .8)
MapSymDiffYX <- tm_shape(Zone)+tm_fill(col = NA, alpha = 0, border.col = NA, lwd = 0)+
tm_shape(Diff.symXY)+
tm_polygons(col = "green", alpha = .5,border.col ="black", lwd = 1)+
tm_layout(main.title = "st_sym_difference(st_union(polysY), st_union(polysX))",
main.title.size = .8)
tmap_arrange(MapZone, MapDiffXY, MapDiffYX, MapSymDiffYX, ncol=2, nrow=2)
```
### Fonctions de mesures géométriques et de récupération des coordonnées géographiques {#sec-0124}
Les principales fonctions de mesures géométriques et de coordonnées géographiques sont :
- `st_area(x)` calcule la superficie des polygones ou des multipolygones d'une couche `sf` .
- `st_length(x)` calcule la longueur des lignes ou des polylignes d'une couche `sf` .
- `st_distance(x, y)` calcule la distance 2D entre deux objets `sf`, exprimée dans le système de coordonnées de référence.
- `st_coordinates(x)` renvoie les coordonnées géographiques de géométries.
Ci-dessous, nous affichons les superficies des quatre arrondissements, puis nous enregistrons les superficies en m^2^ et en km^2^ dans deux nouveaux champs dénommés `SupM2` et `SupKm2.`
```{r}
#| echo: true
#| message: false
#| warning: false
## Superficie des polygones des arrondissements
st_area(Arrondissements)
## Ajout de champs de superficie dans la table attributaire
Arrondissements$SupM2 <- as.numeric(st_area(st_transform(Arrondissements, crs = 2949)))
Arrondissements$SupKm2 <- as.numeric(st_area(st_transform(Arrondissements, crs = 2949)))/1000000
head(Arrondissements, n=2)
```
De manière très semblable, calculons la longueur de géométries étant des lignes ou des multilignes.
```{r}
#| echo: true
#| message: false
#| warning: false
## Longueurs en mètres
PistesCyclables$longMetre <- as.numeric(st_length(st_transform(PistesCyclables, crs = 2949)))
PistesCyclables$longKm <- as.numeric(st_length(st_transform(PistesCyclables, crs = 2949)))/10000