5 Metody interpolacji

Odtworzenie obliczeń z tego rozdziału wymaga załączenia poniższych pakietów oraz wczytania poniższych danych:

Przez przejściem do interpolacji geostatystycznych warto zdać sobie sprawę, że nie jest to jedyna możliwa droga postępowania do stworzenia predykcji przestrzennych. Można wyróżnić dwie główne grupy modeli przestrzennych - modele deterministyczne oraz modele statystyczne.

5.1 Tworzenie siatek

W większości przypadków analiz geostatystycznych konieczne jest stworzenie siatki interpolacyjnej (pustego rastra). Istnieją dwa podstawowe rodzaje takich siatek - siatki regularne oraz siatki nieregularne.

5.1.1 Siatki regularne

Siatki regularne mają kształt prostokąta obejmującego cały analizowany obszar. Określenie granic obszaru można wykonać na podstawie zasięgu danych punktowych za pomocą funkcji bbox() z pakietu sp.

##        min      max
## x 745546.9 756937.4
## y 712618.8 721192.6

Do stworzenia siatki można użyć funkcji expand.grid(). Wymaga ona określenia dwóch argumentów - x oraz y (taka ich nazwa nie jest obowiązkowa). Oba argumenty przyjmują trzy wartości: (i) from oznaczający wartość początkową współrzędnej, (ii) to określający wartość końcową współrzędnej, oraz (iii) by określający rozdzielczość. Przy ustalaniu wartości początkowej i końcowej konieczne jest ich rozszerzenie względem wartości z funkcji bbox() lub extent(), aby wszystkie analizowane punkty znalazły się na badanym obszarze.

Utworzony w ten sposób obiekt wymaga określenia współrzędnych (funkcja coordinates()), potwierdzenia że dane mają być siatką (funkcja gridded()), oraz przypisania układu współrzędnych z obiektu punktowego (funkcja proj4string()).

Alternatywnie, do stworzenia siatki można wykorzystać funkcję makegrid(). Tworzy ona nowy obiekt na podstawie istniejącego obiektu punktowego oraz zadanej rozdzielczości.

5.1.2 Siatki regularne - wizualizacja

Sprawdzenie, czy uzyskana siatka oraz dane punktowe się na siebie nakładają można sprawdzić za pomocą funkcji plot(). W poniższym przykładzie, pierwszy wiersz służy wyświetleniu siatki, a drugi dodaje dane punktowe z użyciem argumentu add.

5.1.3 Siatki nieregularne

Siatki nieregularne mają zazwyczaj kształt wieloboku obejmującego analizowany obszar. Mogą one powstać, np. w oparciu o wcześniej istniejące granice. Siatki nieregularne w R mogą być reprezentowane poprzez klasę SpatialPixelsDataFrame.

W poniższym przypadku odczytywana jest granica badanego obszaru z pliku GeoPackage. Taki obiekt można np. stworzyć za pomocą oprogramowania gisowego takiego jak QGIS. Następnie tworzony jest nowy obiekt nowa_siatka_n poprzez wybranie tylko tych oczek siatki, które znajdują się wewnątrz zadanych granic.

## OGR data source with driver: GPKG 
## Source: "/home/travis/build/Nowosad/geostat_book/dane/granica.gpkg", layer: "granica.gpkg"
## with 1 features
## It has 3 fields

Wynik przetworzenia można zobaczyć z użyciem funkcji plot.

5.2 Modele deterministyczne

Modele deterministyczne charakteryzują się tym, że ich parametry są zazwyczaj ustalane w oparciu o funkcję odległości lub powierzchni. W tych modelach brakuje szacunków na temat oceny błędu modelu. Zaletą tych modeli jest ich prosta oraz krótki czas obliczeń. Do modeli deterministycznych należą, między innymi:

  • Metoda diagramów Voronoi’a (ang. Voronoi diagram)
  • Metoda średniej ważonej odległością (ang. Inverse Distance Weighted - IDW)
  • Funkcje wielomianowe (ang. Polynomials)
  • Funkcje sklejane (ang. Splines)

5.2.1 Voronoi

Metoda diagramów Voronoi’a polega na stworzeniu nieregularnych poligonów na podstawie analizowanych punktów, a następnie wpisaniu w każdy poligon wartości odpowiadającego punktu. Na poniższym przykładzie ta metoda stosowana jest z użyciem funkcji voronoi() z pakietu dismo.

5.2.2 IDW

Metoda średniej ważonej odległością (IDW) wylicza wartość dla każdej komórki na podstawie wartości punktów obokległych ważonych odwrotnością ich odległości. W efekcie, czym bardziej jest punkt oddalony, tym mniejszy jest jego wpływ na interpolowaną wartość. Wagę punktów ustala się z użyciem argumentu wykładnika potęgowego (ang. power). W pakiecie gstat istnieje do tego celu funkcja idw(), która przyjmuje analizowaną cechę (temp~1), zbiór punktowy, siatkę, oraz wartość wykładnika potęgowego (argument idp).

## [inverse distance weighted interpolation]

5.2.3 Funkcje wielomianowe

Stosowanie funkcji wielomianowych w R może odbyć się z wykorzystaniem funkcji gstat() z pakietu gstat. Wymaga ona podania trzech argumentów: formula określającego naszą analizowaną cechę (temp~1 mówi, że chcemy interpolować wartość temperatury zależnej od samej siebie), data określający analizowany zbiór danych, oraz degree określającą stopień wielomianu. Następnie funkcja predict() przenosi nowe wartości na wcześniej stworzoną siatkę.

## [ordinary or weighted least squares prediction]

## [ordinary or weighted least squares prediction]

## [ordinary or weighted least squares prediction]

5.2.4 Funkcje sklejane

Interpolacja z użyciem funkcji sklejanych (funkcja Tps() z pakietu fields) dopasowuje krzywą powierzchnię do wartości analizowanych punktów.

5.2.5 Porównanie modeli deterministycznych

5.3 Modele statystyczne

Modele statystyczne charakteryzują się tym, że ich parametry określane są w oparciu o teorię prawdopodobieństwa. Dodatkowo wynik estymacji zawiera także oszacowanie błędu, jednak te metody zazwyczaj wymagają większych zasobów sprzętowych. Do modeli statystycznych należą, między innymi:

  • Kriging
  • Modele regresyjne
  • Modele bayesowe
  • Modele hybrydowe

W kolejnych rozdziałach można znaleźć omówienie kilku podstawowych typów pierwszej z tych metod - krigingu.

5.4 Zadania

  1. Stwórz siatkę interpolacyjną o rozdzielczości 200 metrów dla obszaru Suwalskiego Parku Krajobrazowego.
  2. Korzystając z danych punkty wykonaj interpolację zmiennej srtm używając:
  • Poligonów Voronoi’a
  • Metody IDW
  • Funkcji wielomianowych
  • Funkcji sklejanych
  1. Porównaj uzyskane wyniki poprzez ich wizualizację. Czym różnią się powyższe metody?
  2. Wykonaj interpolację zmiennej temp metodą IDW sprawdzając różne parametry argumentu idp. W jaki sposób wpływa on na uzyskaną interpolację?