8 Estymacje jednozmienne

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

8.1 Kriging

8.1.1 Kriging | Interpolacja geostatystyczna

Kriging (interpolacja geostatystyczna) to grupa metod estymacji zaproponowana w latach 50. przez Daniela Krige. Główna zasada mówi, że prognoza w danej lokalizacji jest kombinacją obokległych obserwacji. Waga nadawana każdej z obserwacji jest zależna od stopnia (przestrzennej) korelacji - stąd też bierze się istotna rola semiwariogramów.

8.1.2 Metod krigingu

Istnieje szereg meteod krigingu, w tym:

  • Kriging prosty (ang. Simple kriging)
  • Kriging zwykły (ang. Ordinary kriging)
  • Kriging z trendem (ang. Kriging with a trend)
  • Kriging danych kodowanych (ang. Indicator kriging)
  • Kriging stratyfikowany (ang. Kriging within strata – KWS)
  • Kriging prosty ze zmiennymi średnimi lokalnymi (ang. Simple kriging with varying local means - SKlm)
  • Kriging z zewnętrznym trendem/Uniwersalny kriging (ang.Kriging with an external trend/Universal kriging)
  • Kokriging (ang. Co-kriging)
  • Inne

8.2 Kriging prosty

8.2.1 Kriging prosty (ang. Simple kriging)

Kriging prosty zakłada, że średnia jest znana i stała na całym obszarze. W poniższym przykładzie po stworzeniu semiwariogramu empritycznego, dopasowano model semiwariogramu składający się z funkcji sferycznej o zasięgu 4000 metrów i wartości nuggetu równej 0,5.

vario <- variogram(temp~1, punkty)
model <- vgm(10, model = "Sph", range = 4000, nugget = 0.5)
model
##   model psill range
## 1   Nug   0.5     0
## 2   Sph  10.0  4000
fitted <- fit.variogram(vario, model)
plot(vario, model = fitted)

Następnie nastąpiła estymacja wartości z użyciem metody kriginu prostego. W funkcji krige() z pakietu gstat, użycie tej metody wymaga ustalenia średniej wartości cechy za pomocą argumentu beta.

sk <- krige(temp~1, punkty, siatka, model = fitted, beta = 15.32)
## [using simple kriging]

Wynik krigingu prostego, jak i każdy inny uzyskany z użyciem pakietu gstat, można podejrzeć używając funkcji summary(). Szczególnie ważne są dwie, nowe zmienne - var1.pred oraz var1.var. Pierwsza z nich oznacza wartość estymowaną dla każdego oczka siatki, druga zaś mówi o wariancji estymacji.

summary(sk)
## 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.083   Min.   :1.062  
##  1st Qu.:13.208   1st Qu.:1.875  
##  Median :15.224   Median :2.183  
##  Mean   :15.862   Mean   :2.253  
##  3rd Qu.:18.081   3rd Qu.:2.574  
##  Max.   :25.218   Max.   :4.615

Obie uzyskane zmienne można wyświetlić z użyciem funkcji spplot().

spplot(sk, "var1.pred")
spplot(sk, "var1.var")

8.3 Kriging zwykły

8.3.1 Kriging zwykły (ang. Ordinary kriging)

W krigingu zwykłym średnia traktowana jest jako wartość nieznana. Metoda ta uwzględnia lokalne fluktuacje średniej poprzez stosowanie ruchomego okna. Parametry ruchomego okna można określić za pomocą jednego z dwóch argumentów:

  • nmax - użyta zostanie określona liczba najbliższych obserwacji
  • maxdist - użyte zostaną jedynie obserwacje w zadanej odległości
# ok <- krige(temp~1, punkty, siatka, model=fitted, nmax=30)
ok <- krige(temp~1, punkty, siatka, model = fitted, maxdist = 1500)
## [using ordinary kriging]

Podobnie jak w przypadku krigingu prostego, można przyjrzeć się wynikom estymacji używając funkcji summary() oraz wyświetlić je używając funkcji spplot().

summary(ok)
## 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.096   Min.   :1.062  
##  1st Qu.:13.248   1st Qu.:1.880  
##  Median :15.225   Median :2.194  
##  Mean   :15.880   Mean   :2.278  
##  3rd Qu.:18.028   3rd Qu.:2.602  
##  Max.   :25.210   Max.   :4.847
spplot(ok, "var1.pred")
spplot(ok, "var1.var")

8.4 Kriging z trendem

8.4.1 Kriging z trendem (ang. Kriging with a trend)

Kriging z trendem, określany również jako kriging z wewnętrznym trendem, do estymacji wykorzystuje (oprócz zmienności wartości wraz z odległością) położenie analizowanych punktów. Na poniższym przykładzie w funkcji variogram() pierwszy z argumentów przyjął postać temp~x+y, co oznacza, że uwzględniamy liniowy trend zależny od współrzędnej x oraz y.

vario_kzt <- variogram(temp~x + y, data = punkty)
plot(vario_kzt)

Dalszym etapem jest dopasowanie modelu semiwariancji, a następnie wyliczenie estymowanych wartości z użyciem funkcji krige(). Należy tutaj pamiętać, aby wzór (w przykładzie temp~x+y) był taki sam podczas budowania semiwariogramu, jak i interpolacji.

model_kzt <- vgm(psill = 5, model = "Sph", range = 2500, nugget = 1)
fitted_kzt <- fit.variogram(vario_kzt, model_kzt)
fitted_kzt
##   model     psill    range
## 1   Nug 0.5308819    0.000
## 2   Sph 6.6620981 2705.967
plot(vario_kzt, fitted_kzt)

kzt <- krige(temp~x + y, punkty, siatka, model = fitted_kzt)
## [using universal kriging]
summary(kzt)
## 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.14   Min.   :0.8141  
##  1st Qu.:13.20   1st Qu.:1.6472  
##  Median :15.20   Median :1.9750  
##  Mean   :15.86   Mean   :2.0481  
##  3rd Qu.:18.10   3rd Qu.:2.3847  
##  Max.   :25.35   Max.   :4.5134
spplot(kzt, "var1.pred")
spplot(kzt, "var1.var")

8.5 Porównanie wyników SK, OK i KZT

Poniższe porównanie krigingu prostego (SK), zwykłego (OK) i z trendem (KZT) wykazuje niewielkie różnice w uzyskanych wynikach. W rozdziałach 10 oraz 11 pokazane będą uzyskane wyniki interpolacji temperatury powietrza korzystając z innych metod krigingu.