Kartenerstellung mit interpolierten Werte

Die Interpolation der Temperatur

Zuerst erstellen wir ein gleichmäßiges Grid was für für die Vorhersage der Temperatur brauchen und dann wandeln wir dieses grid und unsere Originaldaten in sog. SpatialPointsDataFrame Objekte, welche für die Interpolationsfunktion ebraucht wird. Anstelle des Ladens beider Pakete werden ich beim Aufruf der Funktion, das entsprechende Paket mitangeben, um Euch zu zeigen welche Funktion aus welchem Paket ist.

# Um den Code allgemeiner zu halten werden die Daten unter anderem Namen gespeichert:
dat <- temp_bot

# Für die weiteren Schritten müssen die Koordinaten als 'x' und 'y' bezeichnet werden
dat$x <- dat$long
dat$y <- dat$lat

# Den räumliche Ausdehnung der Stationen ermitteln (mit etwas Rand)
x_min <- min(dat$x) - min(dat$x) / 100
x_max <- max(dat$x) + max(dat$x) / 100
y_min <- min(dat$y) - min(dat$y) / 100
y_max <- max(dat$y) + max(dat$y) / 100

# Ein gleichmäßiges grid erstellen, basierend auf meinem definierten Samplinggebiet
grd <- expand.grid(
  x = seq(from = x_min, to = x_max, by = .1),
  y = seq(from = y_min, to = y_max, by = .1))

# The following converts dat into a SpatialPointsDataFrame object (sp package)
# Die folgende Funktion aus sp wandelt den grid in das sog. SpatialPointsDataFrame 
# Objekt (wird bereits automatisch abgespeichert!!) 
sp::coordinates(grd) <- ~x + y
sp::gridded(grd) <- TRUE
# Das gleiche nun mit den Originaldaten
sp::coordinates(dat) = ~x + y

Nun kommt die Interpolation mit der ìdw() Funktion aus dem gstat Paket. Diese Funktion erstellt das Modell und berechnet die vohergesagten Werte gleich mit.

idw_mod <- gstat::idw(formula = temp ~ 1, locations = dat,
  newdata = grd)  # apply idw model for the data

## [inverse distance weighted interpolation]

class(idw_mod)

## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"

Wir Ihr sehen könnt, ist das entstandene Objekt kein SpatialPointsDataFrame mehr sondern ein SpatialPixelsDataFrame. Dieses muss zurück in einen data frame umgewandelt werden, damit wir es mit ggplot plotten können:

idw_mod <- as.data.frame(idw_mod)  
head(idw_mod)

##          x       y var1.pred var1.var
## 1 13.14403 53.8565  5.823782       NA
## 2 13.24403 53.8565  5.825647       NA
## 3 13.34403 53.8565  5.830258       NA
## 4 13.44403 53.8565  5.837742       NA
## 5 13.54403 53.8565  5.848246       NA
## 6 13.64403 53.8565  5.862020       NA

Einbau der Interpolation in die Karte

Die Reihenfolge des Codes aus Schritt 1-3 muss leicht umgeändert werden, wenn die Interpolierten mit geplottet werden sollen:

• Zuerst wird ein leerer Plot erstellt, da wir in den geom_ Funktionen auf unterschiedliche Dateätze greifen will.
• Zu diesem Plot kommt dann ein geom_tile layer für die interpolierten Werte.
• Aschl. ändern wir die Farbskala der tiles bzw. Kacheln. Wichtig ist dass wir das fill aesthetic jetzt ansprechen wollen und nicht mehr color - daher nehmen wir auch scale_fill_gradient2 anstelle von scale_color_gradient2.
• Anschl. fügen wir die Weltkarte ein (als geom_polygon) und ändern das Koordinatensystem wie vorhin.
• Die Punkte könnte man wenn als hohle Kreise mit einzeichnen
• Zum Abschluss wird das Layout geändert:

p2 <- ggplot() + 
  geom_tile(data = idw_mod, 
    mapping = aes(x, y, fill = var1.pred)) +
  scale_fill_gradient2(midpoint = 6,
    low = "blue", mid = "yellow", high = "red", 
    na.value = "black", limits = c(0, 12)) +
  geom_point(data = temp_bot, aes(x = long, y = lat), shape = 21,
    colour = "grey80", alpha = 0.4) +
  geom_polygon(data = world, aes(x = long, y = lat, group = group),
    fill = "black", colour = "black") +
  coord_map("ortho",  xlim = c(12, 28), ylim = c(54,66)) +
  labs(fill = "Temp (°C)", y = "Längengrad", x = "Breitengrad",
    title = "Interpolierte Bodentemperatur im Februar 2015") +
  theme_bw() +
  theme(panel.grid = element_blank())
p2

Um die Abstufungen deutlicher zu machen, kann man auch die vorhergesagten Temperaturen runden und in einen Faktor umwandeln, um dann eine diskrete Skala zu verwenden:

idw_mod$var1.pred <- as.factor(round(idw_mod$var1.pred, 0))
p2 <- ggplot() + 
  geom_tile(data = idw_mod, 
    mapping = aes(x, y, fill = idw_mod$var1.pred)) +
  # hier muss nun die scale_fill_manual für diskrete Werte genommen werden;
  # und um eine ähnliche Farbskala zu verwenden können, wir die Matlab-Farben
  # nehmen: dafür gibt es eine Funktion im Paket 'colorRamps' in der man n
  # (Anzahl der Kategorien) spezifiziert:
  scale_fill_manual(values = colorRamps::matlab.like(11)) +
  # damit die hohen Zahlen weiterhin in der Legende oben stehen drehen wir die Skala um:
  guides(fill =  guide_legend(reverse = TRUE)) +
  geom_point(data = temp_bot, aes(x = long, y = lat), shape = 21,
    colour = "grey80", alpha = 0.4) +
  geom_polygon(data = world, aes(x = long, y = lat, group = group),
    fill = "black", colour = "black") +
  coord_map("ortho",  xlim = c(12, 28), ylim = c(54,66)) +
  labs(fill = "Temp (°C)", y = "Längengrad", x = "Breitengrad",
    title = "Interpolierte Bodentemperatur im Februar 2015") +
  theme_bw() +
  theme(panel.grid = element_blank())
p2

Et voila, fertig ist die Karte!