library("terra")

Wstęp

Interpolacja przestrzenna jest to metoda wykorzystywana do oszacowania nieznanych wartości zmiennej (cechy) na badanym obszarze na podstawie zmierzonych wartości w określonych lokalizacjach. Dzięki niej możliwe jest wypełnienie luk w zbiorze danych, tj. obszarów nieobjętych pomiarem.

Istnieją dwa zasadnicze podejścia do interpolacji przestrzennej:

  • modele deterministyczne – zakładają, że istnieje prosta zależność przestrzenna bez uwzględnienia losowości i niepewności:
    • Naturalna interpolacja sąsiadów (Natural neighbor interpolation)
    • Interpolacja wielomianowa (Polynomial interpolation)
    • Odwrotne ważenie odległości (Inverse Distance Weighting)
    • Interpolacja funkcjami sklejanymi (Spline interpolation)
    • Interpolacja radialną funkcją bazową (Radial Basis Function interpolation)
  • modele geostatystyczne – uwzględniają właściwości statystyczne rozkładu przestrzennego danych (autokorelację), dzięki czemu możliwe jest zbadanie wariancji w zależności od kierunku i odległości:
    • Prosty kriging (Simple Kriging)
    • Zwyczajny kriging (Ordinary Kriging)
    • Uniwersalny kriging (Universal Kriging)
    • Kriging blokowy (Block Kriging)
    • Co-kriging

Interpolacja

Przygotowanie danych

W katalogu dane znajdziesz plik dane_meteo.csv z pomiarami średniej temperatury dobowej (°C) oraz sumy dobowej opadów (mm) pochodzących ze stacji meteorologicznych Instytutu Meteorologii i Gospodarki Wodnej (https://danepubliczne.imgw.pl/).

Wczytajmy go używając funkcji read.csv() i wyświetlmy jego strukturę za pomocą funkcji str().

meteo_df = read.csv("../dane/dane_meteo.csv")
str(meteo_df)
## 'data.frame':    110 obs. of  6 variables:
##  $ X    : num  19 20 20 20.7 21.8 ...
##  $ Y    : num  49.8 49.3 49.2 49.6 49.7 ...
##  $ NAZWA: chr  "BIELSKO-BIAŁA" "ZAKOPANE" "KASPROWY WIERCH" "NOWY SĄCZ" ...
##  $ TEMP : num  8.8 6.9 -2.4 10.1 10.1 8.6 9 -0.5 8.9 10.1 ...
##  $ OPAD : num  3.7 0.2 0.9 1.6 0.1 5.5 2.8 1.4 2.9 7.2 ...
##  $ TYP  : chr  "synop" "synop" "synop" "synop" ...

Do obliczenia podstawowych statystyk opisowych temperatury możemy wykorzystać funkcję summary().

summary(meteo_df$TEMP)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -2.400   6.325   8.450   7.559   9.200  10.800

Dodatkowo, możemy sprawdzić rozkład wartości temperatury używając funkcji hist().

hist(meteo_df$TEMP, main = NULL, breaks = 10, xlab = "Temperatura [°C]",
     ylab = "Częstość" )

W kolejnym kroku dokonajmy konwersji ramki danych do obiektu przestrzennego, tj. SpatVector. W tym celu należy użyć funkcji vect() określając nazwy kolumn z długością i szerokością geograficzną (argument geom) oraz układ współrzędnych (argument crs).

meteo = vect(meteo_df, geom = c("X", "Y"), crs = "EPSG:4326")
meteo
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 110, 4  (geometries, attributes)
##  extent      : 14.242, 24.036, 49.214, 54.754  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :           NAZWA  TEMP  OPAD   TYP
##  type        :           <chr> <num> <num> <chr>
##  values      :   BIELSKO-BIAŁA   8.8   3.7 synop
##                       ZAKOPANE   6.9   0.2 synop
##                KASPROWY WIERCH  -2.4   0.9 synop

Następnie możemy wyświetlić dane z uwzględnieniem typu stacji meteorologicznej.

plot(meteo, "TYP", col = c("purple", "blue"), alpha = 0.7,
     main = "Stacje meteorologiczne")

Oraz sprawdzić przestrzenny rozkład temperatur definiując wcześniej wybraną paletę kolorów, np. z funkcji hcl.colors().

paleta = hcl.colors(10, palette = "RdYlBu", rev = TRUE)
plot(meteo, "TEMP", type = "continuous", col = paleta,
     main = "Temperatura [°C]")

Zasadniczym celem interpolacji jest stworzenie ciągłej powierzchni (rastra) na podstawie danych dyskretnych (punktowych). Wymaga to zdefiniowania siatki, dla której każdej komórce zostanie przypisana estymowana wartość. Do stworzenia takiej siatki można zastosować funkcję rast() określając podstawowe parametry takie jak liczba wierszy i kolumn, zakres przestrzenny oraz układ przestrzenny.

W naszym przypadku możemy skopiować metadane z obiektu wektorowego meteo i wskazać rozdzielczość jako 0.01 stopnia.

r = rast(meteo, resolution = 0.01)
r
## class       : SpatRaster 
## dimensions  : 554, 979, 1  (nrow, ncol, nlyr)
## resolution  : 0.01, 0.01  (x, y)
## extent      : 14.242, 24.032, 49.214, 54.754  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)

Walidacja wyników

Tworzenie modeli zawsze wymaga sprawdzenia ich skuteczności na niezależnym zbiorze danych. Proces ten nazywany jest walidacją (inaczej testowaniem). Celem walidacji jest ocena, jak model będzie działał na nowych, wcześniej niewidzianych danych, co ma kluczowe znaczenie dla określenia jego poprawności i przydatności. Wyróżnić można kilka metod walidacji, jednak dwie najważniejsze to:

  • walidacja podzbiorem (hold-out validation) – zbiór danych wejściowych dzielony jest na dwa podzbiory: zbiór treningowy (używany do trenowania modelu) oraz zbiór testowy (używany do oceny skuteczności wytrenowanego modelu)
  • walidacja krzyżowa (k-fold cross-validation) – zbiór danych wejściowych dzielony jest na kilka mniejszych podzbiorów (fold), np. 5 lub 10. Każdy podzbiór używany jest jednokrotnie jako testowy, a pozostałe podzbiory używane są do trenowania modelu. Proces ten powtarzany jest \(k\) razy, zapewniając, że wszystkie podzbiory wykorzystywane są zarówno do szkolenia, jak i testowania. Umożliwia to bardziej wiarygodne oszacowanie skuteczności modelu.

Na potrzeby niniejszej analizy wykorzystamy pierwszą prostszą metodę, tj. walidację podzbiorem. Ogólnie przyjmuje się, że dane testowe powinny stanowić około 30% wejściowego zbioru danych, natomiast dane treningowe około 70%. Podziału na podzbiory można dokonać używając funkcji sample() w celu wylosowania indeksów obserwacji, które trafią do podzbioru treningowego. Pozostałe obserwacje można przypisać do podzbioru testowego.

Zauważ, że losowanie za każdym razem zwraca różne indeksy. Oznacza to, że podczas ponownego wykonywania skryptu otrzymujemy inny wynik, co tym samym powoduje, że analiza przestaje być powtarzalna / odtwarzalna. Aby temu zapobiec należy ustawić ziarno losowości przy pomocy funkcji set.seed().

set.seed(1)
n = nrow(meteo_df)
indeksy = sample(n, size = 0.7 * n)
trening = meteo_df[indeksy, ]
test = meteo_df[-indeksy, ]

Sprawdźmy podział punktów pomiarowych na mapie.

plot(meteo[indeksy], col = "green", alpha = 0.8)
plot(meteo[-indeksy], col = "blue", alpha = 0.8, add = TRUE)
add_legend("bottomleft", legend = c("Treningowe", "Testowe"),
           col = c("green", "blue"), pch = 19, cex = 0.9)

Walidacja modelu obejmuje wykorzystanie różnych wskaźników (metryk) do oceny jego skuteczności. W zależności od problemu badawczego i celu oceny można zastosować różne metryki. W naszej analizie wykorzystamy najpopularniejszą metrykę stosowaną w problemach z zakresu regresji, tj. spierwiastkowany błąd średniokwadratowy (Root Mean Squared Error), niemniej nie jest to jedyny wskaźnik.

RMSE = function(obserwowane, predykcja) {
  sqrt(mean((predykcja - obserwowane) ^ 2))
}

Im większa wartość tej metryki, tym większy błąd.

Metody interpolacji

Pakiet gstat oferuje szereg metod do modelowania przestrzennego, które są dostępne za pomocą funkcji o tej samej nazwie, tj. gstat().

library("gstat")

Naturalna interpolacja sąsiadów

Naturalna interpolacja sąsiadów wymaga zdefiniowana pięciu argumentów:

  1. wzoru modelu – formula. W naszym przykładzie nie wykorzystujemy żadnych dodatkowych zmiennych (tj. zmiennych wyjaśniających takich jak wysokość terenu czy odległości od zbiorników wodnych). Modelujemy wyłącznie rozkład temperatury. Zapis wzoru wygląda w następujący sposób: formula = TEMP ~ 1 (można to rozumieć, że temperatura modelowana jest od samej siebie; uwzględniając jedynie rozkład przestrzenny).
  2. nazwy współrzędnych geograficznych – locations. W naszym przykładzie: locations = ~X + Y.
  3. zbiór danych – data.
  4. liczba najbliższych obserwacji, które zostaną wykorzystane do predykcji – nmax.
  5. wykładnik odwrotności odległości (inverse distance power). W przypadku tej metody jego wartość należy ustawić na 0, żeby wszystkie obserwacje miały równą wagę: set = list(idp = 0).
mdl = gstat(formula = TEMP ~ 1, locations = ~X + Y, data = trening, nmax = 10,
            set = list(idp = 0))

W ten sposób opracowaliśmy pierwszy model predykcyjny. Tak wytrenowany model możemy zastosować dla całego obszaru analizy używając funkcji interpolate(), w której określimy raster, model, nazwy kolumn ze współrzędnymi (xyNames).

nn = interpolate(r, mdl, xyNames = c("X", "Y"), debug.level = 0)
nn = subset(nn, 1) # wybiera tylko pierwszą warstwę
plot(nn, col = paleta)

Walidacji wyników można dokonać za pomocą funkcji predict(). Predykcję przeprowadzimy jedynie dla punktów testowych, które nie były wykorzystane do wytrenowania modelu. W przeciwieństwie do funkcji interpolate() najpierw należy zdefiniować model, a następnie zbiór danych.

nn_test = predict(mdl, test, debug.level = 0)
# wybierz kolumnę z prognozowanymi wartościami
nn_test = nn_test$var1.pred
nn_test
##  [1] 7.68 7.65 9.35 8.79 9.02 8.99 9.09 9.26 9.00 9.42 8.37 8.72 9.10 6.95 7.44
## [16] 7.96 5.04 4.22 6.87 7.69 9.00 7.68 7.68 8.54 8.79 7.75 7.75 9.09 9.26 8.19
## [31] 7.72 6.15 8.86

W ten sposób otrzymaliśmy wartości prognozowane przez model dla punktów ze zbioru testowego. W kolejnym kroku wyliczmy błąd prognozy za pomocą wcześniej zdefiniowanej funkcji RMSE().

RMSE(test$TEMP, nn_test)
## [1] 2.665025

Interpolacja wielomianowa

Podobnie jak w poprzednim przykładzie obligatoryjnie musimy zdefiniować argumenty formula, locations oraz data. Jednak, aby wykonać interpolację za pomocą funkcji wielomianowych, dodatkowo trzeba zdefiniować stopień wielomianu (argument degree) z przedziału od 1 do 3. Pozostałe kroki pozostają bez zmian.

mdl = gstat(formula = TEMP ~ 1, locations = ~X + Y, data = trening,
            degree = 3) # wielomian trzeciego stopnia
poly = interpolate(r, mdl, xyNames = c("X", "Y"), debug.level = 0)
poly = subset(poly, 1)
plot(poly, col = paleta)

poly_test = predict(mdl, test, debug.level = 0)$var1.pred
RMSE(test$TEMP, poly_test)
## [1] 2.588132

Odwrotne ważenie odległości

Metoda odwrotnej ważonej odległości zakłada, że otaczające punkty wpływają na przewidywaną wartość komórki na podstawie odwrotności ich odległości, czyli punkty, które są położone dalej, mają mniejszą wagę w predykcji wartości. Wykładnik odwrotności odległości ustawia się za pomocą argumentu set = list(idp = 1).

mdl = gstat(formula = TEMP ~ 1, locations = ~X + Y, data = trening,
            set = list(idp = 2))
idw = interpolate(r, mdl, xyNames = c("X", "Y"), debug.level = 0)
idw = subset(idw, 1)
plot(idw, col = paleta)

idw_test = predict(mdl, test, debug.level = 0)$var1.pred
RMSE(test$TEMP, idw_test)
## [1] 2.300736

Kriging

W porównaniu do zaprezentowanych metod interpolacji, modele geostatystyczne są bardziej skomplikowane, ponieważ wymagają zdefiniowania modelu na podstawie wariogramu, który służy do oceny autokorelacji w ujęciu przestrzennym. W tym celu najpierw należy stworzyć obiekt gstat, a następnie wykorzystać funkcję variogram(). Szerokość przedziałów odległości (argument width) jest dobierana automatycznie, jednak najczęściej wartość optymalna wymaga samodzielnego wyboru metodą prób i błędów.

gst = gstat(formula = TEMP ~ 1, locations = ~X + Y, data = trening)
v = variogram(gst, width = 0.4)
plot(v)

Cała sztuka modelowania geostatystycznego polega na dopasowaniu modeli matematycznych do punktów na wariogramie. Dostępne modele w pakiecie gstat można wyświetlić za pomocą funkcji show.vgms(). Do manualnego zdefiniowania modelu służy funkcja vgm(), dla której należy określić następujące parametry:

  • psill (próg) – wartość, przy której wariogram się wyrównuje, reprezentując maksymalną wariancję lub całkowitą wariancję danych. Poza odległością, dla której próg jest osiągnięty, punkty są uważane za nieskorelowane.
  • range (zasięg) – zakres, w którym punkty są skorelowane przestrzennie.
  • nugget – wariancja przy zerowej odległości związana z błędem pomiaru lub krótkozasięgową zmiennością przestrzenną występująca w odległościach mniejszych niż interwał próbkowania.
# model sferyczny
mdl_sph = vgm(psill = 3, model = "Sph", range = 2, nugget = 1)
plot(v, model = mdl_sph)

Dodatkowo, istnieje możliwość automatycznej optymalizacji parametrów modelu przy pomocy funkcji fit.variogram().

fv = fit.variogram(v, mdl_sph)
fv
##   model    psill    range
## 1   Nug 2.117132  0.00000
## 2   Sph 8.756669 21.61317
plot(v, model = fv)

Po zdefiniowaniu modelu (i optymalizacji jego parametrów) możemy przejść do predykcji używając krigingu zwyczajnego, który zakłada, że średnia modelowanej zmiennej jest stała, ale nieznana w lokalnym otoczeniu interpolacji. Procedura predykcji wygląda identycznie jak w poprzednich przypadkach, jednak warto zauważyć, że zwracane są dwa obiekty. Pierwszy reprezentuje estymowane wartości modelowanej zmiennej, natomiast drugi wariancje estymacji (wariancja rośnie wraz z odległością od punktu pomiarowego).

mdl = gstat(formula = TEMP ~ 1, locations = ~X + Y, data = trening, model = fv)
kr = interpolate(r, mdl, xyNames = c("X", "Y"), debug.level = 0)
names(kr) = c("Predykcja", "Wariancja")
par(mfrow = c(1, 2))
plot(kr[[1]], col = paleta, main = "Predykcja")
plot(kr[[2]], col = gray.colors(n = 10, rev = TRUE), main = "Wariancja")

kr_test = predict(mdl, test, debug.level = 0)$var1.pred
RMSE(test$TEMP, kr_test)
## [1] 2.50999

Jeśli chcesz rozszerzyć swoją wiedzę z zakresu statystyki przestrzennej, to koniecznie sprawdź podręcznik “Geostatystyka w R” autorstwa Jakuba Nowosada.

Zadanie

8. Porównaj wymienione metody dla dobowej sumy opadów (meteo_df$OPAD). Dodatkowo wykorzystaj metodę “cienkiej płytki” (thin plate spline) z pakietu fields (funkcja Tps()). Zwróć uwagę, że prognozowana wartość opadu nie powinna być ujemna.