Rauman 3: Übung B

In dieser Übung geht es darum, zwei verschiedene Interpolationsverfahren in R umzusetzen. Im ersten Interpolationsverfahren verwenden wir die inverse distance weighted interpolation, später verwenden wir die nearest neighbour Methode. Dazu braucht ihr die folgenden Packages:

library("sf")
library("dplyr")
library("ggplot2")
library("gstat")

Weiter benötigt ihr die nachstehenden Datensätze:

Importiere diese in R. Prüfe, ob das CRS korrekt gesetzt wurde und setze es wenn nötig. Mache dich mit den Daten vertraut (visualieren, durchscrollen usw).

luftqualitaet <- read_sf("datasets/rauman/luftqualitaet.gpkg")
schweiz <- read_sf("datasets/rauman/schweiz.gpkg")

Die Library gstat bietet verschiedene Möglichkeiten, Datenpunkte zu interpolieren, unter anderem auch die inverse distance weighted Methode. Leider ist das Package noch nicht so benutzerfreundlich wie sf: Das Package wird aber aktuell überarbeitet und in mittlerer Zukunft sollte es ebenso einfach zugänglich sein. Damit Ihr Euch nicht mit den Eigenheiten dieser Library herumschlagen müsst, haben wir eine Function vorbereitet, die Euch die Verwendung der IDW-Interpolation erleichtern soll.

Wir nehmen Euch damit etwas Komplexität weg und liefern Euch ein pfannenfertiges Werkzeug. Das hat auch Nachteile und wir ermutigen alle, die dafür Kapazität haben, unsere Funktion eingehend zu studieren und allenfalls ganz auf die Funktion zu verzichten und stattdessen direkt gstat zu verwenden. Wenn ihr mit unserer Funktion arbeiten möchtet, müsst ihr den unten stehenden Code in euer Skript kopieren und ausführen.

my_idw <- function(groundtruth, column, cellsize, nmax = Inf, maxdist = Inf, idp = 2, extent = NULL){
  library("gstat")
  library("sf")
  
  if(is.null(extent)){
    extent <- groundtruth
  }
  
  samples <- st_make_grid(extent, cellsize, what = "centers")
  my_formula <- formula(paste(column,"~1"))
  idw_sf <- gstat::idw(formula = my_formula, groundtruth, newdata = samples, nmin = 1, nmax = nmax, maxdist = maxdist, idp = idp) 
  
  idw_matrix <- cbind(as.data.frame(st_coordinates(idw_sf)), pred = st_drop_geometry(idw_sf)[,1])
  idw_matrix
}

Nun könnt Ihr mit my_idw() den Datensatz luftqualitaet folgendermassen interpolieren.

my_idw(groundtruth = luftqualitaet, column = "value", cellsize = 10000, extent = schweiz)

Folgende Parameter stehen Euch zur Verfügung:

Beim Output handelt sich hier um einen Raster-ähnlichen Datentyp (siehe Vorlesung Spatial DataScience 1). Diesen können wir mit geom_raster mit ggplot visualisieren. Dafür müsst ihr in aes die X und Y Koordinaten angeben und der interpolierte Wert mit fill einfärben.

Aufgabe 1: Räumliche Interpolation mit IDW

Berechnet IDW-Interpolationen für die Luftqualitätsmessungen mit verschiedenen Parametern und visualisiert jeweils die Resultate. Experimentiert mit nmax sowie maxdist. Was stellt ihr fest?

Tips:

  • Wählt am Anfang eine etwas konservative (grosse) cellsize und verringert diesen nur, wenn euer Rechner damit gut klar kommt.
  • Da der Output aus der Interpolation im gleichen Koordinatenbezugssystem ist wie schweiz.gpkg, kann man diese beiden Datensätze im gleichen ggplot darstellen. Dafür müsst ihr die aesthetics (aes()) für jeden Layer einzeln setzen und nicht auf der Ebene von ggplot().
Abbildung 22.1: Stickstoffdioxid (NO2) in μg/m3, interpoliert über die ganze Schweiz mit der Inverse Distance Weighted Methode. Die verschiedenen Plots zeigen die Veränderung der Interpolation bei steigendem IDP-Wert

Aufgabe 2: Interpolation mit Nearest Neighbour

Eine weitere einfache Möglichkeit zur Interpolation bietet die Erstellung eines Voronoi-Diagrammes, auch als Thiessen Polygone oder Dirichlet-Zerlegung bekannt. sf liefert dazu die Funktion st_voronoi(), die einen Punktdatensatz annimmt und eben um die Punkte die Thiessen Polygone konstruiert. Dazu braucht es lediglich einen kleinen Vorverarbeitungsschritt: sf möchte für jedes Feature, also für jede Zeile in unserem Datensatz, ein Voronoidiagramm. Das macht bei uns wenig Sinn, weil jede Zeile nur aus einem Punkt besteht. Deshalb müssen wir vorher luftqualitaet mit st_union() von einem POINT- in ein MULTIPOINT-Objekt konvertieren, in welchem alle Punkte in einer Zeile zusammengefasst sind.

st_voronoi hat die Thiessen Polygone etwas weiter gezogen als wir sie wollen. Dies ist allerdings eine schöne Illustration der Randeffekte von Thiessen Polygonen, die zum Rand hin (wo es immer weniger Punkte hat) sehr gross werden können. Wir können die Polygone auf die Ausdehnung der Schweiz mit st_intersection() clippen. Vorher konvertieren wir mit st_collection_extract() die GEOMETRYCOLLECTION in POLYGON.

Jetzt müssen wir nur noch den jeweiligen Wert für jedes Polygon ermitteln. Dies erreichen wir wieder durch st_join. Auch hier ist noch ein kleiner Vorverarbeitungsschritt nötig: Wir konvertieren das sfc-Objekt (nur Geometrien) in ein sf-Objekt (Geometrien mit Attributtabelle).

Abbildung 22.2: Stickstoffdioxid (NO2) in μg/m3, interpoliert über die ganze Schweiz nach der Nearest Neighbour Methode.