Ćwiczenie 2: Geostatystyczna analiza danych

Modele złożone. Anizotropia.

Author

Anna Dmowska

1 Przygotowanie danych

# Wczytanie bibliotek do analiz przestrzennych, statystycznych i wizualizacji
library(sf)
library(tmap)
library(gstat)
library(terra)
library(tidyverse)

# Wczytanie danych surowych
df = read.csv("data/cw_dane.csv")
# Konwersja tabeli na obiekt przestrzenny w układzie UTM (strefa 33N)
punkty = st_as_sf(df, coords = c("Xutm", "Yutm"), crs = "EPSG:32633")
# Usuwamy wiersze, gdzie wartość parametru Ba (Bar) jest nieznana (NA)
punkty = punkty[!is.na(punkty$Ba),]
# odfiltrowanie wartości anomalnych 
punkty = punkty |> filter(!Lp%in%c(434, 557, 690, 731, 973,974,975,976,977))

1.1 Modele złożone struktury przestrzennej

Gdy pojedynczy model matematyczny nie opisuje wystarczająco dobrze zmienności, stosuje się modele złożone (nazywane także modelami zagnieżdżonymi). Tworzenie modelu złożonego polega to na łączeniu (sumowaniu) kilku modeli (np. sferycznego i Gaussa), co pozwala lepiej dopasować krzywą do semiwariogramu empirycznego.

1.1.1 Model prosty

# Obliczenie semiwariogramu empirycznego dla Baru
vario = variogram(Ba ~ 1, locations = punkty, cutoff = 12000, width = 400)

# Definicja i wizualizacja prostego modelu sferycznego
model = vgm(model = "Sph", psill = 1000, range = 1900, nugget = 40)
plot(vario, model)

1.1.2 Tworzenie modelu złożonego

W R model złożony definiuje się za pomocą funkcji vgm(), używając argumentu add.to, który pozwala na dołączenie kolejnej struktury do modelu podstawowego. Dla każdej części składowej modelu złożonego określa się indywidualne parametry, takie jak typ modelu (model), wariancję progową (psill) oraz zasięg (range), a także wspólny efekt samorodkowy (nugget). Po wstępnym określeniu parametrów poszczególnych składowych modelu złożonego można użyć funkcji fit.variogram() aby lepiej dopasować parametry zasięgu oraz wartości progowej.

W modelu złożonym wyróżnia się tzw. zasięgi cząstkowe (ang. partial ranges) oraz odpowiadające im wartości semiwariancji, czyli progi cząstkowe (ang. partial sills). Suma progów cząstkowych wraz z wariancją nuggetową stanowi całkowitą wariancję progową (sill) modelu.

# Tworzenie modelu złożonego: połączenie modelu sferycznego z modelem gausowskim 
# add.to pozwala na sumowanie struktur o różnych zasięgach
model_zl = vgm(model = "Sph", range = 1500, psill = 1000, nugget = 30,
               add.to = vgm(model = "Gau", psill = 600, range = 10000))

# Automatyczne dopasowanie parametrów modelu złożonego do danych empirycznych
fitted_zl = fit.variogram(vario, model_zl)
plot(vario, fitted_zl)

2 Anizotropia struktury przestrzennej

Anizotropia struktury przestrzennej występuje gdy podobieństwo wartości cechy zmienia się w przestrzeni nie tylko wraz z odległością, ale również z kierunkiem. Kierunek ten można zidentyfikować za pomocą tzw. powierzchni semiwariogramu, nazywanej także mapą semiwariogramu (ang. semivariogram surface).

2.1 Mapa semiwariogramu

Mapa semiwariogramu tworzona jest za pomocą funkcji variogram() z pakietu gstat. Utworzenie mapy wymaga podania zmiennej, zbioru punktowego oraz argumentów cutoff i width. Argument map = TRUE oznacza wygenerowanie mapy zamiast semiwariogramu empirycznego

Argument cutoff określa maksymalną odległość, dla której obliczana jest semiwariancja. Za pomocą tego argumentu definiujemy jak duży obszar ma być pokazany na mapie. Np. ustawiając cutoff = 12000 mapa będzie kwadratem o boku 24 km (od -12 do +12 km).
Argument width określa wielkość boku pojedyńczej komórki. Każda komórka na mapie semiwariogramu reprezentuje pewien przedział do którego wpadają pary punktów o określonej odległości od siebie. Dla każdej komórki obliczana jest średnia wartość semiwariancji. Np. width = 1200 oznacza, że komórki będą miały wielkość 1,2 x 1,2 km. Taka mapa będzie bardziej wygładzonadza, ponieważ uśrednia dane z większego obszaru, co ułatwia dostrzeżenie głównego trendu anizotropii.

# Generowanie mapy semiwariogramu w celu identyfikacji kierunków anizotropii
vario_map = variogram(Ba ~ 1, locations = punkty, 
                      cutoff = 12000, width = 1200, map = TRUE)

Wizualizację mapy przeprowadza się funkcją plot(). Dodatkowo można użyć parametr threshold, który ogranicza wyświetlanie wartości wyliczonych na małej liczbie punktów.

# Wyświetlenie mapy semiwariogramu z zastosowaniem palety "Temps"
plot(vario_map, 
     col.regions = hcl.colors(40, palette = "Temps", rev = TRUE))

2.1.1 Semiwariogramy kierunkowe

Gdy zjawisko wykazuje anizotropię przestrzenną, możliwe jest stworzenie semiwariogramów dla różnych wybranych kierunków. Kierunki definiuje się w funkcji variogram() poprzez dodanie argumentu alpha. Przykładowo, alpha = c(0, 45, 90, 135) tworzy wykresy dla tych czterech kątów. Każdy kierunek obejmuje pary punktów znajdujące się w określonym zakresie kątowym (np. dla kąta 45° jest to przedział 22,5°–67,5°).

#Obliczenie semiwariogramów kierunkowych
vario_kier = variogram(Ba ~ 1, 
                       locations = punkty,
                       alpha = c(0, 45, 90, 135),
                       cutoff = 12000, width = 600)
plot(vario_kier, as.table = TRUE)

Zestawienie semiwariogramów kierunkowych na jednym wykresie pozwala określić jaki typ anizotropii występuje w naszych danych. W geostatystyce wyróżnia się dwa główne rodzaje anizotropii:

Anizotropia geometryczna: Występuje wtedy, gdy semiwariogramy wyznaczone dla różnych kierunków mają taki sam kształt i tę samą wariancję progową (próg, ang. sill), ale różnią się zasięgiem (ang. range)
Anizotropia stefowa: Ma miejsce, gdy semiwariogramy dla poszczególnych kierunków różnią się wartością wariancji progowej (ang. sill).

W praktyce często spotyka się kombinacje obu powyższych typów: w takim modelu zarówno zasięg, jak i wariancja progowa zmieniają się wraz z kierunkiem.

Do lepszego zwizualizowania semiwariogramów kierunkowych użyjemy pakietu ggplot2.

ggplot(vario_kier, aes(x = dist, y = gamma, color = as.factor(dir.hor))) +
  geom_point(size = 2) +
  geom_line(size = 1) +
  scale_color_manual(name = "Kierunek", values = c("0" = "blue", "45" = "orange", "90" = "darkgreen", "135" = "red")) + 
  labs(
    title = "Semiwariogramy kierunkowe dla zmiennej Ba",
    subtitle = "Główny kierunek anizotropii: 135",
    x = "Odstęp - h [m]",
    y = "Semiwariancja - gamma",
    color = "Ustawienia:"
  ) +
  theme_bw()

2.1.2 Modelowanie anizotropowe

Uwzględnienie anizotropii w modelu wymaga zdefiniowania trzech parametrów: kąta dominującego kierunku zmienności, zasięgu w tym kierunku oraz proporcji anizotropii (tj. relacji pomiędzy zasięgiem w przeciwległym kierunku do zasięgu w dominującym kierunku). Przy budowie modelu do funkcji vgm() należy dodać argument anis, który przyjmuje dwie wartości:

Pierwsza wartość anis to dominujący kierunek anizotropii.
Druga wartość to proporcja anizotropii (relacja zasięgu prostopadłego do dominującego).

W poniższym przykładzie utworzono prosty model sferyczny. Argument range oznacza tutaj zasięg w dominującym kierunku zmienności. Argument anis = c(135, 0.5) oznacza, że dominującym kierunkiem zmienności to 135 stopni. Zasięg 5000 m i proporcja 0,5 oznaczają zasięg 2500 m w kierunku prostopadłym.

# Definicja modelu anizotropowego (kierunek 135, proporcja 0.5)
model_kier = vgm(psill = 1000, model = "Sph", range = 5000, 
                     nugget = 0.5, anis = c(135, 0.5))

plot(vario_kier, model_kier, as.table = TRUE)

Funkcja fit.variogram() zastosowana do modelowania semiwariogramów kierunkowych dopasowuje wartości zasięgu (range) oraz wariancji progowej (psill). Wartość kierunku anizotropii oraz proporcji anizotropii nie jest automatycznie modyfikowana.

model_kier_fit = fit.variogram(vario_kier, model_kier)
plot(vario_kier, model_kier_fit, as.table = TRUE)

model_kier_fit

  model     psill    range ang1 anis1
1   Nug  812.6082     0.00    0   1.0
2   Sph 1123.7919 36973.69  135   0.5

2.1.3 Model złożony z anizotropią

Modele złożone mogą również uwzględniać anizotropię struktury przestrzennej poprzez dodanie argumentu anis do poszczególnych komponentów modelu. Wymaga to zdefiniowania kąta dominującego kierunku, zasięgu w tym kierunku oraz proporcji anizotropii (relacji zasięgów w kierunkach prostopadłych). Poniżej znajduje się przykład zaawansowanego modelu złożonego, który uwzględnia parametry anizotropii dla obu składowych struktur.

# Tworzenie modelu złożonego z parametrami anizotropii
model_zl_kier = vgm(model = "Sph", range = 1500, psill = 850, nugget = 30, anis = c(135, 0.6),
               add.to = vgm(model = "Gau", psill = 500, range = 8000, anis = c(135, 0.6)))

plot(vario_kier, model_zl_kier, as.table = TRUE)

# Ostateczne dopasowanie modelu złożonego anizotropowego
model_zl_kier_fit = fit.variogram(vario_kier, model_zl_kier)
plot(vario_kier, model_zl_kier_fit)

3 Dodatkowe informacje

Geostatystyka w R: