8 Estymacje jednozmienne

Odtworzenie obliczeń z tego rozdziału wymaga załączenia poniższych pakietów oraz wczytania poniższych danych:

8.1 Kriging

8.1.1 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) (sekcja 8.2)
  • Kriging zwykły (ang. Ordinary kriging) (sekcja 8.3)
  • Kriging z trendem (ang. Kriging with a trend) (sekcja 8.4)
  • Kriging stratyfikowany (ang. Kriging within strata – KWS)
  • Kriging prosty ze zmiennymi średnimi lokalnymi (ang. Simple kriging with varying local means - SKlm) (sekcja 9.1)
  • Kriging z zewnętrznym trendem/Uniwersalny kriging (ang.Kriging with an external trend/Universal kriging) (sekcja 9.2)
  • Kokriging (ang. Co-kriging) (sekcja 10.1)
  • Kriging danych kodowanych (ang. Indicator kriging) (sekcja 11.1)
  • 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 empirycznego, dopasowano model semiwariogramu składający się z funkcji sferycznej o zasięgu 4000 metrów i wartości nuggetu równej 0,5 (rycina 8.1).

##   model psill range
## 1   Nug   0.5     0
## 2   Sph  10.0  4000
Model złożony z modelu nuggetowego i sferycznego dla zmiennej temp.

Rycina 8.1: Model złożony z modelu nuggetowego i sferycznego dla zmiennej temp.

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.

## [1] 15.22251
## [using simple kriging]

Wynik krigingu prostego, jak i każdy inny uzyskany z użyciem pakietu gstat, można podejrzeć wpisując nazwę wynikowego obiektu. 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.

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##    var1.pred        var1.var      
##  Min.   : 8.473   Min.   :0.8211  
##  1st Qu.:12.706   1st Qu.:1.6133  
##  Median :14.961   Median :1.9748  
##  Mean   :15.488   Mean   :2.0411  
##  3rd Qu.:17.712   3rd Qu.:2.3550  
##  Max.   :24.363   Max.   :6.2675  
##  NA's   :1270     NA's   :1270    
## dimension(s):
##   from  to offset delta               refsys point values    
## x    1 127 745542    90 ETRS89 / Poland CS92    NA   NULL [x]
## y    1  96 721256   -90 ETRS89 / Poland CS92    NA   NULL [y]

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

Estymacja i wariancja estymacji używając metody krigingu prostego (SK).

Rycina 8.2: Estymacja i wariancja estymacji używając metody krigingu prostego (SK).

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.
## [using ordinary kriging]

Podobnie jak w przypadku krigingu prostego, można przyjrzeć się wynikom estymacji podając nazwę wynikowego obiektu oraz wyświetlić je używając funkcji plot() (rycina 8.3).

## stars object with 2 dimensions and 2 attributes
## attribute(s):
##    var1.pred        var1.var      
##  Min.   : 8.476   Min.   :0.8215  
##  1st Qu.:12.717   1st Qu.:1.6206  
##  Median :14.967   Median :1.9910  
##  Mean   :15.531   Mean   :2.0820  
##  3rd Qu.:17.709   3rd Qu.:2.3925  
##  Max.   :24.327   Max.   :9.8822  
##  NA's   :1270     NA's   :1270    
## dimension(s):
##   from  to offset delta               refsys point values    
## x    1 127 745542    90 ETRS89 / Poland CS92    NA   NULL [x]
## y    1  96 721256   -90 ETRS89 / Poland CS92    NA   NULL [y]
Estymacja i wariancja estymacji używając metody krigingu zwykłego (OK).

Rycina 8.3: Estymacja i wariancja estymacji używając metody krigingu zwykłego (OK).

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. W pierwszym kroku konieczne jest dodanie współrzęnych do używanego obiektu i do używanej siatki.

Współrzędne dodawane są do całej (regularnej) siatki. Możemy przyciąć je do badanego obszaru poprzez wpisanie wartości NA w miejscach poza naszym obszarem zainteresowań.

Następnie pierwszy z argumentów w funkcji variogram() musi przyjąć postać temp ~ x + y, co oznacza, że uwzględniamy liniowy trend zależny od współrzędnej x oraz y (rycina 8.4).

Semiwariogram zmiennej temp uwzględniający współrzędne x i y.

Rycina 8.4: Semiwariogram zmiennej temp uwzględniający współrzędne x i y.

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 estymacji (rycina 8.5).

##   model     psill    range
## 1   Nug 0.2568762    0.000
## 2   Sph 6.7107464 2087.465
Model semiwariogramu zmiennej temp uwzględniający współrzędne x i y.

Rycina 8.5: Model semiwariogramu zmiennej temp uwzględniający współrzędne x i y.

## [using universal kriging]

Wyświelenie wyników odbywa się używając funkcji plot().

Estymacja i wariancja estymacji używając metody krigingu z trendem (KZT).

Rycina 8.6: Estymacja i wariancja estymacji używając metody krigingu z trendem (KZT).

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 (rycina 8.7). W rozdziałach 10 oraz 9 pokazane będą uzyskane wyniki interpolacji temperatury powietrza korzystając z innych metod krigingu.

Porównanie wyników estymacji używając metody krigingu prostego (SK), krigingu zwykłego (OK) i krigingu z trendem (KZT).

Rycina 8.7: Porównanie wyników estymacji używając metody krigingu prostego (SK), krigingu zwykłego (OK) i krigingu z trendem (KZT).

Porównanie wariancji estymacji używając metody krigingu prostego (SK), krigingu zwykłego (OK) i krigingu z trendem (KZT).

Rycina 8.8: Porównanie wariancji estymacji używając metody krigingu prostego (SK), krigingu zwykłego (OK) i krigingu z trendem (KZT).

8.6 Zadania

Zadania w tym rozdziale są oparte o dane z obiektu punkty_pref. Możesz go wczytać używając poniższego kodu:

Na jego podstawie stwórz trzy obiekty - punkty_pref1 zawierający wszystkie punkty, punkty_pref2 zawierający losowe 100 punktów, oraz punkty_pref3 zawierający losowe 30 punktów.

  1. Zbuduj optymalne modele semiwariogramu zmiennej srtm dla trzech zbiorów danych - punkty_pref1, punkty_pref2, punkty_pref3. Porównaj graficznie uzyskane modele.
  2. W oparciu o uzyskane modele stwórz estymacje zmiennej srtm dla trzech zbiorów danych - punkty_pref1, punkty_pref2, punkty_pref3 używając krigingu prostego. Porównaj graficznie zarówno mapy estymacji jak i mapy wariancji. Opisz zaobserwowane różnice.
  3. W oparciu o uzyskane modele stwórz estymacje zmiennej srtm dla zbioru danych punkty_pref3 używając krigingu zwykłego. Sprawdź jak wygląda wynik estymacji uwzględniając (i) 10 najbliższych obserwacji, (ii) 30 najbliższych obserwacji, (iii) obserwacje w odległości do 2 km.
  4. Używając krigingu z trendem, stwórz optymalne modele zmiennej srtm dla dwóch zbiorów danych - punkty_pref1 oraz punkty_pref3.
  5. Porównaj graficznie zarówno mapy estymacji jak i mapy wariancji dla krigingu prostego, zwykłego oraz z trendem dla danych punkty_pref3. Jakie można zauważyć podobieństwa a jakie różnice?
  6. Dla zmiennej temp z obiektu punkty_pref1 stwórz mapę semiwariogramu. Czy ta zmienna wykazuje anizotropię przestrzenną? Jeżeli tak to stwórz semiwariogramy kierunkowe i ich modele, a następnie estymację zmiennej temp.