-
Notifications
You must be signed in to change notification settings - Fork 3
/
08-les-reprojections.Rmd
176 lines (116 loc) · 4.99 KB
/
08-les-reprojections.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
# Les reprojections {#les-reprojections}
Dans ce chapitre, nous allons utiliser les librairies suivantes.
```{r pkg_chap_8_reproj}
# CRAN
library(cowplot)
library(mapview)
library(sf)
library(tidyverse)
library(tmap)
```
Pour rappel, il existe deux types de CRS, les CRS géographiques (longitude/latitude avec pour unité de compte des degrés) et les CRS projetés (avec un datum et une unité en mètre par exemple).
Autrefois, la plupart des fonctions de `sf` présupposaient de travailler avec un CRS projeté, car les fonctions de GEOS sur lesquelles elles se basent, le nécessitaient aussi.
## Un premier exemple de reprojection
Prenons les coordonnées de Nantes en WGS 84:
```{r}
nantes <- data.frame(lon = -1.553621, lat = 47.218371) %>%
st_as_sf(coords = c("lon", "lat")) %>%
st_set_crs(4326)
```
On peut visualiser nos données
```{r}
mapview(nantes)
```
`st_is_longlat()` est une fonction de `sf` qui permet de faire un test sur la famille de CRS à laquelle on a affaire.
```{r}
st_is_longlat(nantes)
```
Essayons de créer un buffer de 1 km autour de Nantes.
```{r, warning=TRUE}
nantes_buffer <- st_buffer(nantes, dist = 1000)
mapview(list(nantes, nantes_buffer))
```
Le buffer est pixelisé en longitude / latitude.
Tentons une reprojection. La fonction permettant une reprojection est `st_transform()`.
On va ici passer en Lambert 93 nos données.
```{r}
nantes_proj <- st_transform(nantes, 2154)
```
Le CRS lambert 93 est bien un CRS projeté :
```{r}
st_is_longlat(nantes_proj)
```
```{r, warning=TRUE}
nantes_proj_buffer <- st_buffer(nantes_proj, dist = 1000)
mapview(list(nantes_proj, nantes_proj_buffer))
```
Les contours du buffer sont cette fois plus nets, également les mesures de surfaces sont plus précises en utilisant un crs projeté, adapté à la zone observée.
```{r}
st_area(nantes_buffer)
st_area(nantes_proj_buffer)
pi*1000^2
```
## Quand reprojeter ?
Quelques cas usuels qui peuvent vous amener à reprojeter vos données :
- la manipulation de données fournies dans des CRS différents
- l'usage du package leaflet impose des données spécifiées en WGS 84
- le besoin de visualiser vos données suivant la conversion de certaines propriétés des objets à la surface de la terre.
- l'usage de fonctions recommandant d'utiliser des CRS projetés (comme `st_buffer()`, `st_area()` ci-dessus)
Un exemple d'usage : la distance de Rennes à Nantes
Prenons les coordonnées WGS 84 de Rennes
```{r}
rennes <- data.frame(lon = -1.6777926, lat = 48.117266) %>%
st_as_sf(coords = c("lon", "lat")) %>%
st_set_crs(4326)
```
```{r}
mapview(rennes)
```
Tentons de calculer la distance de Rennes à Nantes.
Avec la données en Lambert 93, la fonction `st_distance()` renvoie un message d'erreur.
```{r, error = TRUE}
st_distance(rennes, nantes_proj)
```
Avec la données en WGS 84, la fonction `st_distance()` renvoie bien un résultat.
```{r, error = TRUE}
st_distance(rennes, nantes)
```
Avec la données en Lambert 93, la fonction `st_distance()` renvoie un résultat plus précis :
```{r}
st_distance(rennes %>% st_transform(2154), nantes_proj)
```
## Quel CRS utiliser ?
A cette question, il y a rarement une seule bonne réponse.
En ce qui concerne les CRS géographiques, le plus simple est d'utiliser le WGS 84 (EPSG 4326), qui est de loin le plus populaire, avec lequel beaucoup de données sont fournies.
En ce qui concerne les CRS projetés, il existe des CRS officiels à utiliser préférentiellement pour les données françaises.
En métropole, on utilise le Lambert 93 (EPSG 2154).
Pour les DROM, on utilise les CRS :
* RGAF09 / UTM zone 20N, code EPSG 5490, pour les Antilles françaises,
* RGFG95 / UTM zone 22N, code EPSG 2972, pour la Guyane,
* RGR92 / UTM zone 40S, code EPSG 2975, pour la Réunion,
* RGM04 / UTM zone 38S, code EPSG 4471, pour Mayotte.
Ensuite votre choix va dépendre des propriétés que vous souhaitez conserver.
## Comment projeter ?
### Projeter des vecteurs
Reprojeter des données vecteur se fait à l'aide de la fonction `st_transform()`, que nous avons vue précédemment.
### Modifier la projection d'une carte.
Parfois on souhaite pouvoir aller plus loin dans les reprojections, en adaptant le centre de la projection, pour cela on peut utiliser une `proj4string` ad hoc.
Pour cela, on va modifier l'argument `+proj` de notre crs avec `st_transform`.
Tentons par exemple de reprojeter notre carte du globe en utilisant la projection azimutale équivalente de Lambert, centrée sur Pékin.
```{r}
data("World")
World_pekin <- st_transform(World, crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=116 +lat_0=40")
```
Le paramètre `+proj=laea` permet de redéfinir la projection, les paramètres `+lon_0` et `lat_0` permettent de définir le centre de la projection. `x_0` et `y_0` définissent le centre du plan pour les coordonnées.
Qu'est ce qui a changé entre nos deux cartes ?
```{r}
st_crs(World)
st_crs(World_pekin)
```
```{r, echo=FALSE}
p1 <- ggplot() +
geom_sf(data = World)
p2 <- ggplot() +
geom_sf(data = World_pekin)
plot_grid(p1, p2)
```