Geostatystyczna analiza danych

Wpływ danych odstajacych na wyniki analizy

Author

Anna Dmowska

1 Pakiety i fukcje pomocnicze

# Wczytanie niezbędnych bibliotek
library(sf)         # Dane wektorowe
library(tmap)       # Tworzenie map
library(gstat)      # Geostatystyka (kriging, variogramy)
library(terra)      # Dane rastrowe
library(tidyverse)  # Przetwarzanie danych i wykresy (ggplot2)
library(patchwork)  # Łączenie wykresów

# Funkcja pomocnicza do interpolacji gstat na obiektach SpatRaster
# Pozwala na użycie gstat bezpośrednio na siatce rastrowej z pakietu terra
interpolate_gstat = function(model, x, ...) {
  # Konwersja współrzędnych rastra na punkty (sf) z zachowaniem układu współrzędnych (CRS)
  v = st_as_sf(x, coords = c("x", "y"), crs = st_crs(model$data[[1]]$data))
  # Wykonanie predykcji (krigingu) dla tych punktów
  p = predict(model, v, ...)
  # Usunięcie geometrii, aby zwrócić czyste wartości tabelaryczne do rastra
  st_drop_geometry(p)
}

# Definicja siatki interpolacyjnej (obszar badań, rozdzielczość 100m, układ UTM 33N)
siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150, res = 100, crs = "EPSG:32633")

2 Dane odstające

Przed wykonaniem analizy geostatystycznej należy zidentyfikować, czy w zbiorze danych nie występują globalne lub lokalne wartości odstające.

  • Dane globalnie odstające to bardzo niskie lub bardzo wysokie wartości względem pozostałych w zbiorze danych. Można je zidentyfikować poprzez wykonanie histogramu oraz obejrzenie rozkładu danych.

  • Dane lokalnie odstające, to wartości które odbiegają jedynie od swoich sąsiadów w przestrzeni lub czasie, nie stanowiąc ekstremów dla całego zbioru. Wartości lokalnie odstające identyfikuje się za pomocą mapy lokalizacyjnej.

3 Wpływ wartości odstających

Wartości odstające i ekstremalne mają bardzo silny i zazwyczaj negatywny wpływ na semiwariancję.

  • Semiwariancja obliczana jest jako połowa średniej kwadratów różnic między wartościami. Ze względu na podnoszenie tych różnic do kwadratu, nawet pojedyncze, bardzo wysokie lub niskie wyniki pomiarów drastycznie zawyżają wartość semiwariancji dla danego odstępu.

  • Pojedyncze anomalie mogą zniekształcić wynik analizy opartej na setkach lub tysiącach poprawnych danych. Powoduje to, że wykresy semiwariogramów eksperymentalnych stają się chaotyczne i przestają wiarygodnie odzwierciedlać rzeczywiste relacje przestrzenne w badanej populacji.

  • Negatywne oddziaływanie wartości odstających jest najbardziej znaczące dla pierwszych trzech przedziałów odległości (lagów). Zawierają one zazwyczaj znacznie mniej par punktów niż kolejne odstępy, przez co jedna para z wartością odstającą ma w nich proporcjonalnie większy udział w końcowym wyniku.

3.1 Zmniejszenie wpływu wartości odstających

Aby zminimalizować wpływ wartości odstających na analizę i modelowanie struktury przestrzennej, w geostatystyce stosuje się kilka metod:

  1. Transformacja danych: np. transformacja logarytmiczna, która zmniejsza skalę zmienności i osłabia wpływ wysokich wartości.

  2. Interaktywne czyszczenie semiwariogramu: Wykorzystanie „chmury semiwariogramu” (ang. variogram cloud) lub wykresów rozrzutu z przesunięciem (h-scatterplots) do identyfikacji i czasowego maskowania (wyłączania z obliczeń) par punktów generujących anomalne różnice.

  3. Wyłączenie z analizy: Jeśli wartości ekstremalne tworzą wyraźne skupienie w przestrzeni, można je uznać za osobną populację i analizować oddzielnie.

4 Przykład: Tymaczasowe wyłączenie z analizy

Jeśli wartości odstające nie są określone jako błędy pomiarowe, należy je wykluczyć wyłącznie na etapie tworzenia semiwariogramu oraz modelowania struktury przestrzennej. Wartości odstające będą brane pod uwagę podczas wykonywania estymacji.

Dane

# Wczytanie danych z pliku CSV
df = read.csv("data/cw2_dane.csv")

# Konwersja ramki danych na obiekt przestrzenny (Simple Features)
punkty = st_as_sf(df, coords = c("Xutm", "Yutm"), crs = "EPSG:32633")
  • Sprawdzenie współrzędnych
## Sprawdzenie, czy w danych nie występują punkty o tych samych współrzędnych 
# Znajdź punkty leżące w odległości mniejszej niż 0.01 m)
bliskie_punkty = st_is_within_distance(punkty, dist = 0.01)

# Wyświetl te, które mają sąsiada (długość listy > 1, bo punkt zawsze widzi samego siebie)
ktore_blisko = which(sapply(bliskie_punkty, length) > 1)
punkty[ktore_blisko, ]
Simple feature collection with 2 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 633341.4 ymin: 5809455 xmax: 633341.4 ymax: 5809455
Projected CRS: WGS 84 / UTM zone 33N
     Lp   Ba Ba_pref                 geometry
434 434 1180    1180 POINT (633341.4 5809455)
973 973  700     700 POINT (633341.4 5809455)

Punkt 434 oraz 973 mają takie same współrzędne. Algorytmy estymacji wymagają, by punkty nie miały takich samych wspołrzędnych. Istnieje kilka metod poradzenia sobie z takimi samymi współrzędnymi: uśrednienie wartości obu punktów (jeśli mają podobne wartości), usunięcie jednego z punktów, rozsunięcie punktów poprzez dodanie niewielkiego lowowego szumu, np. 1mm do współrzędnych punktów.

W przypadku puntków 434 oraz 970 uśrednimy ich wartości:

# Usuń punkty o identycznych współrzędnych, zostawiając np. tylko pierwszy z nich
#punkty = punkty %>%  distinct(geometry, .keep_all = TRUE)

punkty = punkty %>%
  group_by(geometry) %>%
  summarise(Ba = mean(Ba, na.rm = TRUE),
            Lp = first(Lp)) # pobiera numer Lp pierwszego z grupy punktów

4.1 Identyfikacja wartości globalnie odstających

# Wizualizacja rozkładu zmiennej Ba - szukanie asymetrii i wartości ekstremalnych
ggplot(df, aes(x = Ba)) + geom_histogram() + theme_bw()

W analizowanym zbiorze danych niektóre próby zostały określone jako wartości anomalne. Dla zmiennej Ba są to próby nr (Lp) 434, 557, 690, 731, 973,974,975,976,977. Poniżej utworzymy drugi zbiór danych bez wartości anaomalnych.

punkty_bez_anomalii <- punkty[!punkty$Lp%in%c(434, 557, 690, 731, 973,974,975,976,977),]
ggplot(punkty_bez_anomalii, aes(x = Ba)) + geom_histogram() + 
  labs(title = "Dane bez anomalii") + 
  theme_bw()

4.2 Wpływ wartości odstających na semiwariogram empiryczny

  • Semiwariogram empiryczny uwzględniający wszystkie punkty.
vario = variogram(Ba ~ 1, locations = punkty, width = 400, cutoff = 12000)
plot(vario)

  • Semiwariogram empiryczny na podstawie zbioru danych po wykluczeniu wartości anomalnych.
vario_bez_anomalii = variogram(Ba ~ 1, locations = punkty_bez_anomalii, width = 400, cutoff = 12000)
plot(vario_bez_anomalii)

  • Powyższy semiwariogram wykonany na podstawie zbioru danych niezawierającego wartości anomalnych zostanie użyty do modelowania struktury przestrzennej.
model = vgm(model = "Sph", psill = 2000, range = 5000, nugget = 30)
fitted_1 = fit.variogram(vario_bez_anomalii, model)
fitted_1
  model     psill range
1   Nug  815.3383     0
2   Sph 1194.6597 25169
plot(vario_bez_anomalii, fitted_1)

  • Estymację wykonuje się natomiast wykorzystując pełen zbiór danych. W parametrze locations należy wskazać pełny zbiór danych (punkty). W ten sposób algorytm do samej średniej ważonej weźmie wartości wszystkich punktów, w tym tych ekstremalnych.
siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150, res = 100, crs = "EPSG:32633")
sk_param =  gstat(formula = Ba ~ 1, 
            locations = punkty,
            model = fitted_1,
            beta = mean(punkty$Ba))

#interpolacja         
sk = interpolate(siatka, sk_param, fun = interpolate_gstat)
[using simple kriging]
[using simple kriging]
ok_param =  gstat(formula = Ba ~ 1, 
            locations = punkty,
            model = fitted_1,
            nmax = 30)

#interpolacja         
ok = interpolate(siatka, ok_param, fun = interpolate_gstat)
[using ordinary kriging]
[using ordinary kriging]
plot(sk)

plot(ok)