Ćwiczenie 1: Geostatystyczna analiza danych

Author

Anna Dmowska

library(sf)
library(gstat)
library(terra)
library(tmap)
library(tidyverse) #wczytuje także pakiet ggplot2 oraz dplyr

1 Wprowadzenie

Podstawowy przebieg geostatystycznej analiza danych obejmuje 5 etapów:

Etap 1: Eksploracyjna analiza danych.
Etap 2: Analiza i interpretacja struktury przestrzennej.
Etap 3: Matematyczne modelowanie struktury przestrzennej.
Etap 4: Estymacja.
Etap 5: Walidacja i ocena jakości estymacji.

Więcej informacji: Geostatystyka w R (rozdział 1 - wprowadzenie)

Etap 1: Eksploracyjna analiza danych

Eksploracyjna analiza danych (EDA) polega na wstępnej analizie statystycznej i wizualizacji zgromadzonych danych pomiarowych. Na tym etapie bada się rozkład statystyczny cechy (np. poprzez histogramy) oraz identyfikuje i analizuje wpływ wartości ekstremalnych i odstających, które mogą silnie zniekształcać dalsze wyniki.

Wstępna analiza statystyczna obejmuje obliczenie statystyk podstawowych, wykonanie histogramu pozwalającego na określenie rozkładu danych (czy jest to rozkład symetryczny, czy asymetryczny), określenie typu i siły korelacji zmiennych ilościowych oraz określenie istotności różnic zmiennej ilościowej w grupach wyróżnionych względem zmiennej jakościowej.
Wizualizacja danych pomiarowych w postaci mapy punktowej (nazywanej też mapą lokalizacyjną) ma na celu przede wszystkim sprawdzenie poprawności współrzędnych, określenie typu próbkowania oraz sprawdzenie poprawności danych (np. czy istnieją wartości lokalnie odstające).

Identyfikacja wartości odstających globalnie i lokalnie

Przed wykonaniem analizy geostatystycznej należy zidentyfikować, czy w zbiorze danych nie występują globalne lub lokalne wartości odstające. Nawet pojedyńcze wartości odstające i ekstremalne mają bardzo silny i zazwyczaj negatywny wpływ na wyniki analizy geostatystycznej.

Dane globalnie odstające to bardzo niskie lub bardzo wysokie wartości względem pozostałych w zbiorze danych. Można je zidentyfikować, np. poprzez obliczenie minimalnej i maksymalnej wartości zmiennej lub wykonanie histogramu oraz obejrzenie rozkładu danych. Wartości globalnie odstające na histogramie będą widoczne jako odosobnione, pojedyncze słupki znajdujące się daleko od głównej grupy danych.

Dane lokalnie odstające oznaczają nietypowe przestrzennie wartości danej cechy, np. niska wartość otoczona wysokimi wartościami lub też wysoka wartość otoczona niskimi wartościami. Do identyfikacji wartości lokalnie odstających wykorzystuje się mapy punktowe. Warto pamiętać, że nie zawsze możliwe jest określenie wartości lokalnie odstających używając statystyk opisowych czy histogramu. Wartość taka może mieścić się w przedziale danych, jednak być wyraźnie różna od otoczenia.

Przed przystąpieniem do dalszej analizy geostatystycznej należy:

zidentyfikować czy istnieją wartości globalnie lub lokalnie odstające,
zastanowić się nad przyczyną takich wartości (czy są one wynikiem błędu pomiarów?, czy są związane ze zjawiskiem, które analizujemy?),
jeśli istnieją ku temu podstawy, wartości odstające można usunąć z analizy lub czasowo zamaskować np. do obliczania semiwariogramu empirycznego.

Dodatkowe informacje

Więcej informacji nt. Eksploracyjnej analizy danych znajduje się w rozdziale 3 i 4 skryptu Geostatystyka w R opracowanego przez J.Nowosada.

Nieprzestrzenna eksploracyjna analiza danych (rozdział 3)
Przestrzenna eksploracyjna analiza danych (rozdział 4)

Etap 2: Analiza i interpretacja struktury przestrzennej

Jest to etap polegający na wykrywaniu prawidłowości w zmienności badanej cechy w przestrzeni. Wykorzystuje się w tym celu semiwariogramy empiryczne. Semiwariogram empiryczny to wykres utworzony na podstawie danych puntkowych, który pokazuje jak niepodobieństwo wartości cechy (wyrażane za pomocą semiwariancji) zmienia się wraz z odległością.

Na tym etapie można także wykonać mapę semiwariogramu pozwalająca na stwierdzenie, czy istnieją jakieś kierunki w których cecha zmienia się bardziej dynamicznie.

Etap 3: Matematyczne modelowanie struktury przestrzennej

Semiwariogram empiryczny (wyliczony z danych punktowych) jest nieciągły (wartości semiwariancji są średnimi przedziałowymi) oraz chaotyczny. Matematyczne modelowanie struktury przestrzennej polega na zastąpieniu dyskretnych (nieciągłych) miar empirycznych takich jak semiwariogram, ciągłą funkcją matematyczną (np. model sferyczny, wykładniczy lub gaussowski). Model ten stanowi podstawę dla dalszych etapów estymacji i symulacji.

Etap 4: Estymacja

Na tym etapie dokonuje się właściwego szacowania wartości cechy w miejscach, w których nie wykonano pomiarów. Do estymacji wykorzystuje się grupę algorytmów krigingu.

Etap 5: Walidacja i ocena jakości estymacji

Jest to niezbędny krok kontrolny służący do oceny wiarygodności stworzonego modelu. Polega on na porównaniu wartości oszacowanych z rzeczywistymi pomiarami, najczęściej przy użyciu procedur statystycznych takich jak kroswalidacja (cross-validation) lub walidacja podzbiorem (jackknifing).

2 Geostatystyczna analiza danych z wykorzystaniem pakietu gstat

2.1 Pakiet gstat

Pakiet gstat dostarcza funkcji do geostatystycznej analizy danych, uwzględniając m.in:

tworzenie semiwariogramów empirycznych (funkcja variogram()),
matematyczne modelowanie semiwariogramów (funkcje vgm() oraz fit.variogram())
przeprowadzenie estymacji metodą kriging (funkcja gstat(), interpolate())

Poniżej przedstawiono kod do wykonania podstawowej geostatystycznej analizy danych uwzględniającej analizę, modelowanie struktury przestrzennej oraz estymację w siatce.

# Wczytanie niezbędnych bibliotek
library(terra) # Do pracy na danych rastrowych
library(gstat) # Do analiz geostatystycznych (kriging, semiwariogramy)
library(sf)    # Do pracy na danych wektorowych (wymagane przez funkcje wewnątrz gstat)
library(ggplot2)
library(tmap)

# --- FUNKCJA POMOCNICZA ---
# Funkcja interpolate_gstat() integruje pakiet gstat z obiektami klasy SpatRaster z pakietu terra.
# Pozwala na wykonanie interpolacji bezpośrednio na siatce rastrowej. 
# Argument x to siatka interpolacyjna, argument model to parametry estymacji zdefiniowane w funkcji gstat.
interpolate_gstat = function(model, x, ...) {
  # Konwersja punktów siatki na obiekt sf z zachowaniem układu współrzędnych danych wejściowych
  v = st_as_sf(x, coords = c("x", "y"), crs = st_crs(model$data[[1]]$data))
  # Wykonanie predykcji (krigingu) dla punktów siatki
  p = predict(model, v, ...)
  # Usunięcie geometrii, aby zwrócić czyste dane tabelaryczne do obiektu rastrowego
  st_drop_geometry(p)
}

# --- PRZYGOTOWANIE DANYCH ---
# Wczytanie danych z pliku CSV
df = read.csv("data/cw_dane.csv")

# Konwersja tabeli na obiekt przestrzenny (sf) przy użyciu współrzędnych UTM
# EPSG:32633 to układ WGS 84 / UTM zone 33N
punkty = st_as_sf(df, coords = c("Xutm", "Yutm"), crs = "EPSG:32633")

# Usunięcie rekordów, które mają braki danych (NA) w kolumnie Ba
punkty = punkty[!is.na(punkty$Ba),]

# Mapa lokalizacyjna 
tm_shape(punkty) + 
  tm_symbols(fill = 'Ba', size = 0.4, fill.scale = tm_scale(values = "-brewer.spectral", n = 7)) 

# Identyfikacja wartości odstających 
ggplot(punkty, aes(x = Ba)) + geom_histogram() + theme_bw()

# Wykluczenie anomalii ze zbioru danych 
punkty_bez_anomalii = punkty |> filter(!Lp%in%c(434, 557, 690, 731, 973,974,975,976,977))

ggplot(punkty_bez_anomalii, aes(x = Ba)) + geom_histogram() + theme_bw()

# --- ANALIZA I MODELOWANIE STRUKTURY PRZESTRZENNEJ NA PODSTAWIE SEMIWARIOGRAMU ---
# Obliczenie semiwariogramu empirycznego dla zmiennej Ba
# cutoff: maksymalna odległość, width: szerokość odstępu
vario = variogram(Ba ~ 1, locations = punkty_bez_anomalii, cutoff = 12000, width = 400)
plot(vario, main = "Semiwariogram empiryczny zmiennej Ba")

# Definicja teoretycznego modelu semiwariogramu (model sferyczny "Sph")
# psill: wariancja progowa, range: zasięg, nugget: efekt samorodkowy
model = vgm(model = "Sph", psill = 1000, range = 1900, nugget = 40)

# Dopasowanie modelu teoretycznego do semiwariogramu empirycznego z wykorzystaniem funkcji fit.variogram()
fitted = fit.variogram(vario, model)
print(fitted) # Wyświetlenie parametrów dopasowanego modelu
plot(vario, fitted, main = "Dopasowany model semiwariogramu")

# --- PRZYGOTOWANIE SIATKI ---
# Utworzenie pustej siatki rastrowej (SpatRaster), dla której zostanie wykonana interpolacja
siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150,
              res = 100,  # Rozdzielczość komórki: 100m x 100m
              crs = "EPSG:32633")

# --- KRIGING PROSTY (Simple Kriging - SK) ---
# Zdefiniowanie parametrów krigingu prostego. 
# Wymaga on podania znanej średniej populacji(argument beta).
sk_param = gstat(formula = Ba ~ 1, 
            locations = punkty_bez_anomalii,
            model = fitted,
            beta = mean(punkty_bez_anomalii$Ba, na.rm = TRUE))   

# Wykonanie interpolacji metodą SK na zdefiniowanej siatce
sk = interpolate(siatka, sk_param, fun = interpolate_gstat)

# Wyświetlenie mapy wynikowej krigingu prostego z wykorzystaniem podstawowej funkcji plot()
plot(sk, main = "Kriging Prosty (SK)")

# --- KRIGING ZWYKŁY (Ordinary Kriging - OK) ---
# Zdefiniowanie parametrów krigingu zwykłego.
# nmax = 30 ogranicza obliczenia do 30 najbliższych punktów sąsiedztwa dla każdej komórki.
ok_param = gstat(formula = Ba ~ 1, 
            locations = punkty_bez_anomalii,
            model = fitted,
            nmax = 30)

# Wykonanie interpolacji metodą OK
ok = interpolate(siatka, ok_param, fun = interpolate_gstat)

# Wyświetlenie mapy wynikowej krigingu zwykłego. W wyniku otrzymujemy 2 mapy: Estymację oraz wariancję krigingową. 
plot(ok, main = "Kriging Zwykły (OK)")

# Opcjonalnie: Wykorzystanie pakietu tmap do wizualizacji danych 
library(tmap)
m1 = tm_shape(ok) +
  tm_raster(col = c("var1.pred"), col.scale = tm_scale_continuous(values = "-spectral")) + tm_title("Estymacja - kriging OK")

m2 = tm_shape(ok) +
  tm_raster(col = c("var1.var"), col.scale = tm_scale_continuous(values = "-spectral")) + tm_title("Wariancja krigingowa - kriging OK")

tmap_arrange(m1, m2)

3 Przykład: Analiza struktury przestrzennej zmiennej Ba

3.1 Etap 1: Eksploracyjna analiza danych

Dane wejściowe do analizy przechowywane są w pliku tekstowym (cw1_dane.csv), w którym znajduje się 4 zmienne:

Xutm oraz Yutm to współrzędne punktów wyrażone w układzie WGS 84/UTM zone 33N (EPSG:32633)
Lp oznaznacza numer porządkowy próbki.
Ba - zawartość Baru w próbce gleb.

Dane trzeba wczytać do obiektu ramki danych, a następnie przekształcić na obiekt przestrzenny.

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

## Określenie, w których kolumnach znajdują się współrzędne 
str(df)

'data.frame':   979 obs. of  19 variables:
 $ Xutm     : num  616894 617723 618930 619547 620602 ...
 $ Yutm     : num  5821091 5821319 5821238 5821335 5821420 ...
 $ Lp       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Sample_Id: int  141244123 141244124 141244125 141244126 141244127 141244128 141244139 141244140 141244141 141244142 ...
 $ Building : int  1 2 1 1 1 1 1 2 2 5 ...
 $ Land_use : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Soil     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Ba       : int  43 24 32 35 23 18 21 44 39 25 ...
 $ Sr       : int  10 5 25 16 4 4 4 37 12 17 ...
 $ Fe       : num  0.54 0.54 0.56 0.45 0.44 0.34 0.43 0.85 0.77 0.79 ...
 $ Ni       : int  4 4 5 3 3 3 3 6 6 7 ...
 $ Pb       : int  11 9 11 10 8 7 9 10 11 6 ...
 $ V        : int  9 9 12 8 8 7 8 13 15 12 ...
 $ PC1      : num  0.0734 0.0855 1.0629 0.0823 -0.4375 ...
 $ PC2      : num  0.5081 0.1148 1.2812 0.3883 -0.0968 ...
 $ CLS      : int  4 4 4 4 2 2 4 3 3 3 ...
 $ Ba_los   : int  43 NA 32 NA NA 18 NA 44 39 NA ...
 $ Ni_los   : int  4 NA 5 NA NA 3 NA 6 6 NA ...
 $ Ba_pref  : int  NA 24 NA 35 23 18 NA 44 39 NA ...

# 2. Konwersja ramki danych na obiekt przestrzenny z wykorzystaniem funkcji st_as_sf z pakietu sf. 
# Przypisujemy układ EPSG:32633 (UTM zone 33N), gdzie jednostką są metry
punkty = st_as_sf(df, coords = c("Xutm", "Yutm"), crs = "EPSG:32633")

## Opcjonalnie: Sprawdzenie wyniku
head(punkty)

Simple feature collection with 6 features and 17 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 616894 ymin: 5821091 xmax: 621578 ymax: 5821420
Projected CRS: WGS 84 / UTM zone 33N
  Lp Sample_Id Building Land_use Soil Ba Sr   Fe Ni Pb  V         PC1
1  1 141244123        1        1    1 43 10 0.54  4 11  9  0.07340467
2  2 141244124        2        1    1 24  5 0.54  4  9  9  0.08551053
3  3 141244125        1        1    1 32 25 0.56  5 11 12  1.06292200
4  4 141244126        1        1    1 35 16 0.45  3 10  8  0.08228017
5  5 141244127        1        1    1 23  4 0.44  3  8  8 -0.43752300
6  6 141244128        1        1    1 18  4 0.34  3  7  7 -1.37647333
          PC2 CLS Ba_los Ni_los Ba_pref                 geometry
1  0.50810033   4     43      4      NA   POINT (616894 5821091)
2  0.11476630   4     NA     NA      24 POINT (617722.6 5821319)
3  1.28115750   4     32      5      NA POINT (618930.5 5821238)
4  0.38831942   4     NA     NA      35 POINT (619546.5 5821335)
5 -0.09679255   2     NA     NA      23 POINT (620601.7 5821420)
6 -0.74025100   2     18      3      18   POINT (621578 5821284)

Obliczenie statystyk podstawowych dla analizowanej zmiennej

summary(punkty$Ba)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4.00   19.00   27.00   41.27   39.00 1180.00

Identyfikacja wartości odstających globalnie

ggplot(punkty, aes(x = Ba)) + geom_histogram() + theme_bw()

Histogram pokazuje kilka wartości odstających. Autorzy kartowania geochemicznego niektóre pomiary oznaczyli jako anomalne. Dla zmiennej Ba są to punkty o numerze porządkowym (Lp) 434, 557, 690, 731, 973,974,975,976,977. Lokalizacje te zostaną wykluczone z obliczeń.

punkty_bez_anomalii = punkty |> filter(!Lp%in%c(434, 557, 690, 731, 973,974,975,976,977))

ggplot(punkty_bez_anomalii, aes(x = Ba)) + geom_histogram() + theme_bw()

Mapa lokalizacyjna: Identyfikacja wartości odstających lokalnie

granice = read_sf("data/poznan_okolice.gpkg", layer = "granice_administracyjne")

tm_shape(granice)  + 
  tm_borders(col = "black", lwd = 2) + 
tm_shape(punkty) + 
  tm_symbols(fill = 'Ba', size = 0.4, fill.scale = tm_scale(values = "-brewer.spectral", n = 7))

granice_poznan = granice |> filter(TERYT == '306401_1')

tm_shape(punkty) + 
  tm_symbols(fill = 'Ba', size = 0.4, fill.scale = tm_scale(values = "-brewer.spectral", n = 7)) + 
tm_shape(granice)  + 
  tm_borders(col = "black", lwd = 2)

3.2 Etap 2: Analiza i interpretacja struktury przestrzennej

Stworzenie podstawowego semiwariogramu w pakiecie gstat odbywa się z użyciem funkcji variogram(). Należy w niej zdefiniować analizowaną zmienną (w tym przykładzie Ba ~ 1) oraz zbiór punktowy (punkty). Do wyświetlenia semiwariogramu empirycznego 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.

vario = variogram(Ba ~ 1, locations = punkty)
head(vario)

     np      dist    gamma dir.hor dir.ver   id
1   907  862.0113 6146.612       0       0 var1
2  3934 1469.2620 4943.536       0       0 var1
3  6966 2451.5984 4426.256       0       0 var1
4  8988 3434.5079 5194.309       0       0 var1
5 11164 4423.8467 5167.041       0       0 var1
6 13629 5425.8471 5367.089       0       0 var1

plot(vario)

plot(vario, plot.numbers = TRUE)

Szerokość odstępu

Domyślny semiwariogram empiryczny można zmodyfikować za pomocą dwóch parametrów:

cutoff - maksymalny zasięg semiwariogramu domyślnie wyliczany jako 1/3 najdłuższej przekątnej obszaru.
width - szerokość odstępu wyliczana domyślnie jako maksymalny zasięg semiwariogramu (cutoff) podzielony przez 15.

vario = variogram(Ba ~ 1, locations = punkty, cutoff = 12000, width = 800)
plot(vario)

vario = variogram(Ba ~ 1, locations = punkty, cutoff = 12000, width = 400)
plot(vario)

Proszę zwrócić uwagę na wpływ wartości odstających na semiwariogram empiryczny.

Wykluczenie wartości odstających

Poniższy semiwariogram empiryczny zostanie wykonany w oparciu o dane bez wartości anomalnych (zbiór punkty_bez_anomalii).

vario = variogram(Ba ~ 1, locations = punkty_bez_anomalii, cutoff = 12000, width = 400)
plot(vario)

Dodatkowo w funkcji plot() można zdefiniować parametry formatujące osie:

plot(vario, xlab = list(label="odstęp - h", cex=1.5), ylab = list(label="semiwariancja - gamma(h)", cex=1.5), scales = list(cex = 1.5))

head(vario)

    np      dist     gamma dir.hor dir.ver   id
1    3  273.3941  343.5000       0       0 var1
2  151  714.7730  606.7450       0       0 var1
3 1647 1009.5391  915.0000       0       0 var1
4 1687 1397.2074  952.5436       0       0 var1
5 1333 1857.3059 1039.7753       0       0 var1
6 3624 2195.7971  981.9364       0       0 var1

Obiekt vario dostarcza następujących informacji:

np: Liczba par punktów, które znalazły się w danym przedziale odległości. Jeśli w którejś klasie (wierszu) liczba par np jest bardzo mała (np. poniżej 30), punkt ten może być przypadkowy i nie należy “się do niego przywiązywać przy dopasowywaniu modelu”.
dist: Średnia odległość między punktami dla danego przedziału. Te wartości wyświetlane są na osi X wykresu semiwariogramu.
gamma: wartość semiwariancji uśredniona dla danego przedziału odległości. Te wartości wyświetlane są na osi Y wykresu semiwariogramu. Im wyższa wartość gamma, tym mniejsze podobieństwo między punktami.
dir.hor (horizontal direction), dir.ver (vertical direction): definiowany przy analizie anizotropowej.
id: identyfikator zmiennej.

Zmień parametr width na 200 oraz 2000 i zobacz jak zmienia się semiwariogram empiryczny.

3.2.1 Parametry semiwariogramu

Istnieje kilka zasad, które pozwalają na ustalenie parametrów semiwariogramu:

W każdym odstępie powinno się znaleźć co najmniej 30 par punktów.
Maksymalny zasięg semiwariogramu to 1/2 pierwiastka z badanej powierzchni (inne źródła mówią o połowie (lub 1/3) z przekątnej badanego obszaru).
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.

Więcej informacji nt. tworzenia i modyfikacji semiwariogramu znajduje się w rozdziale 6.3 skryptu Geostatystyka w R

3.3 Etap 3: Matematyczne modelowanie struktury przestrzennej

Matematyczne modelowanie struktury przestrzennej polega na dopasowaniu funkcji matematycznej do semiwariogramu empirycznego. Najczęściej stosuje się model (funkcję) sferyczny (Sph), wykładniczy (Exp), gausowski (Gau) uzupełnione o model nuggetowy (Nug).

Model semiwariogramu wymaga określenia trzech podstawowych elementów:

Nugget (efekt samorodka) - szacowana wariancja dla odstępu \(h=0\).
- Zazwyczaj semiwariogram nie zaczyna się w punkcie (0,0). Wartość nuggetu odczytujemy jako przesunięcie wartości na osi y.
Sill (semiwariancja progowa) - to maksymalny poziom zmienności danych w kontekście przestrzennym. Najczęściej, jeśli analizowana cecha jest stacjonarna, jest to wartość zbliżona do wariancji całej próby.
- Wartość semiwariancji progowej odczytujemy z osi y: wartość semiwaiancji po ustabilizowaniu się jej wartości (wartość ta już nie rośnie wraz z odległością).
Range (zasięg) - to odległość na której semiwariogram osiąga wartość semiwariancji progowej, określająca maksymalny zasięg autokorelacji danych.
- Jest to wartość odczytana z osi x (odległość) do której następuje wzrost wartości semiwariancji wraz z odległością.

W R podstawowe modele można wyświetlić używając funkcji show.vgms z pakietu gstat.

show.vgms(models = c("Sph", "Gau", "Pow", "Exp", "Nug" ), 
          range = 1.4, max = 2.5, as.groups = TRUE)

3.3.1 Dopasowanie modelu w R

Dopasowanie modelu w R wymaga:

utworzenia semiwariogramu empirycznego używając funkcji variogram().
zdefiniowania wejściowych parametrów modelu (funkcja vgm()).
dopasowanie modelu do semiwariogramu empirycznego (funkcja fit.variogram())

3.3.1.1 Zdefiniowanie parametrów modelu

W R do określenia typu modelu oraz jego podstawowych parametrów używa się funkcji vgm() z pakietu gstat. Funkcja vgm() pozwala na podanie następujących argumentów: typ modelu (model), nugget (nugget), semiwariancja progowa (psill), zasięg (range). Argumentem obowiązkowym jest typ modelu.

plot(vario)

model = vgm(model = "Sph", psill = 1000, range = 2000, nugget = 40)
model

  model psill range
1   Nug    40     0
2   Sph  1000  2000

Funkcja plot() pozwala na wizualizację dopasowania modelu do danych. Najważniejsze jest dobre dopasowanie modelu do semiwariogramu empirycznego (punkty na wykresie poniżej) w początkowym odcinku semiwariogramu (tam gdzie obserwujemy wzrost semiwariancji wraz z odleglością).

plot(vario, model)

Funkcja fit.variogram() pozwala na automatyczne dopasowanie modelu w oparciu o wstępnie podane parametry.

fitted = fit.variogram(vario, model)
fitted

  model     psill    range
1   Nug  815.3383     0.00
2   Sph 1194.7419 25170.74

plot(vario, model = fitted)

3.4 Etap 4: Estymacja

Do wykonania interpolacji jedną z metod krigingu zostanie wykorzystana funkcja gstat() z pakietu gstat. Funkcja gstat() tworzy obiekt klasy gstat, który zawiera wszystkie informacje potrzebne do wykonania krigingu. Funkcja ta wymaga zdefiniowania następujących argumentów:

formula - zmienna dla której będzie wykonywana estymacja
locations - dane wejściowe zawierające zmierzone wartości zmiennej
model - parametry modelu dopasowanego do semiwariogramu empirycznego (wynik działania funkcji vgm() lub fit.variogram())

Dodatkowo w zależności od wybranej metody należy zdefniować argument:

beta - średnia wartość stosowana przy krigingu prostym (Simple Kriging, SK)
nmax - liczba najbliższych punktów uwzględnianych w krigingu zwykłym (Ordinary kriging, OK)
maxdist - maksymalna odległość, z której mają być uwzględniane punkty w krigingu zwykłym (Ordinary kriging, OK).

W następnym kroku należy wykonać interpolację w siatce wykorzystując funkcję interpolate() z pakietu terra. Funkcja interpolate() wymaga podania 3 argumentów:

obiekt siatki
obiekt klasy gstat z parametrami krigingu
funkcję pozwalająca na wykonanie interpolacji w oparciu o obiekt klasy gstat (funkcja interpolate_gstat())

interpolate_gstat = function(model, x, ...) {
  v = st_as_sf(x, coords = c("x", "y"), crs = st_crs(model$data[[1]]$data))
  p = predict(model, v, ...)
  st_drop_geometry(p)
}

3.4.1 Kriging prosty (SK, ang. Simple Kriging)

Stworzenie siatki wykorzystując funkcję rast() z pakietu terra (utworzony zostanie obiekt klasy SpatRaster);

library(terra)
siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150,
              res = 100,  
              crs = "EPSG:32633")
siatka

class       : SpatRaster 
size        : 363, 273, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 616550, 643850, 5785850, 5822150  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633)

Stworzenie obiektu klasy gstat wykorzystując funkcję gstat() z parametrami krigingu prostego (SK). Argument beta oznacza średnią wartość analizowanej zmiennej.

sk_param = gstat(formula = Ba ~ 1, 
            locations = punkty_bez_anomalii,
            model = fitted,
            beta = mean(punkty_bez_anomalii$Ba, na.rm = TRUE))

sk_param

data:
var1 : formula = Ba`~`1 ; data dim = 970 x 17 beta = 36.38454
variograms:
        model     psill    range
var1[1]   Nug  815.3383     0.00
var1[2]   Sph 1194.7419 25170.74

wykonanie interpolacji w siatce wykorzystując funkcję interpolate() z pakietu terra.
- w funkcji interpolate() trzeba podać 3 argumenty: siatkę, obiekt klasy gstat z parametrami krigingu oraz funkcję pozwalająca na wykonanie interpolacji w oparciu obiekt gstat (funkcja interpolate_gstat()). Należy pamiętać, aby każdorazowo wykonując kriging wczytać poniższą funkcję.

interpolate_gstat = function(model, x, ...) {
  v = st_as_sf(x, coords = c("x", "y"), crs = st_crs(model$data[[1]]$data))
  p = predict(model, v, ...)
  st_drop_geometry(p)
}

sk = interpolate(siatka, sk_param, fun = interpolate_gstat)

[using simple kriging]
[using simple kriging]

Wynikowy obiekt sk zawiera 2 zmienne:
- var1.pred - to wyniki estymacji - szacowane wartości zmiennej Ba dla zdefiniowanej siatki węzłów. Wynik obliczeń przedstawia się w postaci mapy rastrowej.
- var1.var - to wariancja krigingowa.
  - Wariancja krigingowa (wariancja estymacji) służy do względnej oceny błędu (niepewności) estymacji, wynikającej z niepewności związanej z rozmieszczeniem punktów pomiarowych oraz z przyjętym modelem, a nie z konkretnych wartości badanej cechy. Wariancja krigingowa przyjmuje zazwyczaj najniższe wartości w pobliżu punktów pomiarowych, a najwyższe w obszarach rzadko opróbowanych (Im więcej punktów pomiarowych w sąsiedztwie, tym mniejszy błąd, wariancja kriginowa w punktach pomiarowych wynosi 0)

sk

class       : SpatRaster 
size        : 363, 273, 2  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 616550, 643850, 5785850, 5822150  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
source(s)   : memory
names       : var1.pred,  var1.var 
min values  :  13.24224,  907.2156 
max values  : 161.17609, 1096.0021

plot(sk)

3.4.2 Kriging zwykły (OK, ordinary kriging)

Stworzenie siatki wykorzystując funkcję rast() z pakietu terra (utworzony zostanie obiekt klasy SpatRaster);

library(terra)
siatka = rast(xmin=616550, xmax=643850, ymin=5785850, ymax = 5822150,
              res = 100,  
              crs = "EPSG:32633")
siatka

class       : SpatRaster 
size        : 363, 273, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 616550, 643850, 5785850, 5822150  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633)

Stworzenie obiektu klasy gstat wykorzystując funkcję gstat() z parametrami krigingu prostego (SK). Argument nmax, że do obliczeń wykorzystywane jest 30 najbliższych punktów.

ok_param = gstat(formula = Ba ~ 1, 
            locations = punkty_bez_anomalii,
            model = fitted,
            nmax = 30)

W krigingu zwykłym zakłada się, że średnia wartość analizowanej cechy jest stała, ale nieznana. Metoda ta uwzględnia lokalne wahania średniej. Przyjmuje się, że średnia jest stała tylko w lokalnym sąsiedztwie (tzw. ruchomym oknie) wokół punktu, którego wartość szacujemy. Parametry ruchomego okna można zdefiniować podając maksymalną liczbę najbliższych obserwacji, które mogą być użyte w obliczeniach (argument nmax) lub maksymalny możliwy promień uwzględnionego sąsiedztwa (rozmiar ruchomego okna) (argument maxdist).

ok = interpolate(siatka, sk_param, fun = interpolate_gstat)

[using simple kriging]
[using simple kriging]

4 Wykorzystanie pakietu tmap do wizualizacji wyników

Wyniki estymacji można wyświetlić także za pomocą pakietu tmap, który dostarcza więcej możliwości formatowania wizualizacji.

Wizualizacja wyników estymacji dla SK

# Rozpoczęcie tworzenia mapy na podstawie obiektu 'sk' (wynik krigingu)
sk_1 = tm_shape(sk) +
  # Wizualizacja rastra z wynikami estymacji (kolumna var1.pred)
  tm_raster(col = "var1.pred", 
            # Ustawienie skali kolorów: odwrócona paleta spectral
            col.scale = tm_scale_continuous(values = "-spectral"), 
            # Konfiguracja legendy: pozioma, pod mapą (na środku), bez ramki
            col.legend = tm_legend(orientation = "landscape", 
                                   position = tm_pos_out("center", "bottom"), 
                                   frame = FALSE, 
                                   title = "Estymacja - SK",
                                   height = 3)) + 
  # Dodanie podziałki liniowej (skali) w prawym dolnym rogu
  tm_scalebar(position = c("right", "bottom"), 
              text.size = 1.2, 
              breaks = c(0, 2,4)) +
  # Dodanie głównego tytułu mapy
  tm_title("Estymacja - SK", 
           position = tm_pos_out("center", "top"),
           size = 1.5, # Duży rozmiar czcionki tytułu
           fontface = "bold") + 
  # Ustawienia ogólnego układu mapy
  tm_layout(# Ustawienie marginesów wewnętrznych; dolny margines (0.1) robi miejsce na legendę
            inner.margins = c(0.1, 0.02, 0.05, 0.02), 
            legend.text.size = 1.2,  
            legend.title.size = 1.4,
            frame = FALSE) # Usunięcie ramki wokół całej mapy

Wizualizacja wariancji krigingowej dla SK

# Rozpoczęcie tworzenia mapy na podstawie obiektu 'sk' (wynik krigingu)
sk_2 = tm_shape(sk) +
  # Wizualizacja rastra z wynikami estymacji (kolumna var1.var)
  tm_raster(col = "var1.var", 
            # Ustawienie skali kolorów: odwrócona paleta spectral
            col.scale = tm_scale_continuous(values = "-spectral"), 
            # Konfiguracja legendy: pozioma, pod mapą (na środku), bez ramki
            col.legend = tm_legend(orientation = "landscape", 
                                   position = tm_pos_out("center", "bottom"), 
                                   frame = FALSE, 
                                   title = "Wariancja krigingowa - SK",
                                   height = 3)) + 
  # Dodanie podziałki liniowej (skali) w prawym dolnym rogu
  tm_scalebar(position = c("right", "bottom"), 
              text.size = 1.2, 
              breaks = c(0, 2,4)) +
  # Dodanie głównego tytułu mapy
  tm_title("Wariancja krigingowa - SK", 
           position = tm_pos_out("center", "top"),
           size = 1.5, # Duży rozmiar czcionki tytułu
           fontface = "bold") + 
  # Ustawienia ogólnego układu mapy
  tm_layout(# Ustawienie marginesów wewnętrznych; dolny margines (0.1) robi miejsce na legendę
            inner.margins = c(0.1, 0.02, 0.05, 0.02), 
            legend.text.size = 1.2,  
            legend.title.size = 1.4,
            frame = FALSE) # Usunięcie ramki wokół całej mapy

tmap_arrange(sk_1, sk_2)

4.1 Wykorzystanie danych przestrzennych

drogi = read_sf("data/poznan_okolice.gpkg", layer = "drogi_glowne")
granice_poznan = read_sf("data/poznan_okolice.gpkg", layer = "granice_administracyjne") |> filter(TERYT == '306401_1')

tm_shape(sk) +
  # Wizualizacja rastra z wynikami estymacji (kolumna var1.pred)
  tm_raster(col = "var1.pred", 
            # Ustawienie skali kolorów: odwrócona paleta spectral
            col.scale = tm_scale_continuous(values = "-spectral"), 
            # Konfiguracja legendy: pozioma, pod mapą (na środku), bez ramki
            col.legend = tm_legend(orientation = "landscape", 
                                   position = tm_pos_out("center", "bottom"), 
                                   frame = FALSE, 
                                   title = "Estymacja - SK",
                                   height = 3)) + 
tm_shape(granice_poznan) + tm_borders(col = "black", lwd = 1) + 
tm_shape(drogi) + tm_lines(col = "grey40", lwd = 1) + 
  # Dodanie podziałki liniowej (skali) w prawym dolnym rogu
  tm_scalebar(position = c("right", "bottom"), 
              text.size = 1.2, 
              breaks = c(0, 2,4)) +
  # Dodanie głównego tytułu mapy
  tm_title("Estymacja - SK", 
           position = tm_pos_out("center", "top"),
           size = 1.5, # Duży rozmiar czcionki tytułu
           fontface = "bold") + 
  # Ustawienia ogólnego układu mapy
  tm_layout(# Ustawienie marginesów wewnętrznych; dolny margines (0.1) robi miejsce na legendę
            inner.margins = c(0.1, 0.02, 0.05, 0.02), 
            legend.text.size = 1.2,  
            legend.title.size = 1.4,
            frame = FALSE) # Usunięcie ramki wokół całej mapy

5 Wykorzystanie pakietu ggplot2 do wizualizacji semiwariogramów

W niektórych sytuacjach potrzebne jest zwizualizowanie kilku semiwariogramów empirycznych na jednym wykresie. Można w tym celu wykorzystać pakiet ggplot2.

v1 = variogram(Ba ~ 1, locations = punkty_bez_anomalii)

head(v1)

     np      dist     gamma dir.hor dir.ver   id
1   879  875.5992  849.6559       0       0 var1
2  3843 1466.9138  970.2134       0       0 var1
3  6831 2451.2030 1008.5306       0       0 var1
4  8784 3433.1966 1077.3258       0       0 var1
5 10921 4422.9259 1073.4742       0       0 var1
6 13326 5424.7998 1178.3113       0       0 var1

Zmienna dist zawiera informacje o średniej odległości między punktami dla danego przedziału (jest to wartość wyświetlana na osi X). Zmienna gamma to wartość semiwariancji uśredniona dla danego przedziału odległości (wyświetlane są na osi Y wykresu semiwariogramu).

ggplot(v1, aes(x = dist, y = gamma)) + 
  geom_point(color = "steelblue", size = 2) + 
  expand_limits(y = 0) + 
  theme_bw()

library(ggplot2)
ggplot(v1, aes(x = dist, y = gamma)) +
  # Punkty empiryczne - powiększone dla lepszej widoczności
  geom_point(color = "steelblue", size = 2.5) + 
  # Linia łącząca punkty (pomocnicza)
  geom_line(color = "darkgrey", linetype = "dashed", alpha = 0.6) +
  geom_text(aes(label = np), 
            vjust = -1.2,      # Przesunięcie etykiety nad punkt (pionowo)
            size = 3,          # Rozmiar czcionki etykiet
            color = "grey30") + # Kolor czcionki
  labs(
    title = "Semiwariogram empiryczny dla Baru (Ba)",
    x = "Odstęp - h [m]",
    y = "Semiwariancja - gamma",
  ) +
  # Ustawienie początku osi na zero
  expand_limits(y = 0) +
  theme_bw() +
  # Dodatkowe dopracowanie wyglądu
  theme(
    plot.title = element_text(face = "bold", size = 16),
    axis.title = element_text(size = 14),
    panel.grid.minor = element_blank()
  )

5.1 Porównanie semiwariogramów empirycznych

Poniższy kod przedstawia 3 warianty wykonania semiwariogramu empirycznego dla różnych wartości odstępu.

v1 = variogram(Ba ~ 1, locations = punkty_bez_anomalii)
v2 = variogram(Ba ~ 1, locations = punkty_bez_anomalii, cutoff = 12000, width = 800)
v3 = variogram(Ba ~ 1, locations = punkty_bez_anomalii, cutoff = 12000, width = 400)

W celu zwizualizowania kilku semiwariogramów empirycznych na jednym wykresie obiekty v1, v2 oraz v3 połączymy w jedną ramkę danych.

# 1. Dodanie kolumny identyfikującej do każdego obiektu
v1$parametry = "Domyślne"
v2$parametry = "Odstęp: 800"
v3$parametry = "Odstęp: 400"
# 2. Połączenie obiektów w jedną ramkę danych (long format)
v_all = rbind(v1, v2, v3)

# 3. Przykładowy wykres porównawczy w ggplot2
ggplot(v_all, aes(x = dist, y = gamma, color = parametry)) +
  geom_point(size = 2) +
  geom_line(linewidth = 1) +
  labs(
    title = "Porównanie semiwariogramów empirycznych dla zmiennej Ba",
    x = "Odstęp - h [m]",
    y = "Semiwariancja - gamma",
    color = "Ustawienia:"
  ) +
  theme_bw()

Zwróć uwagę na różnice między semiwariogramami w zależności od wybranych parametrów.

5.2 Wizualizacja modelu w ggplot2

v3 = variogram(Ba ~ 1, locations = punkty_bez_anomalii, cutoff = 12000, width = 400)
# Dopasowanie modelu sferycznego 
fitted_sph = fit.variogram(v3, model = vgm(model = "Sph", psill = 1000, range = 2000, nugget = 40))
# Dopasowanie modelu wykładniczego
fitted_exp = fit.variogram(v3, model = vgm(model = "Exp", psill = 1000, range = 2000, nugget = 40))

Aby zwizualizować model używając pakietu ggplot2, w pierwszym kroku na podstawie parametrów modelu trzeba wygenerować linię używając funkcji variogramLine().

# maxdist powinien odpowiadać maksymalnej odległości z Twojego semiwariogramu empirycznego
## n oznacza liczbę punktów dla których zostaną obliczone wartości semiwariancji na podstawie parametrów modelu
v_model_sph_line = variogramLine(fitted_sph, maxdist = max(v3$dist), n = 400)

head(v_model_sph_line)

          dist    gamma
1   0.01181441 815.3391
2  29.62182469 817.4473
3  59.23183497 819.5555
4  88.84184526 821.6637
5 118.45185554 823.7718
6 148.06186583 825.8799

Sprawdź co robi funkcja variogramLine.

v_model_exp_line = variogramLine(fitted_exp, maxdist = max(v3$dist), n = 400)

ggplot() + 
  # 1. Punkty empiryczne
  geom_point(data = v3, aes(x = dist, y = gamma, color = "Semiwariogram empiryczny"), size = 2.5) +
  # 2. Model sferyczny
  geom_line(data = v_model_sph_line, aes(x = dist, y = gamma, color = "Model sferyczny"), linewidth = 1) +
  # 3. Model wykładniczy
  geom_line(data = v_model_exp_line, aes(x = dist, y = gamma, color = "Model wykładniczy"), linewidth = 1) +
  
  # Ręczne przypisanie kolorów do etykiet z aes()
  scale_color_manual(
    values = c( name = "",
      "Semiwariogram empiryczny" = "black", 
      "Model sferyczny" = "darkred", 
      "Model wykładniczy" = "darkgreen"
    )
  ) +
  
  # Estetyka wykresu
  labs(
    title = "Porównanie modeli teoretycznych",
    x = "Odstęp - h [m]",
    y = "Semiwariancja - gamma"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")