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)
  • 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.

##   model psill range
## 1   Nug   0.5     0
## 2   Sph  10.0  4000

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.95036
## [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.

## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min      max
## x 745541.7 756971.7
## y 712616.2 721256.2
## Is projected: TRUE 
## proj4string :
## [+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.   : 8.934   Min.   :0.7755  
##  1st Qu.:13.246   1st Qu.:1.6193  
##  Median :15.264   Median :1.9533  
##  Mean   :15.877   Mean   :2.0225  
##  3rd Qu.:17.900   3rd Qu.:2.3675  
##  Max.   :25.369   Max.   :4.4515

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

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 używając funkcji summary() oraz wyświetlić je używając funkcji spplot().

## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min      max
## x 745541.7 756971.7
## y 712616.2 721256.2
## Is projected: TRUE 
## proj4string :
## [+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.   : 8.958   Min.   :0.7757  
##  1st Qu.:13.243   1st Qu.:1.6266  
##  Median :15.241   Median :1.9691  
##  Mean   :15.884   Mean   :2.0562  
##  3rd Qu.:18.072   3rd Qu.:2.4076  
##  Max.   :25.339   Max.   :4.7311

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.

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     psill    range
## 1   Nug 0.5308661    0.000
## 2   Sph 6.6620797 2705.928

## [using universal kriging]
## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min      max
## x 745541.7 756971.7
## y 712616.2 721256.2
## Is projected: TRUE 
## proj4string :
## [+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.9749  
##  Mean   :15.86   Mean   :2.0481  
##  3rd Qu.:18.10   3rd Qu.:2.3847  
##  Max.   :25.35   Max.   :4.5134

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 9 pokazane będą uzyskane wyniki interpolacji temperatury powietrza korzystając z innych metod krigingu.

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?