11 Wykorzystanie do estymacji danych uzupełniających

library('gstat')
library('sp')
library('geostatbook')
data(punkty)
data(siatka)

W wielu przypadkach, oprócz konkretnych pomiarów, istnieje również informacja na temat zmienności innych cech na analizowanym obszarze. W sytuacji, gdy dodatkowe zmienne są skorelowane ze zmienną analizowaną można wykorzystać jedną z metod krigingu wykorzystującą dane uzupełniające, tj. kriging stratyfikowany, prosty kriging ze zmiennymi średnimi lokalnymi, czy kriging uniwersalny.

11.1 Kriging stratyfikowany

11.1.1 Kriging stratyfikowany (ang. Kriging within strata)

Kriging stratyfikowany zakłada, że zmienność badanego zjawiska zależy od cechy jakościowej (kategoryzowanej). Przykładowo, wartość badanej zmiennej jest różna w zależności od pokrycia terenu. Kriging stratyfikowany wymaga posiadania danych zmiennej jakościowej (kategoryzowanej) na całym badanym obszarze.

W poniższym przykładzie zmienną jakościową jest uproszczone pokrycie terenu ze zmiennej clc. Przyjmuje ono jedno z trzech wartości. 1 oznacza obszary rolnicze, 2 oznacza obszary leśne, a 4 oznacza wody powierzchniowe.

siatka$clc <- as.factor(siatka$clc)
spplot(siatka, 'clc')

Kriging stratyfikowany polega na niezależnym tworzeniu i modelowaniu semiwariogramów dla każdej z kategorii.

vario_kws1 <- variogram(temp~1, punkty[punkty$clc==1, ])
plot(vario_kws1)

fitted_kws1 <- fit.variogram(vario_kws1, vgm(10, model = 'Sph', range = 4500, nugget = 0.5))
plot(vario_kws1, fitted_kws1)

vario_kws2 <- variogram(temp~1, punkty[punkty$clc==2, ])
plot(vario_kws2)

fitted_kws2 <- fit.variogram(vario_kws2, vgm(5, model = 'Gau', range = 4500, nugget = 0.1))
plot(vario_kws2, fitted_kws2)

vario_kws4 <- variogram(temp~1, punkty[punkty$clc==4, ])
plot(vario_kws4)

fitted_kws4 <- fit.variogram(vario_kws4, vgm(0.5, model = 'Nug', range = 0))
plot(vario_kws4, fitted_kws4)

Następnie dla każdego obszaru przeprowadzona jest niezależna estymacja wartości analizowanej cechy. Należy jedynie wcześniej zadbać, by w siatce nie było elementów NA dotyczących zmiennych jakościowych. W przykładzie tworzona jest nowa siatka (siatka2) nie zawierająca braków wartości dla zmiennej clc.

siatka2 <- siatka[!is.na(siatka$clc), ]
kws1 <- krige(temp~1, punkty[punkty$clc==1, ], siatka2[na.omit(siatka2$clc==1), ], model = fitted_kws1)
## [using ordinary kriging]
spplot(kws1, 'var1.pred')

kws2 <- krige(temp~1, punkty[punkty$clc==2, ], siatka2[na.omit(siatka2$clc==2), ], model = fitted_kws2)
## [using ordinary kriging]
spplot(kws2, 'var1.pred')

kws4 <- krige(temp~1, punkty[punkty$clc==4, ], siatka2[na.omit(siatka2$clc==4), ], model = fitted_kws4)
## [using ordinary kriging]
spplot(kws4, 'var1.pred')

Ostatnim etapem jest połączenie cząstkowych wyników w jeden obiekt klasy SpatialPixelsDataFrame.

kws <- rbind(as.data.frame(kws1), as.data.frame(kws2), as.data.frame(kws4))
coordinates(kws) <- ~x+y
kws <- as(kws, 'SpatialPixelsDataFrame')

Uzyskane w ten sposób wyniki znacząco różnią się od estymacji krigingem prostym czy zwykłym, wykazując odrębność zmienności w poszczególnych kategoriach pokrycia/użytkowania terenu.

spplot(kws, 'var1.pred', sp.layout=punkty)
spplot(kws, 'var1.var', sp.layout=punkty)

11.2 Prosty kriging ze zmiennymi średnimi lokalnymi (LVM)

11.2.1 Prosty kriging ze zmiennymi średnimi lokalnymi (LVM) (ang. Simple kriging with varying local means)

Prosty kriging ze zmiennymi średnimi lokalnymi zamiast znanej (stałej) stacjonarnej średniej wykorzystuje zmienne średnie lokalne uzyskane na podstawie innej informacji. Lokalna średnia może być uzyskana za pomocą wyliczenia regresji liniowej pomiędzy zmienną badaną a zmienną dodatkową. W takiej sytuacji konieczne jest użycie funkcji lm(). W poniższym przykładzie budowany jest model liniowy relacji pomiędzy temperaturą powietrza (temp), a wysokością nad poziomem morza (srtm).

coef <- lm(temp~srtm, punkty)$coef
coef
##  (Intercept)         srtm 
## 17.506469957 -0.007291269

Wykorzystując relację pomiędzy tymi dwoma zmiennymi tworzony jest semiwariogram empiryczny, który następnie jest modelowany.

vario <- variogram(temp~srtm, punkty)
model_sim <- vgm(10, model = 'Sph', range = 4000, nugget = 1)
model_sim
##   model psill range
## 1   Nug     1     0
## 2   Sph    10  4000
fitted_sim <- fit.variogram(vario, model_sim)
fitted_sim
##   model      psill    range
## 1   Nug  0.7705839    0.000
## 2   Sph 13.1155524 5474.628
plot(vario, model=fitted_sim)

Ostatnim krokiem jest estymacja geostatystyczna, w której oprócz czterech podstawowych argumentów, definiujemy także parametr beta. W tym wypadku jest to wypadku obiekt uzyskany na podstawie regresji liniowej.

sk_lvm <- krige(temp~srtm, punkty, siatka, model=fitted_sim, beta = coef)
## [using simple kriging]
summary(sk_lvm)
## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min      max
## x 745586.7 756926.7
## y 712661.2 721211.2
## Is projected: TRUE 
## proj4string :
## [+init=epsg:2180 +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993
## +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs]
## Number of points: 10993
## Grid attributes:
##    cellcentre.offset cellsize cells.dim
## s1          745586.7       90       127
## s2          712661.2       90        96
## Data attributes:
##    var1.pred         var1.var    
##  Min.   : 9.098   Min.   :1.108  
##  1st Qu.:13.215   1st Qu.:1.931  
##  Median :15.237   Median :2.241  
##  Mean   :15.864   Mean   :2.312  
##  3rd Qu.:18.049   3rd Qu.:2.637  
##  Max.   :25.240   Max.   :4.705  
##  NA's   :43       NA's   :43
spplot(sk_lvm, 'var1.pred')
spplot(sk_lvm, 'var1.var')

11.3 Kriging uniwersalny

11.3.1 Kriging uniwersalny (ang. Universal kriging)

Kriging uniwersalny, określany również jako kriging z trendem (ang. Kriging with a trend model) zakłada , że nieznana średnia lokalna zmienia się stopniowo na badanym obszarze. W krigingu uniwersalnym możemy stosować zarówno zmienne jakościowe, jak i ilościowe.

W pierwszym przykładzie, kriging uniwersalny służy stworzeniu semiwariogramu, modelowaniu oraz estymacji temperatury powietrza z użyciem zmiennej pokrycia terenu.

punkty$clc <- as.factor(punkty$clc)
vario_uk1 <- variogram(temp~clc, punkty)
vario_uk1
##      np      dist     gamma dir.hor dir.ver   id
## 1   127  204.2244  1.513128       0       0 var1
## 2   271  494.5572  3.120279       0       0 var1
## 3   522  798.7911  3.706435       0       0 var1
## 4   587 1112.7783  4.494623       0       0 var1
## 5   856 1428.7866  5.451369       0       0 var1
## 6   920 1743.7864  6.054877       0       0 var1
## 7  1068 2062.3041  6.140934       0       0 var1
## 8  1106 2378.9333  6.241295       0       0 var1
## 9  1184 2691.5206  7.665422       0       0 var1
## 10 1215 3011.7729  7.647765       0       0 var1
## 11 1362 3324.6705  7.973022       0       0 var1
## 12 1379 3642.5616  8.374149       0       0 var1
## 13 1405 3963.6776  8.681657       0       0 var1
## 14 1468 4276.0078  9.723868       0       0 var1
## 15 1482 4598.4144 11.789961       0       0 var1
plot(vario_uk1)

model_uk1 <- vgm(8, model = 'Sph', range = 3000, nugget = 1)
vario_fit_uk1 <- fit.variogram(vario_uk1, model=model_uk1)
vario_fit_uk1
##   model    psill    range
## 1   Nug 1.261541    0.000
## 2   Sph 8.472224 4742.211
plot(vario_uk1, vario_fit_uk1)

siatka$clc <- as.factor(siatka$clc)
spplot(siatka, 'clc')

uk1 <- krige(temp~clc, locations = punkty, newdata=siatka, model=vario_fit_uk1)
## [using universal kriging]
spplot(uk1, 'var1.pred')
spplot(uk1, 'var1.var')

W kolejnym przykładzie zastosowane są już dwie zmienne uzupełniające - wartość wskaźnika wegetacji (ndvi) oraz wysokość nad poziomem morza (srtm).

vario_uk2 <- variogram(temp~ndvi+srtm, punkty)
vario_uk2
##      np      dist     gamma dir.hor dir.ver   id
## 1   127  204.2244  1.324959       0       0 var1
## 2   271  494.5572  2.906522       0       0 var1
## 3   522  798.7911  3.964307       0       0 var1
## 4   587 1112.7783  5.058684       0       0 var1
## 5   856 1428.7866  6.088010       0       0 var1
## 6   920 1743.7864  7.331287       0       0 var1
## 7  1068 2062.3041  7.409558       0       0 var1
## 8  1106 2378.9333  8.222238       0       0 var1
## 9  1184 2691.5206  9.404511       0       0 var1
## 10 1215 3011.7729  9.495409       0       0 var1
## 11 1362 3324.6705  9.403714       0       0 var1
## 12 1379 3642.5616 11.182666       0       0 var1
## 13 1405 3963.6776 11.241615       0       0 var1
## 14 1468 4276.0078 13.783975       0       0 var1
## 15 1482 4598.4144 15.536012       0       0 var1
plot(vario_uk2)

model <- vgm(8, model = 'Sph', range = 3000, nugget = 1)
vario_fit_uk2 <- fit.variogram(vario_uk2, model=model)
vario_fit_uk2
##   model      psill    range
## 1   Nug  0.8757491    0.000
## 2   Sph 13.3016445 5789.379
plot(vario_uk2, vario_fit_uk2)

uk2 <- krige(temp~ndvi+srtm, locations = punkty, newdata=siatka, model=vario_fit_uk2)
## [using universal kriging]
spplot(uk2, 'var1.pred')
spplot(uk2, 'var1.var')