10 Estymacje wielozmienne

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

10.1 Kokriging

10.1.1 Kokriging (ang. co-kriging)

Kokriging pozwala na wykorzystanie dodatkowej zmiennej (ang. auxiliary variable), zwanej inaczej kozmienną (ang. co-variable), która może być użyta do prognozowania wartości badanej zmiennej w nieopróbowanej lokalizacji. Zmienna dodatkowa może być pomierzona w tych samych miejscach, gdzie badana zmienna, jak też w innych niż badana zmienna. Możliwa jest też sytuacja, gdy zmienna dodatkowa jest pomierzona w dwóch powyższych przypadkach. Kokriging wymaga, aby obie zmienne były istotnie ze sobą skorelowane. Najczęściej kokriging jest stosowany w sytuacji, gdy zmienna dodatkowa jest łatwiejsza (tańsza) do pomierzenia niż zmienna główna. W efekcie, uzyskany zbiór danych zawiera informacje o badanej zmiennej oraz gęściej opróbowane informacje o zmiennej dodatkowej. Jeżeli informacje o zmiennej dodatkowej są znane dla całego obszaru wówczas bardziej odpowiednią techniką będzie kriging z zewnętrznym trendem (KED).

10.1.2 Wybór dodatkowej zmiennej

Wybór zmiennej dodatkowej może opierać się na dwóch kryteriach:

  • Teoretycznym
  • Empirycznym

10.2 Krossemiwariogramy

10.2.1 Krossemiwariogramy (ang. crossvariogram)

Metoda kokrigingu opiera się nie o semiwariogram, lecz o krossemiwariogramy. Krossemiwariogram jest to wariancja różnicy pomiędzy dwiema zmiennymi w dwóch lokalizacjach. Wyliczając krossemiwariogram otrzymujemy empiryczne semiwariogramy dla dwóch badanych zmiennych oraz krosswariogram dla kombinacji dwóch zmiennych.

W poniższym przykładzie istnieją dwie zmienne, savi ze zbioru punkty pomierzona w 242 lokalizacjach oraz ndvi ze zbioru punkty_ndvi pomierzona w 992 punktach (ryciny 10.1 i 10.1).

Mapa wartości zmiennej savi.

Rycina 10.1: Mapa wartości zmiennej savi.

Mapa wartości zmiennej ndvi.

Rycina 10.2: Mapa wartości zmiennej ndvi.

Tworzenie krossemiwariogramów odbywa się z użyciem funkcji gstat(). Na początku definiujemy pierwszy obiekt g. Składa się on z obiektu pustego (NULL), nazwy pierwszej zmiennej (nazwa może być dowolna), wzoru (w przykładzie savi ~ 1), oraz pierwszego zbioru punktowego. Następnie do pierwszego obiektu g dodajemy nowe informacje również poprzez funkcję gstat(). Jest to nazwa obiektu (g), nazwa drugiej zmiennej, wzór, oraz drugi zbiór punktowy.

## data:
## SAVI : formula = savi`~`1 ; data dim = 242 x 5
## NDVI : formula = ndvi`~`1 ; data dim = 992 x 1

Z uzyskanego w ten sposób obiektu tworzymy krossemiwariogram (funkcja variogram()), a następnie go wizualizujemy używając funkcji plot() (rycina 10.3).

Krossemiwariogram zmiennych savi i ndvi.

Rycina 10.3: Krossemiwariogram zmiennych savi i ndvi.

10.3 Modelowanie krossemiwariogramów

Modelowanie krossemiwariogramów, podobnie jak ich tworzenie, odbywa się używając funkcji gstat(). Podaje się w niej wcześniejszy obiekt g, model, oraz argument fill.all = TRUE. Ten ostatni parametr powoduje, że model dodawany jest do wszystkich elementów krossemiwariogramu.

## data:
## SAVI : formula = savi`~`1 ; data dim = 242 x 5
## NDVI : formula = ndvi`~`1 ; data dim = 992 x 1
## variograms:
##              model psill range
## SAVI[1]        Nug 0.001     0
## SAVI[2]        Sph 0.006  1200
## NDVI[1]        Nug 0.001     0
## NDVI[2]        Sph 0.006  1200
## SAVI.NDVI[1]   Nug 0.001     0
## SAVI.NDVI[2]   Sph 0.006  1200

W przypadku semiwariogramów funkcja fit.variogram() służyła dopasowaniu parametrów modelu do semiwariogramu empirycznego. Podobną rolę w krossemiwariogramach spełnia funkcja fit.lmc() - dopasowuje ona liniowy model koregionalizacji do semiwariogramów wielozmienych (rycina 10.4). Funkcja fit.lmc() oczekuje co najmniej dwóch elementów, krossemiwariogramu oraz modelów krossemiwariancji. W poniższym przykładzie dodatkowo użyto parametru correct.diagonal = 1.01, z uwagi na to że analizowane zmienne wykazywały bardzo silną korelację oraz parametru fit.method. Ten ostatni parametr określa, która z użytych metod automatycznego dopasowania semiwariogramów jest używana. Pakiet gstat ma domyślnie ustawione fit.method = 7, co daje największe znaczenie parom punktów o najmniejszym zasięgu. Ta opcja nie jest zazwyczaj optymalna w przypadku krossemiwariogramów - tutaj częściej stosuje się fit.method = 6 (wagi nie są nadawane) lub fit.method = 1 (wagi proporcjonalne do liczby par punktów w każdym przedziale).

## data:
## SAVI : formula = savi`~`1 ; data dim = 242 x 5
## NDVI : formula = ndvi`~`1 ; data dim = 992 x 1
## variograms:
##              model        psill range
## SAVI[1]        Nug 0.0009873135     0
## SAVI[2]        Sph 0.0050465345  1200
## NDVI[1]        Nug 0.0023006100     0
## NDVI[2]        Sph 0.0070299532  1200
## SAVI.NDVI[1]   Nug 0.0014922022     0
## SAVI.NDVI[2]   Sph 0.0055018027  1200
Modele krossemiwariogramu zmiennych savi i ndvi.

Rycina 10.4: Modele krossemiwariogramu zmiennych savi i ndvi.

10.4 Kokriging

Posiadając dopasowane modele oraz siatkę można uzyskać wynik używając funkcji predict() (rycina 10.5).

## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]

W efekcie otrzymujemy pięć zmiennych:

  1. SAVI.pred - estymacja zmiennej savi
  2. SAVI.var - wariancja zmiennej savi
  3. NDVI.pred - estymacja zmiennej ndvi
  4. NDVI.var - wariancja zmiennej ndvi
  5. cov.SAVI.NDVI - kowariancja zmiennych savi oraz ndvi
## stars object with 2 dimensions and 5 attributes
## attribute(s):
##    SAVI.pred        SAVI.var         NDVI.pred        NDVI.var      
##  Min.   :0.0772   Min.   :0.0015   Min.   :0.2055   Min.   :0.0031  
##  1st Qu.:0.2909   1st Qu.:0.0024   1st Qu.:0.4686   1st Qu.:0.0039  
##  Median :0.3212   Median :0.0026   Median :0.5081   Median :0.0042  
##  Mean   :0.3157   Mean   :0.0027   Mean   :0.5012   Mean   :0.0043  
##  3rd Qu.:0.3480   3rd Qu.:0.0029   3rd Qu.:0.5421   3rd Qu.:0.0046  
##  Max.   :0.4371   Max.   :0.0048   Max.   :0.6585   Max.   :0.0074  
##  NA's   :1270     NA's   :1270     NA's   :1270     NA's   :1270    
##  cov.SAVI.NDVI   
##  Min.   :0.0019  
##  1st Qu.:0.0026  
##  Median :0.0029  
##  Mean   :0.0029  
##  3rd Qu.:0.0032  
##  Max.   :0.0054  
##  NA's   :1270    
## dimension(s):
##   from  to offset delta                       refsys point values    
## x    1 127 745542    90 +proj=tmerc +lat_0=0 +lon...    NA   NULL [x]
## y    1  96 721256   -90 +proj=tmerc +lat_0=0 +lon...    NA   NULL [y]
Estymacja i wariancja estymacji zmiennej savi używając metody kokrigingu (CK).

Rycina 10.5: Estymacja i wariancja estymacji zmiennej savi używając metody kokrigingu (CK).

10.5 Zadania

Zadania w tym rozdziale są oparte o dane z meuse.all z pakietu gstat.

Na jego podstawie wydziel dwa obiekty - meuse164 zawierający tylko zmienną cadmium (kadm) dla 164 punktów, oraz meuse60 zawierający tylko zmienną copper (miedź) dla 60 punktów.

  1. Stwórz siatkę interpolacyjną o rozdzielczości 100 jednostek dla obszaru, w którym znajdują się punkty meuse.
  2. Zbuduj optymalne modele semiwariogramu zmiennej cadmium dla obiektu meuse164 oraz zmiennej copper dla obiektu meuse60. Porównaj graficznie uzyskane modele.
  3. Korzystając z obiektów meuse164 oraz meuse60 stwórz krossemiwariogram.
  4. Zbuduj ręczny model uzyskanego krossemiwariogramu. Następnie stwórz model automatyczny. Porównaj uzyskane wyniki.
  5. Stwórz estymację zmiennej copper w nowo utworzonej siatce korzystając z kokrigingu.
  6. Dodatkowe: stwórz estymację zmiennej savi z obiektu punkty używając krigingu zwykłego. Porównaj uzyskaną estymację z estymacją przedstawioną w tym rozdziale, stworzoną używając kokrigingu.