Geostatystyka czasoprzestrzenna

library(spacetime)
library(gstat)
library(tidyverse)
library(sf)
library(sp)

load("dane/r5to10.RData")

data("air", package = "spacetime")

1 Semiwariancja czasoprzestrzenna

Semiwariancja czasoprzestrzenna obliczana jest dla odstępów określanych poprzez odległość oraz czas. Np. obliczamy semiwariancję dla przesunięcia przestrzennego \(h = 10\) km i czasowego \(u = 1\) h.

\[ \gamma(\mathbf{h}, u) = \frac{1}{2N(\mathbf{h}, u)} \sum_{i=1}^{N(\mathbf{h}, u)} [Z(\mathbf{s}_i + \mathbf{h}, t_i + u) - Z(\mathbf{s}_i, t_i)]^2 \] , gdzie

• \(\gamma(\mathbf{h}, u)\): wartość semiwariancji dla wektora odległości \(\mathbf{h}\) oraz opóźnienia czasowego \(u\).

• \(\mathbf{h}\): separacja przestrzenna - odstęp (odległość i kierunek między punktami).

• \(u\): separacja czasowa - opóźnienie czasowe (różnica czasu między pomiarami).

• \(N(\mathbf{h}, u)\): liczba par obserwacji oddzielonych od siebie o odstęp \(h\) i czas \(u\)

• \(Z(\mathbf{s}_i, t_i)\): wartość cechy zmierzona w punkcie \(s\) i czasie \(t\)

• \(Z(\mathbf{s}_i + \mathbf{h}, t_i + u)\): wartość cechy w punkcie oddalonym o odstęp \(\mathbf{h}\) i czasie przesuniętym o \(u\).

1.1 Obliczenie semiwariancji czasoprzestrzennej

Dane wejściowe

Dysponujemy pomiarami temperatury w dwóch stacjach (A i B) oddalonych od siebie o \(h = 10\) km. Pomiary wykonano w dwóch momentach czasu (\(t_1\) i \(t_2\)) w odstępie \(u = 1\) godziny.

  Stacja  Czas  z
1      A 10:00 15
2      B 10:00 17
3      A 11:00 18
4      B 11:00 21

Zadanie

Oblicz semiwariancję czasoprzestrzenną dla przesunięcia przestrzennego \(h = 10\) km i czasowego \(u = 1\) h.

Rozwiązanie

Krok 1: Wybór par punktów spełniających warunek \((h, u)\)

• Musimy znaleźć wszystkie pary pomiarów, które są oddalone o 10 km ORAZ o 1 godzinę.

• W naszym zbiorze są to dwie takie pary:

  • Para 1: Stacja A o 10:00 (\(z=15\)) i Stacja B o 11:00 (\(z=21\)).
  • Para 2: Stacja B o 10:00 (\(z=17\)) i Stacja A o 11:00 (\(z=18\)).

Krok 2: Obliczenie różnicy kwadratów

• Różnica kwadratów dla 1. pary: \((21 - 15)^2 = 6^2 = \mathbf{36}\)

• Różnica kwadratów dla 2. pary: \((18 - 17)^2 = 1^2 = \mathbf{1}\)

• Suma kwadratów różnic = \(36 + 1 = \mathbf{37}\)

Krok 3: Wynik końcowy

\[ \hat{\gamma}(10\text{km}, 1\text{h}) = \frac{1}{2 \times 2} \times 37 = \frac{37}{4} = \mathbf{9,25} \]

Wartość 9,25 to semiwariancja dla tej konkretnej odległości (punkty oddalone o 10km) i czasu (pomiary wykonane w odstępie 1 godziny).

2 Geostatystyka czasoprzestrzenna w R

Pakiet gstat dostarcza funkcji do geostatystycznej analizy danych czasoprzestrzennych w oparciu o obiekty klasy ST (STFDF, STSDF, STIDF) z pakietu spacetime.

Funkcje pakietu gstat:

• variogramST():Obliczanie empirycznego semiwariogramu czasoprzestrzennego

  • https://r-spatial.github.io/gstat/reference/variogramST.html

• vgmST(): Definiowanie struktury teoretycznego modelu semiwariogramu czasoprzestrzennego.

  • https://r-spatial.github.io/gstat/reference/vgmST.html

• fit.StVariogram(): Automatyczne dopasowanie parametrów modelu do danych empirycznych.

  • https://r-spatial.github.io/gstat/reference/fit.StVariogram.html

• krigeST(): Wykonywanie krigingu czasoprzestrzennego.

  • https://r-spatial.github.io/gstat/reference/krigeST.html

2.1 Główne kroki analizy czasoprzestrzennej

Etap 1: Przygotowanie danych

Dane surowe należy przekształcić w ustrukturyzowane obiekty czasoprzestrzenne klas takich jak STFDF (pełna siatka), STSDF (rzadka siatka) lub STIDF (dane nieregularne). Wykorzystuje się do tego funkcje stConstruct() z pakietu spacetime lub funkcje dedykowane dla określonej klasy, np. STFDF().

Etap 2. Analiza eksploracyjna (EDA) i wizualizacja

Do wizualizacji wykorzystuje się funkcję stplot(), która pozwala m.in na tworzenie diagramów Hovmöllera (mode = “xt) oraz paneli map dla różnych terminów. Wizualizacje pozwalają przeanalizować strukturę danych (np. występowanie braków danych), występowanie trendów, autokorelacji.

Etap 3: Obliczanie semiwariogramu empirycznego

Semiwariogram czasoprzestrzenny podsumowuje semiwariancję dla różnych kombinacji odległości w przestrzeni (\(h\)) i opóźnień w czasie (\(u\)).

Do obliczenia semiwariogramu empirycznego służy funkcja variogramST() z pakietu gstat. Jej kluczowe argumenty to: width (szerokość odstępu “przestrzennego”), cutoff (maksymalna odległość), tlags (opóźnienia czasowe).

Etap 4: Budowa modelu teoretycznego oraz optymalizacja parametrów modelu

Na tym etapie do semiwariogramu empirycznego dopasowywany jest model teoretyczny. Dostępne klasy modeli w funkcji vgmST(): separable (iloczyn kowariancji), metric (czas jako dodatkowy wymiar), productSum, sumMetric oraz simpleSumMetric.

Do dopasowania modelu teoretycznego służy funkcja vgmST(), natomiast funkcja fit.StVariogram() optymalizuje parametry modelu. Funkcja estiStAni() pomaga oszacować współczynnik anizotropii czasoprzestrzennej.

Etap 5: Kriging

Na podstawie dopasowanego modelu generowane są prognozy dla niepomierzonych punktów w przestrzeni i czasie. Do predykcji wykorzystuje się funkcje krigeST(). Kriging wymaga stworzenia obiektu z lokalizacjami do predykcji (siatka, punkty) w formacie klasy ST (lokalizacje określone w przestrzeni i czasie).

3 Semiwariogram czasoprzestrzenny

Empiryczny semiwariogram czasoprzestrzenny obliczany jest za pomocą funkcji variogramST() z pakietu gstat. W przeciwieństwie do standardowego semiwariogramu, tutaj analizujemy zmienność nie tylko w funkcji odległości (\(h\)), ale również w funkcji czasu (\(u\)).

Funkcja variogramST() wymaga zdefiniowania następujących parametrów:

• formula: Określa cechę, dla której liczona jest semiwariancja.

• data: Obiekt zawierający dane czasoprzestrzenne. Musi on należeć do jednej z klas zdefiniowanych w pakiecie spacetime, takich jak:

  • STFDF: Pełna siatka (każda stacja ma pomiar w każdym momencie).
  • STSDF: Dane rzadkie (niektóre stacje w niektórych momentach).
  • STIDF: Dane nieregularne (punkty w dowolnych miejscach i dowolnym czasie).

• cutoff: maksymalna odległość, do której pary punktów są uwzględniane w obliczeniach. Pary leżące dalej niż to zostało określone przez parametrt cutoff są ignorowane. Domyślnie jest to 1/3 przekątnej obszaru badań.

• width: Szerokość odstępu “przestrzennego” h (lag).

• tlags: Określa opóźnienia czasowe (\(u\)), dla których ma być liczona semiwariancja. Dla danych regularnych (STFDF) podaje się liczbę kroków czasowych, np. 0:5 oznacza przesunięcie o 0, 1, 2, 3, 4 i 5 jednostek czasu (np. dni).

# Obliczanie empirycznego wariogramu czasoprzestrzennego
vv = variogramST(
  formula = PM10 ~ 1,      # Definicja zmiennej: badamy PM10
  data = r5to10,           # Obiekt klasy STFDF zawierający dane 
  width = 20,              # Szerokość odstępu h przestrzennego
  cutoff = 200,            # Maksymalna odległość: analiza kończy się na odległości 200 km (pary dalsze są pomijane)
  tlags = 0:5              # Opóźnienie czasowe: oblicz zmienność dla różnicy czasu wynoszącej 0, 1, 2, 3, 4 i 5 dni
)

data(vv)

str(vv)

Classes 'StVariogram' and 'data.frame': 66 obs. of  7 variables:
 $ np      : num  0 3456 16554 21257 33547 ...
 $ dist    : num  NA 13.8 29.5 49.1 70.4 ...
 $ gamma   : num  NA 11 23.4 29.6 31.2 ...
 $ id      : chr  "lag0" "lag0" "lag0" "lag0" ...
 $ timelag : num  0 0 0 0 0 0 0 0 0 0 ...
  ..- attr(*, "units")= chr "days"
 $ spacelag: num  0 10 30 50 70 90 110 130 150 170 ...
  ..- attr(*, "units")= chr "km"
 $ avgDist : num  0 13.8 29.5 49.1 70.4 ...
 - attr(*, "boundaries")= num [1:11] 0 20 40 60 80 100 120 140 160 180 ...

head(vv)

     np     dist    gamma   id timelag spacelag  avgDist
1     0       NA       NA lag0       0        0  0.00000
2  3456 13.81028 10.98393 lag0       0       10 13.80791
3 16554 29.50171 23.41684 lag0       0       30 29.50120
4 21257 49.09772 29.57217 lag0       0       50 49.09900
5 33547 70.36967 31.23434 lag0       0       70 70.36825
6 26053 91.59265 34.77332 lag0       0       90 91.59565

3.1 Mapa semiwariancji czasoprzestrzennej

Podstawową formą prezentacji semiwariancji czasoprzestrzennej w pakiecie gstat jest mapa semiwariancji czasoprzestrzennej. Mapa semiwariancji czasoprzestrzennej na osi poziomej (X) przedstawia odległość “przestrzenną” (np. w km). Maksymalna odległość, dla której wykonywana jest mapa określana jest przez parametr cutoff w funkcji variogramST(), a rozdzielczość pojedyńczej komórki przez parametr width. Oś pionowa (Y) reprezentuje opóźnienie czasu zdefiniowane za pomocą parametru tlags (dni). Kolory reprezentują wartość semiwariancji (\(\gamma\)) obliczoną dla określonego odstępu odległości (\(h\)) i opóźnienia czasowego (\(u\)).

plot(vv)

Pierwszy rząd (time lag 0; \(u = 0\)) pokazuje “czystą” zmienność przestrzenną, określając jak semiwariancja zmienia się wraz z odległością (dla pomiarów wykonanych w tym samym czasie). Pierwsza kolumna (distance równy 0, \(h = 0\)) pokazuje “czystą” zmienność czasową.

Liczba wierszy definiowana jest za pomocą parametru tlags (w przykładzie tlags = 0:5). Analizując mapę dla kolejnych opóźnień w czasie (tlag od 0 do 5), można zaobserwować, jak podobieństwo danych zanika z każdym kolejnym krokiem czasowym (np. o ile mniej podobna jest jakość powietrza dzisiaj do tego sprzed 5 dni w tej samej lokalizacji). Liczba kolumn odpowiada kolejnym odstępom odległości (i jest równa cutoff/width).

Kolory na mapie odpowiadają wartościom semiwariancji (\(\gamma\)) dla określonego odstępu wyrażonego odległością i czasem. Im kolor jest cieplejszy, tym większa jest semiwariancja, czyli mniejsze podobieństwo między pomiarami. Miejsce, w którym kolory przestają stawać się coraz cieplejsze i stabilizują się na jednym poziomie, wyznacza tzw. zasięg (range). Poza tą odległością/czasem dane nie wykazują już autokorelacji.

Należy także zwrócić uwagę na szybkość (gradient) zmian. Jeśli kolory zmieniają się bardzo szybko wzdłuż osi X (reprezentującą wymiar przestrzenny), ale wolno wzdłuż osi Y oznacza to, że cecha zmienia się gwałtownie w przestrzeni, ale jest dość stabilna w czasie. Jeśli jest odwrotnie - kolory zmieniają się szybciej wzdłuż osi Y (reprezentującej wymiar “czasowy”) to cecha nie wykazuje znaczących zmian w przestrzeni (między stacjami pomiarowymi), ale dynamicznie zmienia się w czasie (np. nagły skok ciśnienia atmosferycznego na dużym obszarze).

Jak wygląda tempo zmian dla PM10? Czy zmienia się szybciej w czasie czy w przestrzeni? Co to oznacza?

3.2 Semiwariogramy czasoprzestrzenne

Argument map = FALSE w funkcji plot() pozwala na wyświetlenie serii semiwariogramów empirycznych. Każda linia reprezentuje semiwariogram przestrzenny dla konkretnego przesunięcia czasowego (\(u\)). Najniższa linia (zwykle \(u=0\)): Pokazuje “czystą” zmienność przestrzenną. Każda kolejna linia powyżej pokazuje, jak rośnie niepodobieństwo danych, gdy dodajemy do nich wymiar czasowy.

plot(vv, map = FALSE)

3.3 Semiwariogram “przestrzenny”

Semiwariogram dla opóźnienia czasu \(u = 0\) (lag0) pokazuje wyłącznie zmienność zależną od położenia w przestrzeni. Wynikowy obiekt działania funkcji variogramST() dostarcza informację dla którego opóźnienia czasowego zostały wykonane obliczenia. Informacja ta znajduje się w kolumnie id.

v_space = vv[vv$id == "lag0",]

vs = ggplot(v_space, aes(x = dist, y = gamma)) + geom_line() + geom_point() + 
  labs(x = "Odstęp [km]", y = "Semiwariancja", title = "Semiwariogram dla lag0:\n wymiar przestrzenny") + 
  geom_text(aes(label = np), 
            nudge_x = 5, 
            hjust = 0, 
            size = 3) +
  ylim(0, ceiling(max(v_space$gamma))) + 
  theme_bw() 

3.4 Semiwariogram “czasowy”

Semiwariogram dla odstępu odległości \(h = 0\) (dist = 0) pokazuje wyłącznie zmienność zależną od położenia w czasie. Wynikowy obiekt działania funkcji variogramST() dostarcza informację dla którego odstępu \(h\) zostały wykonane obliczenia. Informacja ta znajduje się w kolumnie dist.

v_time = vv[vv$dist == 0,]

vt = ggplot(v_time, aes(x = timelag, y = gamma)) + geom_line() + geom_point() + 
  labs(x = "Odstęp [dni]", y = "Semiwariancja", title = "Semiwariogram dla h = 0:\n wymiar czasowy") + 
  geom_text(aes(label = np), 
            nudge_x = 0.2, 
            hjust = 0, 
            size = 3) +
  ylim(0, ceiling(max(v_time$gamma))) + 
  theme_bw()

library(patchwork)
vs + vt

4 Budowa modelu teoretycznego oraz optymalizacja parametrów modelu

Istnieje kilka modeli teoretycznych stosowanych do analizy danych czasoprzestrzennych, które różnią się między sobą sposobem uwzględniania interakcji między wymiarem przestrzennym i czasowym. W poniższym przykładzie zostanie użyty tzw. model separowalny (ang. separable model).

Jest to jeden z najbardziej podstawowych i często używanych modeli teoretycznych. Model ten zakłada, że kowariancja czasoprzestrzenna \(C_{sep}(h, u)\) jest iloczynem niezależnych struktur przestrzennych (\(C_s(h)\)) i czasowych (\(C_t(u)\)).

\[C_{sep}(h, u) = C_s(h) \cdot C_t(u)\]

gdzie \(h\) oznacza odległość w przestrzeni, a \(u\) opóźnienie w czasie. Wzór na model separowalny semiwariogramu czasowoprzestrzennego, stosowany w pakiecie gstat (funkcja vgmST() z opcją “separable”), przyjmuje następującą postać:

\[\gamma_{\rm sep}(h,u) = {\rm sill} \cdot \left( \bar{\gamma}_s(h)+\bar{\gamma}_t(u)-\bar{\gamma}_s(h)\bar{\gamma}_t(u) \right)\]

gdzie:

• \(h\): odległość w przestrzeni (spacelag),

• \(u\): opóźnienie w czasie (timelag),

• \(\text{sill}\): całkowita wariancja progowa modelu (parametr skalujący),

• \(\bar\gamma_s(h)\): znormalizowany semiwariogram przestrzenny (o progu (sill) równym 1),

• \(\bar\gamma_t(u)\): znormalizowany semiwariogram czasowy (o progu (sill) równym 1).

Dla modelu separowalnego funkcja vgmST() wymaga zdefniowania następujących parametrów:

vgmST("separable", space = spatialVgm, time = temporalVgm, sill = sill)

• “separable”: tym modelu teoretycznego
• space: parametry modelu dla struktury przestrzennej
• time: parametry modelu dla struktury czasowej
• sill: łączna wartość wariancji progowej dla struktury czasowej i przestrzennej.

Parametry struktury przestrzennej oraz czasowej podaje się wykorzystując funkcję vgm().

Model separowalny wymaga podania następujących parametrów:

• Parametry dla struktury przestrzennej (argument space)

  • Typ modelu, np. Sph, Exp, Gau
  • Zasięg przestrzenny (range dla struktury )
  • Proporcjonalny nugget przestrzenny

• Parametry dla struktury czasowej (argument time)

  • Zasięg czasowy
  • Proporcjonalny nugget czasowy

• Wspólny próg całkowity dla struktury przestrzenej i czasowej (argument sill)

Aby ustalić początkowe parametry modelu należy wyświetlić semiwariogram dla struktury przestrzennej (wykonany dla zerowego opóźnienia czasowego: id = lag0) oraz semiwariogram dla struktury czasowej (wykonany dla \(h = 0\); dist = 0).

Dla semiwariogramu przestrzennego widzimy, że: całkowita wartość progowa (sill) = 50, nugget = 10, zasięg = 120 km. Aby znormalizować parametr sill i nugget, tak aby wynosiły 1, obliczamy: 10/50 = 0,2 oraz 40/50 = 0,8.

Łączna wartość progowa (argument sill) wynosi ok. 100. Aby ustalić tą wartość można sprawdzić, jaka jest maksymalna wartość w kolumnie gamma w obiekcie vv zawierającym wynik działania funkcji variogramST().

#separable model: spatial and temporal sill will be ignored
#and kept constant at 1-nugget respectively. A joint sill is used.
separableModel <- vgmST("separable", 
                        space=vgm(psill = 0.8,"Exp", range = 125, nugget = 0.2),
                        time =vgm(psill = 0.8,"Exp", range = 3, nugget = 0.2),
                        sill=100)

Funkcja plot() pozwala na porównanie semiwariogramu empirycznego z modelem.

plot(vv, separableModel, all = TRUE)

Następnie używamy funkcji fit.StVariogram() w celu automatycznego zoptymalizowania parametrów modelu. Argument lower oraz upper nakłada ograniczenia na zakres parametrów modelu: zasięg przestrzenny, nugget przestrzenny, zasięg czasowy, nugget czasowy, łączna wartość progowa - sill.

fitseparableModel <- fit.StVariogram(vv, separableModel,
                                  method="L-BFGS-B",
                                  lower=c(10,0,0.01,0,1),
                                  upper=c(1000,1,20,1,200))

Funkcja fitStVariogram() dopasowała następujące parametry:

fitseparableModel

space component: 
  model     psill    range
1   Nug 0.1410533   0.0000
2   Exp 0.8589467 478.9536
time component: 
  model psill    range
1   Nug     0 0.000000
2   Exp     1 3.042017
sill: 113.299900478399

Dla porównania wartości początkowe wynosiły:

separableModel

space component: 
  model psill range
1   Nug   0.2     0
2   Exp   0.8   125
time component: 
  model psill range
1   Nug   0.2     0
2   Exp   0.8     3
sill: 100

plot(vv, fitseparableModel, all = TRUE)

Dla porównania ponownie wyświetlamy semiwariogramy empiryczne:

plot(vv, map = FALSE)

5 Kriging

Do predykcji wykorzystuje się funkcje krigeST(). Kriging wymaga stworzenia obiektu z lokalizacjami do predykcji (siatka, punkty) w formacie klasy ST (lokalizacje określone w przestrzeni i czasie). Kriging wykonuje się w siatce interpolacyjnej dla wybranych momentów w czasie.

5.1 Tworzenie siatki interpolacyjnej

# Sprawdzenie zakresu przestrzennego danych 
bb = r5to10@sp@bbox
bb

               min      max
coords.x1  6.28107 14.78617
coords.x2 47.80847 54.92497

# --- 1. Przygotowanie siatki ---
# Wybór podzbioru danych dla konkretnego zakresu dat. 
rr <- r5to10[,"2005-03-20/2005-03-25"]

# Konwertujemy obiekt na STSDF (Spatio-Temporal Sparse Data Frame). 
# Jest to format często wymagany przez funkcje obliczeniowe pakietu gstat.
rr <- as(rr, "STSDF")

# --- 2. Tworzenie siatki przestrzennej (Grid) ---
# Definiujemy zakresy współrzędnych geograficznych (X i Y) dla siatki interpolacyjnej.
x1 <- seq(from = bb[1,1], to = bb[1,2], by = 0.25)   # Długość geograficzna (longitude)
x2 <- seq(from = bb[2,1], to = bb[2,2], by = 0.25)  # Szerokość geograficzna (latitude)

# Tworzymy punkty przestrzenne na podstawie wszystkich kombinacji x1 i x2.
# Przypisujemy im ten sam układ współrzędnych (CRS), który posiada nasz zbiór danych.
DE_gridded <- SpatialPoints(cbind(rep(x1, length(x2)), rep(x2, each = length(x1))), 
                            proj4string = CRS(proj4string(rr@sp)))

# Informujemy R, że te punkty tworzą regularną siatkę (grid), co przyspieszy obliczenia
gridded(DE_gridded) <- TRUE

# --- 3. Przygotowanie struktury do predykcji (newdata) ---
# Tworzymy obiekt STF (Space-Time Full grid), który łączy naszą nową siatkę przestrzenną z wymiarem czasowym (tymi samymi dniami, które wybraliśmy do rr). Obiekt STF przechowuje samą geometrię, bez atrybutów (ramki danych)
DE_pred <- STF(sp = as(DE_gridded, "SpatialPoints"), time = rr@time)

str(DE_pred)

Formal class 'STF' [package "spacetime"] with 3 slots
  ..@ sp     :Formal class 'SpatialPoints' [package "sp"] with 3 slots
  .. .. ..@ coords     : num [1:1015, 1:2] 6.28 6.53 6.78 7.03 7.28 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..@ bbox       : num [1:2, 1:2] 6.16 47.68 14.91 54.93
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"
  .. .. .. .. ..$ comment: chr "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.25722"| __truncated__
  ..@ time   :An xts object on 2005-03-20 / 2005-03-25 containing: 
  Data:    integer [6, 1]
  Columns: timeIndex
  Index:   Date [6] (TZ: "UTC")
  ..@ endTime: POSIXct[1:6], format: "2005-03-21" "2005-03-22" ...

5.2 Estymacja

Po utworzeniu obiektu z przestrzennymi i czasowymi lokalizacjami do predykcji można wykonać kriging:

# --- 4. Interpolacja Krigingiem Czasoprzestrzennym ---
# Wykonujemy kriging czasoprzestrzenny dla zmiennej PM10.
# data: dane źródłowe ograniczone do kilku dni dla których wykonujemy predykcję, newdata: siatka predykcyjna, modelList: dopasowany model wariogramu.
DE_kriged <- krigeST(PM10 ~ 1, data = rr, newdata = DE_pred, 
                     modelList = fitseparableModel)

5.3 Wizualizacja wyników

# --- 5. Wizualizacja wyników ---
# Ponownie upewniamy się, że wynikowy komponent przestrzenny jest traktowany jako siatka.
gridded(DE_kriged@sp) <- TRUE

# Rysujemy mapy interpolacji dla każdego kroku czasowego.
stplot(DE_kriged)

Mapy przedstawiają wyniki predykcji dla kilku wybranych dni w marcu.

Dla porównania rozkład wartości w tych dniach w stacjach pomiarowych przedstawia poniższa rycina:

stplot(r5to10[, "2005-03-20::2005-03-25"])

Wizualizację można rozszerzyć o dodatkowe elementy, np. granicę regionów i lokalizację stacji pomiarowych. Można także zmienić domyślną paletę.

# 1. Zdefiniuj układ dodatkowych warstw (sp.layout)
# Każdy element to lista: c("funkcja_rysująca", nazwa_obiektu, parametry...)
my_layout <- list(
  list("sp.polygons", DE_NUTS1, first = FALSE, col = "grey50", lwd = 1),
  list("sp.points", stations, pch = 3, col = "black", cex = 0.5)
)


# 2. Użyj stplot z argumentem sp.layout
stplot(DE_kriged, 
       sp.layout = my_layout,
       main = "Predykcja stężeń pyłów",
       col.regions = terrain.colors(100))

Można także zastosować dowolną paletę zdefiniowaną w pakiecie RColorBrewer.

library(RColorBrewer)
pal = colorRampPalette(brewer.pal(9, "YlOrRd"))(100)
stplot(DE_kriged, 
       sp.layout = my_layout,
       main = "Predykcja stężeń pyłów",
       col.regions = pal)

#rev odwraca kolejność kolorów. 
library(RColorBrewer)
pal = colorRampPalette(brewer.pal(9, "Spectral"))
pal = rev(pal(100))
stplot(DE_kriged, 
       sp.layout = my_layout,
       main = "Predykcja stężeń pyłów",
       col.regions = pal)