-
-
Notifications
You must be signed in to change notification settings - Fork 1
/
_05-ex-ja.Rmd
145 lines (123 loc) · 5.98 KB
/
_05-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
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
```{r 05-ex-e0, message=FALSE}
library(sf)
library(terra)
library(dplyr)
library(spData)
library(spDataLarge)
```
E1. `nz` データセットの簡略版を生成してプロットしなさい。
`ms_simplify()` の `keep` (0.5 から 0.00005 の範囲) と `st_simplify()` の `dTolerance` (100 から 100,000) の値を変えて実験しなさい。
- 各メソッドで結果の形が崩れ始め、New Zealand を認識できなくなるのはどの値からか？
- 発展: `st_simplify()` の結果のジオメトリ型は、`ms_simplify()` のジオメトリ型と何が違うか？また、どのように解決できるか？
```{r 05-ex-e1}
plot(rmapshaper::ms_simplify(st_geometry(nz), keep = 0.5))
plot(rmapshaper::ms_simplify(st_geometry(nz), keep = 0.05))
# Starts to breakdown here at 0.5% of the points:
plot(rmapshaper::ms_simplify(st_geometry(nz), keep = 0.005))
# At this point no further simplification changes the result
plot(rmapshaper::ms_simplify(st_geometry(nz), keep = 0.0005))
plot(rmapshaper::ms_simplify(st_geometry(nz), keep = 0.00005))
plot(st_simplify(st_geometry(nz), dTolerance = 100))
plot(st_simplify(st_geometry(nz), dTolerance = 1000))
# Starts to breakdown at 10 km:
plot(st_simplify(st_geometry(nz), dTolerance = 10000))
plot(st_simplify(st_geometry(nz), dTolerance = 100000))
plot(st_simplify(st_geometry(nz), dTolerance = 100000, preserveTopology = TRUE))
# Problem: st_simplify returns POLYGON and MULTIPOLYGON results, affecting plotting
# Cast into a single geometry type to resolve this
nz_simple_poly = st_simplify(st_geometry(nz), dTolerance = 10000) |>
st_sfc() |>
st_cast("POLYGON")
nz_simple_multipoly = st_simplify(st_geometry(nz), dTolerance = 10000) |>
st_sfc() |>
st_cast("MULTIPOLYGON")
plot(nz_simple_poly)
length(nz_simple_poly)
nrow(nz)
```
E2. 空間データ操作の章の最初の演習で、Canterbury 地方には New Zealand の 101 の高地のうち 70 地点があることがわかった。
`st_buffer()`を使用して、Canterbury から 100 km 以内にある `nz_height` の点はいくつあるか？
```{r 05-ex-e2}
canterbury = nz[nz$Name == "Canterbury", ]
cant_buff = st_buffer(canterbury, 100000)
nz_height_near_cant = nz_height[cant_buff, ]
nrow(nz_height_near_cant) # 75 - 5 more
```
E3. New Zealand の地理的重心を求めなさい。
Canterburyの地理的重心からの距離は？
```{r 05-ex-e3}
cant_cent = st_centroid(canterbury)
nz_centre = st_centroid(st_union(nz))
st_distance(cant_cent, nz_centre) # 234 km
```
E4. ほとんどの世界地図は北を向いている。
オブジェクト `world` のジオメトリの反射 (この章では触れないアフィン変換の 1 つ) によって、南を上にしたワールドマップを作ることができる。
そのコードを書きなさい。
ヒント: この変換には、本章の `rotation()` 関数を使うことができる。
ボーナス: あなたの国の逆さ地図を作ってみなさい。
```{r 05-ex-e4}
rotation = function(a){
r = a * pi / 180 #degrees to radians
matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
}
world_sfc = st_geometry(world)
world_sfc_mirror = world_sfc * rotation(180)
plot(world_sfc)
plot(world_sfc_mirror)
us_states_sfc = st_geometry(us_states)
us_states_sfc_mirror = us_states_sfc * rotation(180)
plot(us_states_sfc)
plot(us_states_sfc_mirror)
```
E5. Section [5.2.6](https://r.geocompx.org/geometry-operations.html#subsetting-and-clipping) のコードを実行しなさい。そのセクションで作成したオブジェクトを参照して、`x` **と** `y` に含まれる `p` の点の部分集合を作成しなさい。
- 基本サブセット演算子を使用する。
- `st_intersection()`\index{べくた@ベクタ!こうさ@交差}で作成した中間オブジェクトを使用する。
```{r 05-ex-e5a, echo=FALSE}
b = st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points
b = st_buffer(b, dist = 1) # convert points to circles
x = b[1]
y = b[2]
bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
```
```{r 05-ex-e5}
p_in_y = p[y]
p_in_xy = p_in_y[x]
x_and_y = st_intersection(x, y)
p[x_and_y]
```
E6. アメリカの州の境界線の長さをメートル単位で計算しなさい。
どの州の境界線が最も長く、どの州の境界線が最も短いか。
ヒント: `st_length` 関数は `LINESTRING` または `MULTILINESTRING` 形状の長さを計算する。
```{r 05-ex-e6}
us_states9311 = st_transform(us_states, "EPSG:9311")
us_states_bor = st_cast(us_states9311, "MULTILINESTRING")
us_states_bor$borders = st_length(us_states_bor)
arrange(us_states_bor, borders)
arrange(us_states_bor, -borders)
```
E7. srtm.tif ファイルを R で読み込みなさい (`srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))`)。
このラスタの解像度は 0.00083 * 0.00083度。
**terra** パッケージで利用可能なすべてのメソッドを使用して、その解像度を 0.01 * 0.01度に変更しなさい。
結果を視覚化しなさい。
リサンプリング方法の結果の違いは何があるか?
```{r 05-ex-e7}
srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
rast_template = rast(ext(srtm), res = 0.01)
srtm_resampl1 = resample(srtm, y = rast_template, method = "bilinear")
srtm_resampl2 = resample(srtm, y = rast_template, method = "near")
srtm_resampl3 = resample(srtm, y = rast_template, method = "cubic")
srtm_resampl4 = resample(srtm, y = rast_template, method = "cubicspline")
srtm_resampl5 = resample(srtm, y = rast_template, method = "lanczos")
srtm_resampl_all = c(srtm_resampl1, srtm_resampl2, srtm_resampl3,
srtm_resampl4, srtm_resampl5)
plot(srtm_resampl_all)
# differences
plot(srtm_resampl_all - srtm_resampl1, range = c(-300, 300))
plot(srtm_resampl_all - srtm_resampl2, range = c(-300, 300))
plot(srtm_resampl_all - srtm_resampl3, range = c(-300, 300))
plot(srtm_resampl_all - srtm_resampl4, range = c(-300, 300))
plot(srtm_resampl_all - srtm_resampl5, range = c(-300, 300))
```