C Przykład analizy geostatystycznej

C.1 Analiza geostatystyczna

Analiza geostatystyczna jest złożonym procesem, często wymagającym sprawdzenia jakości danych i ich korekcji oraz wypróbowania wielu możliwości modelowania. Poniższy appendiks skupia się na pokazaniu przykładu uproszczonej analizy geostatystycznej, w której głównym celem jest estymacja średniej wartości temperatury.

C.2 Przygotowanie danych

Pierwszym krokiem analizy geostatystycznej jest załadowanie pakietów, które zostaną użyte. Brakujące pakiety można także załadować także w trakcie analizy geostatystycznej.

Kolejnym krokiem jest wczytanie danych oraz sprawdzenie ich jakości.

##       tavg             geom    
##  Min.   :-3.463   POINT  :190  
##  1st Qu.: 1.081   epsg:NA:  0  
##  Median : 1.647                
##  Mean   : 1.633                
##  3rd Qu.: 2.181                
##  Max.   : 9.158

Obiekt moje_punkty zawiera tylko jedną zmienną tavg, która ma być użyta do stworzenia estymacji. Warto zwizualizować rozkład wartości tej zmiennej w postaci histogramu oraz mapy (rycina C.1 i C.2).

Rozkład wartości zmiennej tavg.

Rycina C.1: Rozkład wartości zmiennej tavg.

Rozkład przestrzenny wartości zmiennej tavg.

Rycina C.2: Rozkład przestrzenny wartości zmiennej tavg.

Pozwala to na zauważenie, że w badanej zmiennej występują co najmniej dwie wartości odstające.

## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 349744.4 ymin: 320313.8 xmax: 349744.4 ymax: 320313.8
## CRS:            NA
## # A tibble: 1 x 2
##    tavg                geom
##   <dbl>             <POINT>
## 1  9.16 (349744.4 320313.8)
## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 348791.9 ymin: 317639.4 xmax: 348791.9 ymax: 317639.4
## CRS:            NA
## # A tibble: 1 x 2
##    tavg                geom
##   <dbl>             <POINT>
## 1 -3.46 (348791.9 317639.4)

Jedna z nich ma wartość ok. -3.5 °C i jest znacznie niższa od pozostałych, druga natomiast jest znacznie wyższa od pozostałych i ma wartość ok. 9.2 °C. Należy w tym momencie zastanowić się czy te wartości odstające są prawidłowymi wartościami, czy też są one błędne. W tej sytuacji, nie posiadając zewnętrznej informacji, bezpieczniej jest usunąć te dwa pomiary. Można to zrobić wyszukując id punktów za pomocą pakietu mapview.

Teraz id punktów można użyć do ich wybrania i zastąpienia potencjalnie błędnych wartości wartościami NA.

Te punkty nadal istnieją jednak w obiekcie moje_punkty. Można je usunąć korzystając z funkcji is.na oraz indeksowania:

Po tej zmianie powinno się po raz kolejny obejrzeć dokładnie dane w celu stwierdzenia, czy problem został naprawiony i czy nie występują dodatkowe sytuacje problemowe (rycina C.3 i C.4).

Rozkład wartości zmiennej tavg po usunięciu wartości odstających.

Rycina C.3: Rozkład wartości zmiennej tavg po usunięciu wartości odstających.

Rozkład przestrzenny wartości zmiennej tavg po usunięciu wartości odstających.

Rycina C.4: Rozkład przestrzenny wartości zmiennej tavg po usunięciu wartości odstających.

Można dodatkowo stworzyć chmurę semiwariogramu w celu wyszukania potencjalnych wartości lokalnie odstających (rycina C.5).

Chmura semiwariogramu zmiennej tavg.

Rycina C.5: Chmura semiwariogramu zmiennej tavg.

C.3 Tworzenie modeli semiwariogramów

Posiadając już poprawne dane można sprawdzić czy badane zjawisko wykazuje anizotropię przestrzenną poprzez stworzenie mapy semiwariogramu (rycina C.6).

Mapa semiwariogramu zmiennej tavg.

Rycina C.6: Mapa semiwariogramu zmiennej tavg.

Uzyskana mapa nie pozwala na jednoznaczne stwierdzenie kierunkowej zmienności podobieństwa badanej cechy, w związku z tym można skupić się na modelowaniu izotropowym. Kolejnym etapem jest stworzenie semiwariogramu oraz jego modelowanie. Optymalnie tworzy się więcej niż jeden model semiwariogramu, co pozwala na porównanie uzyskanych wyników i wybór lepszego modelu. Do tego przykładu zostały stworzone dwa modele semiwariogramu. Pierwszy z nich używa tylko zmiennej tavg oraz modelu ręcznego o wybranych parametrach (ryciny C.7, C.8, C.9, C.10).

Semiwariogram zmiennej tavg.

Rycina C.7: Semiwariogram zmiennej tavg.

Model semiwariogramu zmiennej tavg.

Rycina C.8: Model semiwariogramu zmiennej tavg.

Drugi model, oprócz zmiennej tavg, używa też wartości współrzędnych oraz modelu o parametrach zmodyfikowanych przez funkcję fit.variogram().

Semiwariogram zmiennej tavg uwzględniający współrzędne x i y.

Rycina C.9: Semiwariogram zmiennej tavg uwzględniający współrzędne x i y.

##   model      psill    range
## 1   Nug 0.01521914    0.000
## 2   Sph 0.73418974 1866.631
Model semiwariogramu zmiennej tavg uwzględniający współrzędne x i y.

Rycina C.10: Model semiwariogramu zmiennej tavg uwzględniający współrzędne x i y.

C.4 Ocena jakości semiwariancji

Aby porównać oba modele należy przyjąć metodę walidacji oraz współczynnik jakości estymacji. W tym przykładzie użyto kroswalidacji metodą LOO (funkcja krige.cv) oraz pierwiastek średniego błędu kwadratowego (RMSE) jako miarę jakości.

Porównanie dwóch wartości RMSE pozwala zdecydowanie stwierdzić, że drugi model charakteryzuje się lepszą jakością estymacji.

## [1] 6.193243
## [1] 0.587956

C.5 Stworzenie siatki

Przedostatnim krokiem jest utworzenie siatki do estymacji. Do niej zostaną wpisane wartości uzyskane z modelu semiwariogramu.

C.6 Stworzenie estymacji

Następnie nowo utworzona siatka może posłużyć do stworzenia estymacji (rycina C.11).

## [using universal kriging]
Estymacja i wariancja estymacji drugiego modelu zmiennej tavg.Estymacja i wariancja estymacji drugiego modelu zmiennej tavg.

Rycina C.11: Estymacja i wariancja estymacji drugiego modelu zmiennej tavg.

Wynikiem uproszczonej analizy geostatystycznej jest mapa estymowanych wartości temperatury dla całego badanego obszaru oraz mapa przedstawiająca wariancję estymacji temperatury.