6 Miary relacji przestrzennych

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

6.1 Geostatystyka

6.1.1 Definicja

Geostatystyka jest to zbiór narzędzi statystycznych uwzględniających w analizie danych ich przestrzenną i czasową lokalizację, a opartych o teorię funkcji losowych.

6.1.2 Funkcje

Istnieją cztery główne funkcje geostatystyki:

  1. Identyfikacja i modelowanie struktury przestrzennej/czasowej zjawiska
  2. Estymacja - szacowanie wartości badanej zmiennej w nieopróbowanym miejscu i/lub momencie czasu
  3. Symulacja - generowanie alternatywnych obrazów, które honorują wyniki pomiarów i strukturę przestrzenną/czasową zjawiska
  4. Optymalizacja próbkowania/sieci pomiarowej

Inaczej mówiąc, celem geostatystyki nie musi być tylko interpolacja (estymacja) przestrzenna, ale również zrozumienie zmienności przestrzennej lub czasowej zjawiska, symulowanie wartości, oraz optymalizacja sieci pomiarowej.

6.1.3 Podstawowe etapy

W przypadku estymacji geostatystycznej, zwanej inaczej interpolacją geostatystyczną, pełna ścieżka postępowania składa się z siedmiu elementów:

  1. Zaprojektowanie sposobu (typu) próbkowania oraz organizacji zadań
  2. Zebranie danych, analiza laboratoryjna
  3. Wstępna eksploracja danych, ocena ich jakości
  4. Modelowanie semiwariogramów na podstawie dostępnych danych
  5. Estymacja badanej cechy
  6. Porównanie i ocena modeli
  7. Stworzenie wynikowego produktu i jego dystrybucja

6.1.4 Dane wejściowe

Jedną z najważniejszych ograniczeń stosowania metod geostatystycznych jest spełnienie odpowiednich założeń dotyczących danych wejściowych. Muszą one:

  1. Zawierać wystarczająco dużą liczbę punktów (minimalnie >30, ale zazwyczaj więcej niż 100/150)
  2. Być reprezentatywne
  3. Być niezależne
  4. Być stworzone używając stałej metodyki
  5. Być wystarczająco dokładne

6.1.5 Badane zjawisko

Oprócz ograniczeń dotyczących danych wejściowych, istnieją również założenia dotyczące analizowanej cechy (analizowanego zjawiska):

  1. Przestrzennej ciągłości - przestrzenna korelacja między zmienny w dwóch lokalizacjach zależy tylko od ich odległości (oraz czasem kierunku), lecz nie od tego gdzie są one położone
  2. Stacjonarności - średnia i wariancja są stałe na całym badanym obszarze

6.1.6 Podstawowe symbole

Podstawowe symbole w określaniu przestrzennej zmienności to:

  • \(s\) - wektor współrzędnych
  • \(z(s)\) - badana zmienna jako funkcja położenia - inaczej określany jako ogon (ang. tail)
  • \(h\) - lag - odstęp pomiędzy dwoma lokalizacjami
  • \(z(s+h)\) - wartość badanej zmiennej odległej o odstęp \(h\) - inaczej określany jako głowa (ang. head)

6.2 Przestrzenna kowariancja i korelacja

6.2.1 Miary podobieństwa

Przestrzenna kowariancja, korelacja i semiwariancja to miary określające przestrzenną zmienność analizowanej cechy.

  • Kowariancja i korelacja to miary podobieństwa pomiędzy dwoma zmiennymi
  • Przenosząc to na aspekt przestrzenny, badamy jedną zmienną, ale pomiędzy dwoma punktami odległymi od siebie o pewien dystans (określany jako lag)
  • W efekcie otrzymujemy miarę podobieństwa pomiędzy wartością głowy i ogona

6.2.2 Wykres rozrzutu z przesunięciem

Wykres rozrzutu z przesunięciem pokazuje korelację pomiędzy wartościami analizowanej cechy w pewnych grupach odległości. Taki wykres można stworzyć używając funkcji hscat() z pakietu gstat. Przykładowo, na poniższym wykresie widać wartość cechy temp z kolejnymi przesunięciami - od 0 do 500 metrów, od 500 metrów do 1000 metrów, itd. W pierwszym przedziale wartość cechy temp z przesunięciem wykazuje korelację na poziomie 0,876, a następnie z każdym kolejnym przedziałem (odległością) ta wartość maleje. W przedziale 3500 do 4000 metrów osiąga jedynie 0,128. Pozwala to na stwierdzenie, że cecha temp wykazuje zmienność przestrzenną - podobieństwo jej wartości zmniejsza się wraz z odległością.

6.2.3 Autokowariancja

Podobną informację jak wykres rozrzutu z przesunięciem daje autokowariancja. Pokazuje ona jak mocno przestrzennie powiązane są wartości par obserwacji odległych od siebie o kolejne przedziały. Jej wyliczenie jest możliwe z użyciem funkcji variogram() z pakietu gstat, gdzie definiuje się analizowaną zmienną, zbiór punktowy, oraz ustala się argument covariogram na TRUE.

6.2.4 Autokorelacja

Kolejną miarą przestrzennego podobieństwa jest autokorelacja. Jej wykres (autokorelogram) pokazuje wartość jednej z miar autokorelacji (np. I Morana lub C Geary’ego) w stosunku do odległości. Na poniższym przykładzie, wartość statystyki I Morana jest wyliczana poprzez funkcję correlog() z pakietu pgirmess.

6.3 Semiwariancja

6.3.1 Miara niepodobieństwa

  • Trzecią miarę charakteryzującą relację między obserwacjami odległymi o kolejne odstępy jest semiwariancja
  • Z praktycznego punktu widzenia, semiwariogram jest preferowaną miarą relacji przestrzennej, ponieważ wykazuje tendencję do lepszego wygładzania danych niż funkcja kowariancji
  • Dodatkowo, semiwariogram jest mniej wymagający obliczeniowo
  • Jednocześnie, dla potrzeb interpretacji relacji, kowarancja i korelacja przestrzenna nadaje się nie gorzej niż semiwariancja

6.3.2 Wzór

Zmienność przestrzenna analizowanej cechy może być określona za pomocą semiwariancji. Jest to połowa średniej kwadratu różnicy pomiędzy wartościami badanej zmiennej w dwóch lokalizacjach odległych o wektor \(h\):

\[ \gamma(h) = \frac{1}{2}[(z(s) - z(s+h))]^2 \]

Przykładowo, aby wyliczyć wartość semiwariancji (gamma) pomiędzy dwoma punktami musimy znać wartość pierwszego z nich (w przykładzie jest to ok. 22,46 stopni Celsjusza) oraz drugiego z nich (ok. 16,04 stopni Celsjusza). Korzystając z wzoru na semiwariację otrzymujemy wartość równą ok. 20,60. Znamy również odległość między punktami (ok. 4576,59 metra), więc możemy w uproszczeniu stwierdzić, że dla tej pary punktów odległych o 4577 metry wartość semiwariancji wynosi około 20,6.

## [1] 20.60122

6.3.3 Chmura semiwariogramu

Jeżeli w badanej próbie mamy \(n\) obserwacji oznacza to, że możemy zaobserwować \(\frac{1}{2}n(n-1)\) par obserwacji. Każda z tych par obserwacji daje nam informację o semiwariancji występującej wraz z odległością. Tę semiwariancję można zaprezentować na wykresie zwanym chmurą semiwariogramu. Do jej wyliczenia służy funkcja variogram() z argumentem cloud = TRUE.

Chmura semiwariogramu pozwala także na zauważenie par punktów, których wartości znacząco odstają. Aby zlokalizować te pary punktów, można zastosować funkcję plot() z argumentem digitize = TRUE, a następnie za pomocą kilku kliknięć określić obszar zawierający nietypowe wartości. Po skończonym wyborze należy nacisnąć klawisz Esc.

Następnie można wyświetlić wybrane pary punktów z użyciem funkcji plot().

6.3.4 Charakterystyka struktury przestrzennej

Semiwariogram jest wykresem pokazującym relację pomiędzy odległością a semiwariancją. Inaczej mówiąc, dla kolejnych odstępów (lagów) wartość semiwariancji jest uśredniana i przestawiania w odniesieniu do odległości.

\[ \hat{\gamma}(h) = \frac{1}{2N(h)}\sum_{i=1}^{N(h)}(z(s_i) - z(s_i+h))^2 \]

,gdzie \(N(h)\) oznacza liczbę par punktów w odstępie \(h\).

W oparciu o semiwariogram empiryczny (czyli oparty na danych) możemy następnie dopasować do niego model/e.

6.3.5 Tworzenie semiwariogramu

Stworzenie podstawowego semiwariogramu w pakiecie gstat odbywa się z użyciem funkcji variogram(). Należy w niej zdefiniować analizowaną zmienną (w tym przykładzie temp~1) oraz zbiór punktowy (punkty).

##      np      dist     gamma dir.hor dir.ver   id
## 1   127  204.2244  1.233999       0       0 var1
## 2   271  494.5572  2.802019       0       0 var1
## 3   522  798.7911  3.797612       0       0 var1
## 4   587 1112.7783  5.037545       0       0 var1
## 5   856 1428.7866  5.967276       0       0 var1
## 6   920 1743.7864  7.411276       0       0 var1
## 7  1068 2062.3041  7.514067       0       0 var1
## 8  1106 2378.9333  8.369087       0       0 var1
## 9  1184 2691.5206  9.539870       0       0 var1
## 10 1215 3011.7729  9.636891       0       0 var1
## 11 1362 3324.6705  9.461429       0       0 var1
## 12 1379 3642.5616 11.301515       0       0 var1
## 13 1405 3963.6776 11.335589       0       0 var1
## 14 1468 4276.0078 13.866888       0       0 var1
## 15 1482 4598.4144 15.353206       0       0 var1

Do wyświetlenia semiwariogramu służy funkcja plot(). Można również dodać informację o liczbie par punktów, jaka posłużyła do wyliczenia semiwariancji dla kolejnych odstępów poprzez argument plot.numbers = TRUE.

6.3.6 Rules of thumb

Przy ustalaniu parametrów semiwariogramu powinno się stosować do kilku utartych zasad (tzw. rules of thumb):

  • W każdym odstępie powinno się znaleźć co najmniej 30 par punktów
  • Maksymalny zasięg semiwariogramu (ang. cutoff distance) to 1/2 pierwiastka z badanej powierzchni (inne źródła mówią o połowie z przekątnej badanego obszaru/jednej trzeciej)
  • Liczba odstępów powinna nie być mniejsza niż 10
  • Optymalnie maksymalny zasięg semiwariogramu powinien być dłuższy o 10-15% od zasięgu zjawiska
  • Optymalnie odstępy powinny być jak najbliżej siebie i jednocześnie nie być chaotyczne
  • Warto metodą prób i błędów określić optymalne parametry semiwariogramu
  • Należy określić czy zjawisko wykazuje anizotropię przestrzenną

6.3.7 Obliczenia pomocnicze

  • Liczba par obserwacji:
## [1] 31125

  • Powierzchnia zajmowana przez jedną próbkę:
## [1] 253506.6
  • Średnia odległość między punktami :
## [1] 503.4944
  • Połowa pierwiastka powierzchni:
## [1] 3980.472

6.3.8 Modyfikacja semiwariogramu

Maksymalny zasięg semiwariogramu (ang. cutoff distance) jest domyślnie wyliczany w pakiecie gstat jako 1/3 z najdłuższej przekątnej badanego obszaru. Można jednak tę wartość zmodyfikować używając argumentu cutoff.

Dodatkowo, domyślnie w pakiecie gstat odległość między przedziałami (ang. interval width) jest wyliczana jako maksymalny zasięg semiwariogramu podzielony przez 15. Można to oczywiście zmienić używając zarówno argumentu cutoff, jak i argumentu width mówiącego o szerokości odstępów.

6.4 Anizotropia

6.4.1 Anizotropia struktury przestrzennej

W wielu sytuacjach, wartość cechy zależy nie tylko od odległości, ale także od kierunku. O takim zjawisku mówimy, że wykazuje ono anizotropię struktury przestrzennej.

6.4.2 Mapa semiwariogramu

Mapa semiwariogramu (zwana inaczej powierzchnią semiwariogramu) służy do określenia czy struktura przestrzenna zjawiska posiada anizotropię, a jeżeli tak to w jakim kierunku. Na podstawie mapy semiwariogramu identyfikuje się także parametry potrzebne do zbudowania semiwariogramów kierunkowych.

Stworzenie mapy semiwariogramu odbywa się poprzez dodanie kilku argumentów do funkcji variogram(): oprócz argumentu zmiennej i zbioru punktowego, jest to cutoff, width, oraz map = TRUE. Następnie można ją wyświetlić używając funkcji plot(). Warto w tym wypadku użyć dodatkowego argumentu threshold, który powoduje, że niepewna wartość semiwariancji (wyliczona na małej liczbie punktów) nie jest wyświetlana.

6.4.3 Semiwariogramy kierunkowe

W przypadku, gdy zjawisko wykazuje anizotropię przestrzenną, możliwe jest stworzenie semiwariogramów dla różnych wybranych kierunków. Dla argumentu alpha = c(0, 45, 90, 135) wybrane cztery główne kierunki to 0, 45, 90 i 135 stopni. Oznacza to, że przykładowo dla kierunku 45 stopni brane pod uwagę będą wszystkie pary punktów pomiędzy 22,5 a 67,5 stopnia.

6.5 Zadania

Przyjrzyj się danym z obiektu punkty_pref. Możesz go wczytać używając poniższego kodu:

  1. Stwórz wykres rozrzutu z przesunięciem dla zmiennej srtm dla odstępów od 0 do 5000 metrów co 625 metry. Co można odczytać z otrzymanego wykresu?
  2. Wylicz odległość oraz wartość semiwariancji dla zmiennej srtm pomiędzy pierwszym i drugim, pierwszym i trzecim, oraz drugim i trzecim punktem. Zwizualizuj te trzy punkty. W jaki sposób można zinterpretować otrzymane wyniki wartości semiwariancji?

  3. Twórz chmurę semiwariogramu dla zmiennej temp. W jaki sposób można zaobserwować na niej wartości odstające? Co one oznaczają?
  4. Jaka jest liczba par obserwacji, średnia powierzchnia zajmowana przez jedną próbkę (w km2) oraz średnia odległość między punktami (w km)?
  5. Stwórz semiwariogram zmiennej srtm. Ile par punktów posłużyło do wyliczenia pierwszego odstępu?
  6. Zmodyfikuj powyższy semiwariogram, żeby jego maksymalny zasięg wynosił 3,5 kilometra.
  7. Stwórz mapę semiwariogramu dla zmiennej srtm. Sprawdź różne wartości cutoff (od 3 do 5 km) oraz width (od 200 do 500 metrów). Zmień domyślną paletę kolorów w wizualizacji mapy semiwariogramu. (Dodatkowe: spróbuj do tego celu użyć pakietu viridisLite).

  8. Co przedstawa stworzona mapa semiwariogramu? Czy badane zjawisko wykazuje izotropia czy anizotropia?
  9. Jeżeli mapa semiwariogramu wykazuje anizotropię struktury przestrzennej to w jakim kierunku? Stwórz semiwariogramy dla wybranych kierunków.