Veröffentlichungsdatum

4. Dezember 2023

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:

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 umschlagen 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 Function eingehend zu studieren und allenfalls ganz auf die Function zu verzichten und stattdessen direkt gstat zu verwenden. Wenn ihr mit unserer Function 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: Raeumliche Interpolation mit IDW

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

Tips:

  • Was für Distanzen bei maxdist Sinn machen, könnt ihr dem Output aus der G-Funktion (vorherige Übung) entnehmen
  • Wählt am Anfang eine etwas Konvervative (grosse) cellsize und verringert diesen nur wenn euer Rechner damit gut klar kommt
  • Da der Output aus der Interpolation im gleichen Koordinatenbezugssystem sind 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().
Code 
library("purrr")
library("tidyr")

plots <- lapply(1:4, function(idp){
  idw <- my_idw(groundtruth = luftqualitaet,column = "value",cellsize = 1000, nmax = Inf,maxdist = Inf,idp = idp,extent = schweiz)
  ggplot() +
    geom_raster(data = idw, aes(X,Y, fill = pred)) +
    geom_sf(data = schweiz, fill = NA) +
    geom_sf(data = luftqualitaet, size = 1, shape = 3, alpha = 0.3) +
    scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11,"RdYlBu")),limits = c(0, 60), na.value = NA) +
    labs(fill = "μg/m3", title = paste("Inverse Distance Power (IDP):",idp)) +
    theme_void() +
    theme(legend.position = "bottom", legend.title = element_blank(),
      legend.key.width = unit(0.10, 'npc'),
      legend.key.height = unit(0.02, 'npc'))
})
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]

cowplot::plot_grid(plotlist = plots)

ggsave("rauman/images/idw.png", height = 18, width = 20, units = "cm")

Abbildung 60.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 Thiessenpolygone 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.

Code 
luftqualitaet_union <- st_union(luftqualitaet)

thiessenpolygone <- st_voronoi(luftqualitaet_union)
Code 
ggplot() +
  geom_sf(data = schweiz) +
  geom_sf(data = thiessenpolygone, fill = NA)

st_voronoi hat die Thiessenpolygone etwas weiter gezogen als wir sie wollen. Dies ist allerdings eine schöne Illustration der Randeffekte von Thiessenpolygonen, 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. Auch hier braucht es zwei kleine Vorverarbeitungsschritte:

  1. wie vorher müssen wir die einzelnen Kantons-Polygone miteinander verschmelzen. Dies erreichen wir mit st_union(). Wir speichern den Output als schweiz, was als Resultat ein einzelnes Polygon der Schweizergrenze retourniert.
  2. für die Thiessen-Polygone machen wir genau das Umgekehrte: st_voronoi() liefert ein einzelnes Feature mit allen Polygonen, welches sich nicht gerne clippen lässt. Mit st_cast() wird die GEOMETRYCOLLECTION in Einzelpolygone aufgeteilt.
Code 
thiessenpolygone <- st_cast(thiessenpolygone)

thiessenpolygone_clip <- st_intersection(thiessenpolygone, schweiz)
Code 
ggplot() +
  geom_sf(data = schweiz) +
  geom_sf(data = thiessenpolygone_clip, fill = NA)

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).

Code 
thiessenpolygone_clip <- st_as_sf(thiessenpolygone_clip)
thiessenpolygone_clip <- st_join(thiessenpolygone_clip, luftqualitaet)
Code 
ggplot() +
  geom_sf(data = schweiz) +
  geom_sf(data = thiessenpolygone_clip, aes(fill = value)) +
  geom_sf(data = luftqualitaet) +
  scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11, "RdYlBu"))) +
  theme_void() +
  theme(
    legend.position = "bottom", legend.title = element_blank(),
    legend.key.width = unit(0.10, "npc"),
    legend.key.height = unit(0.02, "npc")
  )

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