-
-
Notifications
You must be signed in to change notification settings - Fork 1
/
_07-ex-ja.Rmd
89 lines (74 loc) · 3.27 KB
/
_07-ex-ja.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. オブジェクト `nz` を WGS84 CRS に変換した `nz_wgs` というオブジェクトを作成しなさい。
- クラス `crs` のオブジェクトを作成し、CRS を調べなさい。
- オブジェクトの範囲への参照について、CRS によってどの単位を使っているか？
- `nz_wgs` から CRS を削除してプロットしなさい。New Zealand の地図のどこがおかしいか？その理由は？
```{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))
# answer: it is fatter in the East-West direction
# because New Zealand is close to the South Pole and meridians converge there
plot(st_geometry(nz_27700))
par(mfrow = c(1, 1))
```
E2. データセット `world` をユニバーサル横メルカトル図法に変換し (`"+proj=tmerc"`)、っ結果をプロットしなさい。
何が変わったか? その理由は?
WGS 84 に戻してプロットしなさい。
なぜ、このオブジェクトはオリジナルと異なるのか?
```{r 07-ex-e2}
# see 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. 連続色ラスタ (`con_raster`) を、最近傍補間法で NAD83 / UTM zone 12N に変換しなさい。
何が変わったか?
それは、結果にどのように影響するか?
```{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. 影鳥ラスタ (`cat_raster`) を、双一次補間法 (biulinear) で WGS 84 に変換しなさい。
何が変わったか?
それは、結果にどのように影響するか?
```{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) -->
<!-- ``` -->