/
11-dane_kodowane.Rmd
196 lines (152 loc) · 8.26 KB
/
11-dane_kodowane.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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
# Estymacja lokalnego rozkładu prawdopodobieństwa {#estymacja-lokalnego-rozkadu-prawdopodobienstwa}
Odtworzenie obliczeń z tego rozdziału wymaga załączenia poniższych pakietów oraz wczytania poniższych danych:
```{r 11-dane-kodowane-1, message=FALSE, warning=FALSE}
library(sf)
library(stars)
library(gstat)
library(tmap)
library(ggplot2)
library(geostatbook)
data(punkty)
data(siatka)
```
```{r 11-dane-kodowane-2, echo=FALSE}
par(mar = c(rep(0, 4)))
```
## Kriging danych kodowanych
### Kriging danych kodowanych (ang. *Indicator kriging*)
Kriging danych kodowanych to metoda krigingu oparta o dane kategoryzowane lub też dane przetworzone z postaci ciągłej do binarnej.
Jest ona zazwyczaj używana jest to oszacowania prawdopodobieństwa przekroczenia zdefiniowanej wartości progowej, może być również używana do szacowania wartości z całego rozkładu.
Wartości danych wykorzystywane do krigingu danych kodowanych są określone jako 0 lub 1, co reprezentuje czy wartość danej zmiennej jest powyżej czy poniżej określonego progu.
<!--http://geostat-course.org/system/files/geostat13_ind.pdf-->
### Wady i zalety krigingu danych kodowanych
Zalety:
- Możliwość zastosowania, gdy nie interesuje nas konkretna wartość, ale znalezienie obszarów o wartości przekraczającej dany poziom.
- Nie jest istotny kształt rozkładu.
Wady:
- Potencjalnie trudne do modelowania semiwariogramy (szczególnie skrajnych przedziałów).
- Czasochłonność/pracochłonność.
## Przykłady krigingu danych kodowanych
### Binaryzacja danych
Pierwszym krokiem w krigingu danych kodowanych jest stworzenie zmiennej binarnej.
Na poniższym przykładzie tworzona jest nowa zmienna `temp_ind`.
Przyjmuje ona wartość `TRUE` (czyli `1`) dla pomiarów temperatury wyższych niż 20 stopni Celsjusza, a dla pomiarów równych i niższych niż 20 stopni Celsjusza jej wartość wynosi `FALSE` (czyli `0`).
```{r 11-dane-kodowane-3 }
summary(punkty$temp)
punkty$temp_ind = punkty$temp > 20
summary(punkty$temp_ind)
```
W przykładzie, próg został wyznaczony arbitralnie.
Istnieje oczywiście szereg innych możliwości wyznaczania progu.
Można wykorzystać wiedzę zewnętrzną (np. toksyczne stężenie analizowanej substancji) lub też posłużyć się wykresem dystrybuanty do określenia istotnej zmiany wartości (rycina \@ref(fig:11-dane-kodowane-4)).
```{r 11-dane-kodowane-4, fig.cap="Dystrybuanta zmiennej temp."}
ggplot(punkty, aes(temp)) + stat_ecdf()
```
### Modelowanie
Tworzenie i modelowanie semiwariogramu empirycznego w krigingu danych kodowanych wygląda tak samo jak, np. w przypadku krigingu zwykłego (ryciny \@ref(fig:11-dane-kodowane-5) i \@ref(fig:11-dane-kodowane-6)).
Korzystając z funkcji `variogram()` tworzony jest semiwariogram empiryczny, używając `vgm()` tworzony jest model "ręczny", który następnie jest dopasowywany z użyciem funkcji `fit.variogram()`.
```{r 11-dane-kodowane-5, fig.cap="Semiwariogram empiryczny binarnej zmiennej."}
vario_ind = variogram(temp_ind ~ 1, locations = punkty)
plot(vario_ind)
```
```{r 11-dane-kodowane-6, fig.cap = "Model semiwariogramu empirycznego binarnej zmiennej."}
model_ind = vgm(model = "Sph", nugget = 0.01)
fitted_ind = fit.variogram(vario_ind, model_ind)
fitted_ind
plot(vario_ind, model = fitted_ind)
```
### Estymacja
Ostatnim etapem jest stworzenie interpolacji geostatystycznej z pomocą funkcji `krige`.
Wymaga ona czterech argumentów - wzoru (`temp_ind ~ 1`), zbioru punktowego (`punkty`), siatki do interpolacji (`siatka`) oraz modelu (`fitted_ind`).
```{r 11-dane-kodowane-7 }
ik = krige(temp_ind ~ 1,
locations = punkty,
newdata = siatka,
model = fitted_ind)
```
W wyniku estymacji otrzymuje się mapę przestawiającą prawdopodobieństwo przekroczenia zadanej wartości (w tym wypadku jest to 20 stopni Celsjusza) oraz uzyskaną wariancję estymacji (rycina \@ref(fig:plotsy2ok)).
```{r plotsy2ok, fig.height=8, fig.cap = "Estymacja i wariancja estymacji używając metody krigingu danych kodowanych (IK)."}
tm_shape(ik) +
tm_raster(col = c("var1.pred", "var1.var"),
style = "cont",
palette = list("-Spectral", "viridis")) +
tm_layout(legend.frame = TRUE)
```
### Tworzenie mapy binarnej
Mapy przestawiające prawdopodobieństwo można też przetworzyć do postaci map binarnych poprzez użycie wybranej wartości progowej.
Rozkład wartości prawdopodobieństwa jest możliwy do zobaczenia, np. używając histogramu (rycina \@ref(fig:11-dane-kodowane-8)).
```{r 11-dane-kodowane-8, fig.cap = "Rozkład wartości estymacji używając metody krigingu danych kodowanych (IK)."}
ik_df = as.data.frame(ik)
ggplot(ik_df, aes(var1.pred)) + geom_histogram()
```
Następnie można stworzyć nową zmienną binarną w oparciu o wartości estymacji.
Poniższy kod tworzy dwie nowe zmienne - `prog1`, gdzie wartość progowa została arbitralnie ustalona na 0.1 oraz `prog2` z wartością progową 0.75.
```{r 11-dane-kodowane-9}
ik$prog1 = ik$var1.pred > 0.1
ik$prog2 = ik$var1.pred > 0.75
```
W wyniku otrzymuje się binarne mapy, dla których stwierdza się obszary powyżej lub poniższej zadanego prawdopodobieństwa (rycina \@ref(fig:plotsy2ik2)).
```{r plotsy2ik2, fig.height=8, message=FALSE, fig.cap = "Binarne mapy dla progu ustalonego na 0.1 oraz 0.75."}
tm_shape(ik) +
tm_raster(col = c("prog1", "prog2")) +
tm_layout(legend.frame = TRUE)
```
### Alternatywne użycie funkcji
Alternatywnie, zamiast tworzenia nowej zmiennej (takiej jak `temp_ind`), można wykorzystać funkcję `I`.
Z jej użyciem można definiować przyjęte progi bezpośrednio do funkcji `variogram` i `krige`.
Na poniższych przykładach w ten sposób ustalono trzy progi - poniżej 20, poniżej 16, oraz poniżej 12 stopni Celsjusza (rycina \@ref(fig:ploty-trzyik)).
```{r 11-dane-kodowane-10, message=FALSE}
vario_ind20 = variogram(I(temp < 20) ~ 1, locations = punkty)
fitted_ind20 = fit.variogram(vario_ind20,
vgm("Sph", nugget = 0.01))
vario_ind16 = variogram(I(temp < 16) ~ 1, locations = punkty)
fitted_ind16 = fit.variogram(vario_ind16,
vgm("Sph", nugget = 0.03))
vario_ind12 = variogram(I(temp < 12) ~ 1, locations = punkty)
fitted_ind12 = fit.variogram(vario_ind12,
vgm("Sph", nugget = 0.03))
```
```{r 11-dane-kodowane-11, message=FALSE}
ik20 = krige(I(temp < 20) ~ 1,
locations = punkty,
newdata = siatka,
model = fitted_ind20,
nmax = 30)
ik16 = krige(I(temp < 16) ~ 1,
locations = punkty,
newdata = siatka,
model = fitted_ind16,
nmax = 30)
ik12 = krige(I(temp < 12) ~ 1,
locations = punkty,
newdata = siatka,
model = fitted_ind12,
nmax = 30)
```
```{r ploty-trzyik, echo=FALSE, fig.height=12, fig.cap = "Estymacje używając metody krigingu danych kodowanych (IK) dla różnych progów prawdopodobieństwa."}
tm1 = tm_shape(ik20["var1.pred"]) +
tm_raster(breaks = seq(0, 1, by = 0.2)) +
tm_layout(main.title = "Prawdopodobieństwo Temp < 20")
tm2 = tm_shape(ik16["var1.pred"]) +
tm_raster(breaks = seq(0, 1, by = 0.2)) +
tm_layout(main.title = "Prawdopodobieństwo Temp < 16")
tm3 = tm_shape(ik12["var1.pred"]) +
tm_raster(breaks = seq(0, 1, by = 0.2)) +
tm_layout(main.title = "Prawdopodobieństwo Temp < 12")
tmap_arrange(tm1, tm2, tm3, ncol = 1)
```
## Zadania {#z9}
Zadania w tym rozdziale są oparte o dane z obiektu `punkty_ndvi`.
```{r 11-dane-kodowane-12}
data(punkty_ndvi)
```
1. Używając obiektu `punkty_ndvi` stwórz nowe zmienne określające czy zmienna `ndvi` <!-- barren rock, sand, snow, water --><!-- sparse vegetation such as shrubs and grasslands or senescing crops --><!-- dense vegetation --> ma wartość:
* poniżej 0.3
* pomiędzy 0.3 a 0.5
* powyżej 0.5
2. Poczytaj na temat tego wskaźnika.
Co mogą oznaczać powyższe przedziały?
3. Korzystając z trzech powyższych przedziałów stwórz mapy prawdopodobieństwa.
W jaki sposób można zinterpretować trzy uzyskane mapy?
4. Używając wybranego progu prawdopodobieństwa, stwórz trzy mapy binarne.
5. (Dodatkowe) połącz trzy mapy binarne w jedną mapę pokazującą uproszczone pokrycie terenu.