Rauman 3: Übung A

Für die heutige Übung A benötigen wir folgende 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).

Als erstes wollen wir für den Datensätze rotmilan Dichteschätzungen durchführen. Lade vorgängig die dafür notwendigen Packages in deine R-Session.

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

library("spatstat.geom")
library("spatstat.explore")

Aufgabe 1: Rotmilan Bewegungsdaten visualisieren

Die erste Frage, die bei solchen Bewegungsstudien typischerweise gestellt wird, lautet: Wo hält sich das Tier hauptsächlich auf? Um diese Frage zu beantworten, kann man als erstes einfach die Datenpunkte in einer einfachen Karte visualisieren. Erstellt zur Beantwortung dieser Frage nachstehende Karte.

Aufgabe 2: Kernel Density Estimation berechnen

In einer ersten Annäherung funktioniert dies, doch wir sehen hier ein klassisches Problem des “Overplotting”. Das heisst, dass wir durch die Überlagerung vieler Punkte in den dichten Regionen nicht abschätzen können, wie viele Punkte dort effektiv liegen und ggf. übereinander liegen. Es gibt hier verschiedene Möglichkeiten, die Punktdichte klarer zu visualisieren. Eine unter Biologen sehr beliebte Methode ist die Dichteverteilung mit einer Kernel Density Estimation (KDE). Dies v.a. darum, weil mit KDE das Habitat (Streifgebiet) eines Tieres abgeschätzt werden kann. Homeranges werden oft mit KDE95 und Core Areas mit KDE50 definiert (Fleming C., Calabrese J., 2016).

Für die Berechnung der Density verwenden wir die Funktion density.ppp aus spatstat. Diese Library ist etwas komplex in der Anwendung. Damit wir dieses Verfahren aber dennoch auf unsere Rotmilan-Daten anwenden können, haben wir eine eigene KDE-Funktion erstellt, die wir Euch zur Verfügung stellen.

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 spatstat 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_kde <- function(points, cellsize, bandwith, extent = NULL){
  
  points_ppp <- as.ppp(points) # konvertiert sf > ppp

  if(!is.null(extent)){
    # falls ein extent angegeben ist, wird dieser verwendet
    # um das "Beobachtungsfenster" zu setzen
    Window(points_ppp) <- as.owin(st_bbox(extent))
  } 

  # macht eine Dichteschätzung
  points_density <- density.ppp(x = points_ppp, sigma = bandwith, eps = cellsize)

  # konvertiert den Output in ein DataFrame
  points_density_df <- as.data.frame(points_density)

  points_density_df
}

Die Parameter der Funktion sollten relativ klar sein:

  • points: Ein Punktdatensatz aus der Class sf
  • cellsize: Die Zellgrösse des Output-Rasters
  • bandwith: Der Suchradius für die Dichteberechnung
  • extent (optional): Der Perimeter, in dem die Dichteverteilung berechnet werden soll. Wenn kein Perimeter angegeben wird, wird die “bounding box” von points genutzt.

Wenn wir nun mit my_kde() die Dichteverteilung berechnen, erhalten wir ein data.frame mit X und Y Koordinaten sowie eine Spalte value zurück. Nutzt diese drei Spalten mit geom_raster(), um eure Daten mit ggplot zu visualisieren (aes(x = X, y = Y, fill = value).

rotmilan_kde <- my_kde(points = rotmilan, cellsize = 1000, bandwith = 10000, extent = schweiz)

head(rotmilan_kde)
##         x       y        value
## 1 2485909 1075767 5.706506e-24
## 2 2485909 1076766 8.289075e-23
## 3 2485909 1077764 3.029525e-23
## 4 2485909 1078763 6.521282e-23
## 5 2485909 1079761 9.598037e-23
## 6 2485909 1080760 1.182799e-22

Die Kernel Density Estimation ist nun sehr stark von den tiefen Werten dominiert, da die Dichte in den meisten Zellen unseres Untersuchungsgebiets nahe bei Null liegt. Wie erwähnt sind Wissenschaftler häufig nur an den höchsten 95% der Werte interessiert. Folge folgenden Schritten, um das Resultat etwas besser zu veranschaulichen:

  1. Berechne die 95. Perzentile aller Werte mit der Funktion quantile und bennene diesen q95
  2. Erstelle eine neue Spalte in rotmilan_kde, wo alle Werte tiefer als q95 NA entsprechen
  3. (Optional): Transformiere die Werte mit log10, um einen differenzierteren Farbverlauf zu erhalten

Aufgabe 3: Dichteverteilung mit Thiessen Polygonen

Thiessen Polygone bieten eine spannende Alternative, um Unterschiede in der Dichteverteilung von Punktdatensätzen zu visualisieren. Wir wollen dies nun ausprobieren und konstruieren zum Schluss die Thiessen Polygone für die Rotmilan-Daten für das Untersuchungsgebiet Schweiz. Nutze die Anleitung aus der Aufgabe 2: Interpolation mit Nearest Neighbour für das Erstellen von Thiessen Polygonen.

Abbildung 21.1: Wenn wir jetzt die Thiessen Polygone (ohne Punkte) darstellen, wird deutlicher, wie die Dichteverteilung im Innern des Clusters aussieht.