1 Legende

Die folgende Notation wird in diesem Dokument benutzt:

ein_paar_R_befehle() #ein Kommentar dazu, wird von R ignoriert

## [1] "Die Ausgabe eines R-Befehls"

Aufgaben

Berlin (Quelle: Carlos Cilleruelo auf towardsdatascience.com)

Das Arbeiten mit Geodaten ist eine zentrale Komponente der Umweltdatenverbeitung. Außerhalb Geographischer Informationssysteme (GIS) gibt es mächtige Werkzeuge, um direkt aus R heraus Geodaten zu analysieren und zu visualisieren. Im Gegensatz zur Erstellung einfacher Karten im GIS ist das zwar mit verhältnismäßig viel Code verbunden. Dafür hat man aber in R mehr Flexibilität und kann die Erstellung von Karten direkt in weitergehende Analysen einbetten. Auch die automatisierte Erstellung von Karten geht mit R leichter.

In diesem Kapitel werden wir Dir im Wesentlichen “Rezepte” vorstellen, wie Du unterschiedliche Typen von Geodaten in R laden und damit arbeiten kannst. In diesem Kontext werden wir auch ein paar grundlegende Begriffe aus den Bereichen GIS und Kartographie rekapitulieren.

2 Karten erstellen in R…eigentlich gar nichts Neues?

Du hast in diesem Kurs bereits einge Pakete und Funktionen für die visuelle Darstellung von Daten kennengelernt. Diese Pakete ermöglichen es, Punkte, Linien oder Flächen in ein Koordinationsystem zu plotten - und um nichts anderes handelt es sich bei einer Karte.

2.1 Beispiel: Niederschlagsstationskollektiv des Deutschen Wetterdienstes

Für das folgende Beispiel benötigen wir zunächst das Paket sf. Installiere es, falls nötig.

library(sf)

## Warning: package 'sf' was built under R version 4.4.2

## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE

Wir lesen die Metadaten der Niederschlagsstationen des DWD ein. Nettes Feature: Wir müssen die Daten gar nicht runterladen, sondern können direkt die Datei via https ansprechen.

Wir verwenden hier read.fwf statt read.table, weil in der entsprechenden Datei die Spalten feste Breite haben (im Gegensatz zur Nutzung von Trennzeichen). Das Argument fileEncoding kannst Du getrost ignorieren… das müssen wir angeben, damit read.fwf die Umlaute etc. kapiert.

fdwdgauges <- "https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/daily/more_precip/recent/RR_Tageswerte_Beschreibung_Stationen.txt"

dwdgauges <- read.fwf(fdwdgauges, c(5,9,9,15,12,10, 42, 41, 10), 
                     fileEncoding="ISO-8859-1", skip=2,header = FALSE,
                     col.names =  c("id", "from", "to", "elevation", 
                                    "lat", "lon","name","state","policy"))

head(dwdgauges)

##   id     from       to elevation     lat     lon
## 1  1 19120101 19860630       478 47.8413  8.8493
## 2  2 19510101 20061231       138 50.8066  6.0996
## 3  3 18910101 20110331       202 50.7827  6.0941
## 4  4 19510101 19791031       243 50.7683  6.1207
## 5  6 19821101 20241118       455 48.8361 10.0598
## 6  7 19510101 19960131       473 48.8140 10.1285
##                                         name
## 1  Aach                                     
## 2  Aachen (Kläranlage)                      
## 3  Aachen                                   
## 4  Aachen-Brand                             
## 5  Aalen-Unterrombach                       
## 6  Aalen-Unterkochen                        
##                                       state     policy
## 1 Baden-Württemberg                         Frei      
## 2 Nordrhein-Westfalen                       Frei      
## 3 Nordrhein-Westfalen                       Frei      
## 4 Nordrhein-Westfalen                       Frei      
## 5 Baden-Württemberg                         Frei      
## 6 Baden-Württemberg                         Frei

Wir sehen: Die Datei enthält die geographischen Koordinaten der Stationen (lon, lat). Damit können wir die Positionen auf eine Karte plotten.

plot(lat ~ lon, dwdgauges, asp = 1,
     main = "DWD Niederschlagsstationen")

Wow…ganz schön viele Stationen! Und die messen alle Niederschlag?!

Manipulieren wir doch mal die Karte, indem wir z.B. diejenigen Stationen in rot darstellen, die gar nicht mehr aktiv sind. Dazu müssen wir die Zeitkomponente der Metadaten als Datum interpretieren und anschließend den Dataframe mit Hilfe des to-Attributs filtern.

dwdgauges$fromdate <- as.Date(as.character(dwdgauges$from), format="%Y%m%d")
dwdgauges$todate <- as.Date(as.character(dwdgauges$to), format="%Y%m%d")
inactive <- dwdgauges$todate < "2021-07-01"

plot(lat ~ lon, dwdgauges[!inactive,], asp = 1,
     main = "DWD Niederschlagsstationen")
points(lat ~ lon, dwdgauges[inactive,], asp = 1, col = "red")

Aha…also doch nicht so viele? Aber klar: Stationen wurden über die Zeit verlegt, stillgelegt, neu gebaut und wieder stillgelegt. Wenn man die Daten verarbeiten möchte, muss aber jeder Standort eine eindeutige ID haben…da kommt dann eben über die Zeit was zusammen. Apropos Zeit…reichern wir doch mal die Karte mit einer Zeitreihe an: Zahl der Stationen von 1901 bis 2021.

Installiere und lade das Paket lubridate, bevor Du die nächsten Zeilen ausführst. Das macht das Arbeiten mit Datumsdaten leichter.

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# Heute ist welches Jahr?
endyear   <- year(today())

# Vektor mit Jahreszahlen von 1901-heute
years     <- 1901:endyear

# Preallokation eines Vektors in den wir das Zählergebnis für jedes 
# Jahr reinschreiben
numgauges <- rep(0,length(years))

# Um uns die Arbeit zu erleichtern 
# Man erinnere sich: Stay DRY, don’t get WET!
hasdatafun <- function(data,Y){
  hasdata = Y >= year(data$fromdate) & 
            Y <= year(data$todate)
}

# Anhand einer Schleife zählen wir nun für jedes Jahr die Anzahl der 
# Stationen
for (iter in 1:length(years)){
  hasdata = hasdatafun(dwdgauges,years[iter])
  numgauges[iter] <- sum(hasdata)
}

# Nun können wir die Daten plotten
par(mfrow = c(2,2))
hasdata <- hasdatafun(dwdgauges,1920)
plot(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = "1920", pch = 20)

hasdata <- hasdatafun(dwdgauges,1987)
plot(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = "1987",ylab = "", pch = 20)

hasdata <- hasdatafun(dwdgauges,endyear)
plot(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = as.character(endyear),ylab = "", pch = 20)

plot(years,numgauges,type = "l")

Fazit

Wir haben einen Geodatensatz eingelesen und Karten erstellt, ohne ein neues Werkzeug zu nutzen. Obendrein haben wir die Karten ganz lässig mit Zeitreihen angereichert… an sowas müsstest Du in einem GIS lange schrauben!

Trotzdem sind wir noch nicht zufrieden. Es wäre z.B. schön, die Landesgrenzen und ggf. Grenzen der Bundesländer mit darzustellen…

3 Was fehlt uns, um effizient schöne Karten zu gestalten?

Wir haben nun also mit Standardpaketen bereits wesentliche Elemente von Karten erstellt? Was fehlt uns noch?

Einlesen gängiger Geodaten in Vektor- und Rasterformat.
Koordinatentransformation
Effizientes Plotten spezifischer Geodatentypen- und Elemente.
Gängige GIS-Operationen

Aufgabe

Kläre (am besten mit Deinen Nachbar:innen) mal ein paar GIS-Grundbegriffe…

Was sind Vektordaten? Nennt zwei gängige Datenformate!

Was sind Rasterdaten? Nennt zwei gängige Datenformate!

Nennt mindestens drei typische GIS-Operationen mit Vektordaten.

Wozu braucht man Projektionen (oder Spatial Reference Systems)?

4 Raster- vs. Vektordaten

Die folgende Abbildung macht noch einmal deutlich, worum es sich bei Vektordaten handelt.

Vektordaten

(Quelle: Colin Williams, NEON)

Typische Datenformate für Vektordaten:

Shapefile (.shp),
geojson (.json)
GPS Exchange Format (.gpx)
Google Keyhole Markup Language (.KML, .KMZ)
…

Aufgabe
Diskutiere mit Deinem Nachbarn: Für welche Umweltdaten ist das Vektorformat geeignet, für welche nicht?

Vektordaten sind geeignet für:

Verwaltungeinheiten, politische Grenzen

Straßen, andere Infrastrukturelemente

Gebäudeumrisse

Gewässer (Flüsse, Seen, …)

Habitate

Messpunkte

…

Vektordaten sind ungeignet für:

Digitale Geländemodelle

Fernerkundungsdaten

Daten aus Klimamodellen

…

5 Lesen und plotten von Vektordaten mit `sf`

sf hält wesentliche Funktionen bereit, die wir uns in R für Vektordaten wünschen.

Beginnen wir doch damit, unserer Karte Landesgrenzen hinzuzufügen. Das Shapefile mit den Bundesländergrenzen haben wir direkt vom Bundesamt für Kartographie und Geodäsie heruntergeladen (Link).

library(sf)

f_bundeslaender <- "../data/shapefiles/VG250_LAN.shp"
bundeslaender   <- st_read(f_bundeslaender)

## Reading layer `VG250_LAN' from data source 
##   `D:\WDmatlab\GIT\umweltdv\06_geodata_maps\data\shapefiles\VG250_LAN.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 35 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 280371.1 ymin: 5235856 xmax: 921292.4 ymax: 6106244
## Projected CRS: ETRS89 / UTM zone 32N

head(bundeslaender)

## Simple feature collection with 6 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 280371.1 ymin: 5471422 xmax: 674202.3 ymax: 6101444
## Projected CRS: ETRS89 / UTM zone 32N
##   ADE GF BSG ARS AGS      SDV_ARS                 GEN                  BEZ IBZ
## 1   2  4   1  01  01 010020000000  Schleswig-Holstein                 Land  20
## 2   2  4   1  02  02 020000000000             Hamburg Freie und Hansestadt  22
## 3   2  4   1  03  03 032410001001       Niedersachsen                 Land  20
## 4   2  4   1  04  04 040110000000              Bremen     Freie Hansestadt  23
## 5   2  4   1  05  05 051110000000 Nordrhein-Westfalen                 Land  20
## 6   2  4   1  06  06 064140000000              Hessen                 Land  20
##   BEM NBD SN_L SN_R SN_K SN_V1 SN_V2 SN_G FK_S3 NUTS        ARS_0    AGS_0
## 1  --  ja   01    0   00    00    00  000     0  DEF 010000000000 01000000
## 2  --  ja   02    0   00    00    00  000     0  DE6 020000000000 02000000
## 3  --  ja   03    0   00    00    00  000     0  DE9 030000000000 03000000
## 4  --  ja   04    0   00    00    00  000     0  DE5 040000000000 04000000
## 5  --  ja   05    0   00    00    00  000     0  DEA 050000000000 05000000
## 6  --  ja   06    0   00    00    00  000     0  DE7 060000000000 06000000
##          WSK         DEBKG_ID RS       SDV_RS         RS_0
## 1 2012-02-01 DEBKGDL200000029 01 010020000000 010000000000
## 2 1974-01-01 DEBKGDL20000E6GD 02 020000000000 020000000000
## 3 2015-01-01 DEBKGDL20000E6EW 03 032410001001 030000000000
## 4 2010-01-01 DEBKGDL20000E0SF 04 040110000000 040000000000
## 5 2009-11-01 DEBKGDL20000E6GR 05 051110000000 050000000000
## 6 2015-01-01 DEBKGDL20000E3LK 06 064140000000 060000000000
##                         geometry
## 1 MULTIPOLYGON (((464810.7 61...
## 2 MULTIPOLYGON (((578219 5954...
## 3 MULTIPOLYGON (((479451.1 59...
## 4 MULTIPOLYGON (((466930.3 58...
## 5 MULTIPOLYGON (((477607.3 58...
## 6 MULTIPOLYGON (((534242 5721...

Puuh, 26 Spalten - ganz schön viele Attribute…naja, wir intessieren uns ja erstmal nur für die Ländergrenzen. Schnell mal gucken, ob die da sind…

Die Funktion st_geometry legt erstmal alle Attribute beiseite, so dass nur die Geometrie geplottet wird.

plot(st_geometry(bundeslaender), axes = TRUE)

Ok, sieht ja schon gut aus…sogar mit Hoheitsgebieten im Meer - wollten wir die haben? Egal.

Was aber bei den Achsen auffällt, ist die Größenordnung…irgendwas um 600.000 für die x-Achse und zwischen 5 und 6 Millionen für die y-Achse. Längen- und Breitengrade sind das jedenfalls nicht.

st_crs(bundeslaender)

## Coordinate Reference System:
##   User input: ETRS89 / UTM zone 32N 
##   wkt:
## PROJCRS["ETRS89 / UTM zone 32N",
##     BASEGEOGCRS["ETRS89",
##         ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
##             MEMBER["European Terrestrial Reference Frame 1989"],
##             MEMBER["European Terrestrial Reference Frame 1990"],
##             MEMBER["European Terrestrial Reference Frame 1991"],
##             MEMBER["European Terrestrial Reference Frame 1992"],
##             MEMBER["European Terrestrial Reference Frame 1993"],
##             MEMBER["European Terrestrial Reference Frame 1994"],
##             MEMBER["European Terrestrial Reference Frame 1996"],
##             MEMBER["European Terrestrial Reference Frame 1997"],
##             MEMBER["European Terrestrial Reference Frame 2000"],
##             MEMBER["European Terrestrial Reference Frame 2005"],
##             MEMBER["European Terrestrial Reference Frame 2014"],
##             MEMBER["European Terrestrial Reference Frame 2020"],
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]],
##             ENSEMBLEACCURACY[0.1]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["UTM zone 32N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
##         BBOX[38.76,6,84.33,12.01]],
##     ID["EPSG",25832]]

Bevor wir Dich ins offene Messer laufen lassen: Wir müssen also umprojizieren. Die Daten liegen im projizierten Koordinatensystem “UTM Zone 32” vor. Für unseren Zweck ist ein geographisches Koordinatensystem besser. Die “EPSG Nummer 4326” sollte man sich merken. Sie steht für das WGS84-Koordinatensystem.

bundeslaender2 <- st_transform(bundeslaender, crs = 4326)

Nun können wir unsere Stationsdaten zusammen mit den Bundesländern plotten.

plot(st_geometry(bundeslaender2), axes = TRUE, main = "Stationen im Jahr 1920")

hasdata <- hasdatafun(dwdgauges,1920)
points(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = as.character(endyear))

Oder wir projizieren die Stationsdaten um. Dazu müssen wir den Dataframe zunächst in ein sf-Objekt umwandeln.

dwdgaugessf <- st_as_sf(dwdgauges, coords = c("lon","lat"), crs = 4326)
dwdgaugessf <- st_transform(dwdgaugessf, crs = st_crs(bundeslaender))

plot(st_geometry(bundeslaender), axes = TRUE, main = "Stationen im Jahr 1920")

hasdata <- hasdatafun(dwdgauges,1920)
plot(st_geometry(dwdgaugessf[hasdata,]),
     main = as.character(endyear), add = TRUE)

6 Lesen und plotten von Rasterdaten mit `terra`

terra ist ein Paket zum Lesen, Schreiben und Verarbeiten von Rasterdaten. Eine ausführliche Dokumentation des Pakets findest Du hier. Ein weiteres, häufig verwendetes Paket ist stars (Link). Das Paket raster (Link) gibt es auch, es ist aber etwas veraltet, und terra die Nachfolgerin.

Wir werden terra in diesem Kurs v.a. zum Lesen eines Rasters verwenden. Unser Ziel ist es, unsere schöne Deutschlandkarte mit einer Topographie zu hinterlegen. Es gibt für die meisten Bundesländer mittlerweile frei verfügbare hochaufgelöste (1x1 m) Digitale Geländemodelle (DGM; engl. digital elevation model, DEM) . Außerdem gibt es hochauflösende globale DGM mit Auflösungen zwischen 20 und 30 Metern (u.a. SRTM, ASTER). Eine “optimale” Kombination von SRTM und ASTER auf europäischer Ebene gibt es bei der Europäische Umweltagentur EEA hier.

Für eine deutschlandweite Darstellung sind diese Produkte aber zu hoch aufgelöst und damit ineffizient (Speicher, Rechenbedarf). Wir nutzen daher ein Produkt, das auf eine Auflösung von 1x1 km für ganz Europa “resampled” wurde. Es beruht auf dem globalen Datensatz GTOPO30. Auch diesen Datensatz haben wir bei der EEA hier heruntergeladen.

## terra 1.7.83

ELEV <- rast('../data/raster/elevation1x1_new.tif')

Das erhaltene Objekt ELEV ist vom Typ raster und wir können es ein wenig interviewen, ohne überhaupt etwas in den Arbeitsspeicher geladen zu haben.

cat(paste("Beginn des Interviews am", now(),'\n'))

## Beginn des Interviews am 2024-11-19 16:34:34.604517

cat(paste("Wie heißen Sie?:\n", names(ELEV),'\n'))

## Wie heißen Sie?:
##  elevation1x1_new

cat(paste("Wie viele Zeilen haben Sie?\n", nrow(ELEV),'\n'))

## Wie viele Zeilen haben Sie?
##  4251

cat(paste("Wie viele Spalten haben Sie?\n", ncol(ELEV),'\n'))

## Wie viele Spalten haben Sie?
##  4901

cat(paste("Wie lautet Ihre Bounding Box?\n", ext(ELEV),'\n'))

## Wie lautet Ihre Bounding Box?
##  ext(2399843.1035138, 7300843.1035138, 1249900.3833563, 5500900.3833563)

Vor allem aber die Frage, die wir jedem Geodatensatz stellen müssen: Welche Projektion (“coordinate reference system”, CRS)?

cat(paste("Welches Koordinatensystem haben Sie?\n", crs(ELEV),'\n'))

## Welches Koordinatensystem haben Sie?
##  PROJCRS["User_Defined_Lambert_Azimuthal_Equal_Area",
##     BASEGEOGCRS["GCS_User_Defined",
##         DATUM["User_Defined",
##             ELLIPSOID["User_Defined_Spheroid",6378137,0,
##                 LENGTHUNIT["metre",1,
##                     ID["EPSG",9001]]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]],
##     CONVERSION["Lambert Azimuthal Equal Area",
##         METHOD["Lambert Azimuthal Equal Area",
##             ID["EPSG",9820]],
##         PARAMETER["Latitude of natural origin",52,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",20,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",5071000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",3210000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]]]

Wir haben es also mit laea, der Lambert Azimuthal Equal Area, also einer flächentreuen Projektion, zu tun. Nun gut, wir nehmen das einfach mal so hin.

Rasterdatenformate wie GeoTIFF sind in der Fernerkundung beliebt. Daher sind sie so ausgelegt, dass in einem Datensatz unterschiedliche Bänder (engl. “bands”) abgelegt werden können (klassischerweise z.B. rot, grün, blau - RGB). Im unserem Datensatz erwarten wir aber nur die Geländehöhe. Dies können wir verifizieren:

cat(paste("Wie viele Bänder haben Sie?\n", nlyr(ELEV),'\n'))

## Wie viele Bänder haben Sie?
##  1

Ok, werfen wir mal einen Blick auf die Werte, die in dem Raster gespeichert sind. Fragen wir doch nach dem größten Wert.

cat(paste("Welches ist Ihr maximaler Wert?\n", global(ELEV, fun = "max"),'\n'))

## Welches ist Ihr maximaler Wert?
##  187

Aufgabe

Stopp…die maximale Geländehöhe in Europa soll 187 Meter betragen? Überlege, was für ein Missverständnis/Problem hier vorliegen könnte. Lies ggf. die Metadaten (Link). Schlage eine pragmatische Lösung vor.
In den Metadaten finden wir, dass ein ‘z factor’ für den Datensatz genutzt wurde. Wir müssen also das Raster mit dem Wert 10 multiplizieren.
ELEV <- ELEV*10

6.1 Testplot

Die Darstellung der Rasterdaten könnte einfacher nicht sein:

plot(ELEV)

Soweit, so gut: Wir haben einen Rasterdatensatz mit terra eingelesen, wir kennen sein Koordinatensystem (CRS) und können den Datensatz mit plot plotten.

7 Synthese: Stationsmessnetz mit Ländergrenzen und Topographie plotten

Wir plotten alles in das CRS des Rasterdatensatzes, weil wir uns mit dem Thema Projektion bzw. Warpen von Rasterdaten nicht befassen wollen. Wir bringen also zunächst unsere Vektordaten (Bundesländer und Stationsmessnetz) in das CRS des Rasterdatensatzes.

bundeslaender_laea <- st_transform(bundeslaender, crs = crs(ELEV))
dwdgaugessf  <- st_transform(dwdgaugessf, crs = crs(ELEV))

Und weil die Topographie mit einem Hillshade immer schicker aussieht, nutzen wir noch das Hillshade-Raster, welches uns auf der gleich EEA-Seite zur Verfügung gestellt wird (Link). Das Hillshade-Raster plotten wir ebenfalls mit imshow, aber mit einer Tansparenz (Parameter alpha). Das Ganze machen wir zudem mit dem Paket tmap, einem Paket für die Erstellung thematischer Karten.

hillsh = rast('../data/raster/hillshade1x1/hillshade1x1.tif')

library(tmap)

## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

Da wir nicht ganz Europa plotten wollen, geht es schneller, wenn wir die beiden Rasterdaten “clippen”.

hillshc <- crop(hillsh,ext(bundeslaender_laea))/255
ELEVc <- crop(ELEV,ext(bundeslaender_laea))
ELEVc[ELEVc == 0] = NA

tmap_mode("plot")

## tmap mode set to plotting

tmap_options(show.messages = FALSE, show.warnings = FALSE)
## tmap mode set to plotting
tm_shape(ELEVc) +
    tm_raster("elevation1x1_new", palette = terrain.colors(255), 
              legend.show = FALSE) +
tm_shape(hillshc) +
    tm_raster("BinValues", palette = gray.colors(255), 
              alpha = 0.5, legend.show = FALSE) +  
tm_shape(bundeslaender_laea) +
    tm_borders("white", lwd = .5) +
tm_shape(dwdgaugessf[hasdata,]) +
    tm_symbols(col = "black", size = 0.3) +
tm_legend(show = FALSE) +
tm_compass() + 
tm_graticules() +
tm_layout(title = "DWD Niederschlagsstationen 1920", 
          title.bg.color = "white", title.bg.alpha = 0.7)

Alternativ können wir mit tmap auch eine interaktive Webkarte erstellen.

tmap_mode("view")
tm_shape(dwdgaugessf[hasdata,]) +
    tm_symbols(col = "black", size = 0.3)

8 Übung für die Coding-Werkstatt

Fertige auf Basis der heute behandelten Inhalte und Daten eine Karte für Dein Heimatbundesland an, in welcher folgende Informationslayer vorhanden sind:

Topographie
Landesgrenzen
Hauptstadt
wesentliche Fließgewässer
…

Suche Dir ggf. erforderliche Geodaten selbst aus dem Netz!

9 Weiteres Studium

Viele weitere Infos und Beispiele findest Du auf der Website RSpatial und r-spatial.

Projektionen

Quelle: https://xkcd.com/977/ Creative Commons Attribution-NonCommercial 2.5

Geodaten in R - mehr als nur eine Alternative zum GIS

Eine Karte sagt mehr als tausend Worte

Wolfgang Schwanghart

1 Legende

Aufgaben

2 Karten erstellen in R…eigentlich gar nichts Neues?

2.1 Beispiel: Niederschlagsstationskollektiv des Deutschen Wetterdienstes

3 Was fehlt uns, um effizient schöne Karten zu gestalten?

Aufgabe

4 Raster- vs. Vektordaten

Aufgabe

5 Lesen und plotten von Vektordaten mit `sf`

6 Lesen und plotten von Rasterdaten mit `terra`

Aufgabe

6.1 Testplot

7 Synthese: Stationsmessnetz mit Ländergrenzen und Topographie plotten

8 Übung für die Coding-Werkstatt

9 Weiteres Studium

Geodaten in R - mehr als nur eine Alternative zum GIS

Eine Karte sagt mehr als tausend Worte

Wolfgang Schwanghart

1 Legende

Aufgaben

2 Karten erstellen in R…eigentlich gar nichts Neues?

2.1 Beispiel: Niederschlagsstationskollektiv des Deutschen Wetterdienstes

3 Was fehlt uns, um effizient schöne Karten zu gestalten?

Aufgabe

4 Raster- vs. Vektordaten

Aufgabe

5 Lesen und plotten von Vektordaten mit sf

6 Lesen und plotten von Rasterdaten mit terra

Aufgabe

6.1 Testplot

7 Synthese: Stationsmessnetz mit Ländergrenzen und Topographie plotten

8 Übung für die Coding-Werkstatt

9 Weiteres Studium

5 Lesen und plotten von Vektordaten mit `sf`

6 Lesen und plotten von Rasterdaten mit `terra`