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 4000fitted <- 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.615Obie 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 obserwacjimaxdist- 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.847spplot(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.967plot(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.5134spplot(kzt, 'var1.pred')
spplot(kzt, 'var1.var')
