-
-
Notifications
You must be signed in to change notification settings - Fork 7
/
_07-ex.Rmd
89 lines (74 loc) · 3.31 KB
/
_07-ex.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
```{r 07-ex-e0, message=FALSE}
library(sf)
library(terra)
library(spData)
```
E1. Créez un nouvel objet appelé `nz_wgs` en transformant l'objet `nz` avec le CRS WGS84.
- Créez un objet de la classe `crs` pour les deux et utilisez-le pour interroger leurs CRS respectifs.
- En se référant au rectangle de délimitation de chaque objet, quelles unités chaque CRS utilise-t-il ?
- Retirez le CRS de `nz_wgs` et tracez le résultat : qu'est-ce qui ne va pas avec cette carte de la Nouvelle-Zélande et pourquoi ?
```{r 07-ex-e1}
st_crs(nz)
nz_wgs = st_transform(nz, "EPSG:4326")
nz_crs = st_crs(nz)
nz_wgs_crs = st_crs(nz_wgs)
nz_crs$epsg
nz_wgs_crs$epsg
st_bbox(nz)
st_bbox(nz_wgs)
nz_wgs_NULL_crs = st_set_crs(nz_wgs, NA)
nz_27700 = st_transform(nz_wgs, "EPSG:27700")
par(mfrow = c(1, 3))
plot(st_geometry(nz))
plot(st_geometry(nz_wgs))
plot(st_geometry(nz_wgs_NULL_crs))
# réponse: il est plus gras dans le sens est-ouest
# parce que la Nouvelle-Zélande est proche du pôle Sud et que les méridiens y convergent.
plot(st_geometry(nz_27700))
par(mfrow = c(1, 1))
```
E2. Transformez le jeu de données `world` en projection transversale de Mercator (`"+proj=tmerc"`) et tracez le résultat.
Qu'est-ce qui a changé et pourquoi ?
Essayez de le retransformer en WGS 84 et tracez le nouvel objet.
Pourquoi le nouvel objet est-il différent de l'objet original ?
```{r 07-ex-e2}
# cf. https://github.com/r-spatial/sf/issues/509
world_tmerc = st_transform(world, "+proj=tmerc")
plot(st_geometry(world_tmerc))
world_4326 = st_transform(world_tmerc, "EPSG:4326")
plot(st_geometry(world_4326))
```
E3. Transformez le raster de valeurs continues (`con_raster`) en NAD83 / zone UTM 12N en utilisant la méthode d'interpolation du plus proche voisin.
Qu'est-ce qui a changé ?
Comment cela influence-t-il les résultats ?
```{r 07-ex-e3}
con_raster = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
con_raster_utm12n = project(con_raster, "EPSG:32612", method = "near")
con_raster_utm12n
plot(con_raster)
plot(con_raster_utm12n)
```
E4. Transformer le raster avec de catégories (`cat_raster`) en WGS 84 en utilisant la méthode d'interpolation bilinéaire.
Qu'est-ce qui a changé ?
Comment cela influence-t-il les résultats ?
```{r 07-ex-e4}
cat_raster = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
cat_raster_wgs84 = project(cat_raster, "EPSG:4326", method = "bilinear")
cat_raster_wgs84
plot(cat_raster)
plot(cat_raster_wgs84)
```
<!--toDo:jn-->
<!--improve/replace/modify the following q-->
<!-- E5. Create your own proj-string. -->
<!-- It should have the Lambert Azimuthal Equal Area (`laea`) projection, the WGS84 ellipsoid, the longitude of projection center of 95 degrees west, the latitude of projection center of 60 degrees north, and its units should be in meters. -->
<!-- Next, subset Canada from the `world` object and transform it into the new projection. -->
<!-- Plot and compare a map before and after the transformation. -->
<!-- ```{r 06-reproj-40} -->
<!-- new_p4s = "+proj=laea +ellps=WGS84 +lon_0=-95 +lat_0=60 +units=m" -->
<!-- canada = dplyr::filter(world, name_long == "Canada") -->
<!-- new_canada = st_transform(canada, new_p4s) -->
<!-- par(mfrow = c(1, 2)) -->
<!-- plot(st_geometry(canada), graticule = TRUE, axes = TRUE) -->
<!-- plot(st_geometry(new_canada), graticule = TRUE, axes = TRUE) -->
<!-- ``` -->