Sitzung 5 Geodaten

5.1 Lernziele dieser Sitzung

Sie können…

  • Pipes benutzen
  • einfache dplyr-Befehle ausführen
  • Koordinaten visualisieren

5.2 Voraussetzungen

Wir laden erstmal tidyverse:

library(tidyverse)

5.3 Exkurs: Pipes

Teil vom tidyverse ist auch das Paket magrittr, das einen besonderen Operator enthält: %>%

Der Operator %>% heißt Pipe und setzt das Ergebnis der vorherigen Funktion als ersten Parameter in die nächste Funktion ein. Zur Veranschaulichung:

anzahl_buchstaben <- length(letters)
sqrt(anzahl_buchstaben)

…ist das gleiche wie…

sqrt(length(letters))

…ist das gleiche wie…

length(letters) %>%
  sqrt()

…ist das gleiche wie…

letters %>%
  length %>%
  sqrt()

So können beliebig viele Funktionen aneinandergereiht werden. Und mit -> kann eine Variable „in die andere Richtung“ zugewiesen werden

letters %>%
  length() %>%
  sqrt() %>%
  round() %>%
  as.character() ->
  my_var

Gerade bei komplizierteren Zusammenhängen wird der Code so oft lesbarer, weil die Logik von links nach rechts, bzw. von oben nach unten gelesen werden kann.

5.4 Daten importieren

Beim Open-Data-Portal der Stadt Frankfurt steht ein Baumkataster zur Verfügung.

Die Datei im CSV-Format (comma separated values) kann entweder heruntergeladen und durch klicken importiert werden, oder direkt über den Befehl:

baumkataster <- read_csv2("http://offenedaten.frankfurt.de/dataset/73c5a6b3-c033-4dad-bb7d-8783427dd233/resource/7a73520b-961a-4aad-a582-449e676c247c/download/gprojekteopen-datadatenamt-67datenbaumauswahl_veroffentlichung_4baumauswahl_veroffentlichung_4.csv")

5.5 Überblick verschaffen

Mit summary() lässt sich eine Zusammenfassung der Werte generieren:

summary(baumkataster)
##  Gattung/Art/Deutscher Name  Baumnummer           Objekt         
##  Length:118403              Length:118403      Length:118403     
##  Class :character           Class :character   Class :character  
##  Mode  :character           Mode  :character   Mode  :character  
##                                                                  
##                                                                  
##                                                                  
##    Pflanzjahr   Kronendurchmesser    HOCHWERT         RECHTSWERT    
##  Min.   :1645   Min.   : 2.000    Min.   :5545117   Min.   :463163  
##  1st Qu.:1970   1st Qu.: 4.000    1st Qu.:5550428   1st Qu.:472715  
##  Median :1982   Median : 6.000    Median :5552601   Median :475219  
##  Mean   :1979   Mean   : 6.688    Mean   :5552953   Mean   :475244  
##  3rd Qu.:1995   3rd Qu.: 9.000    3rd Qu.:5555165   3rd Qu.:478201  
##  Max.   :2017   Max.   :63.000    Max.   :5563639   Max.   :485361

Genauere Infos über diese Merkmale gibt es auf dem Datenportal.

5.6 Visualisieren

Wie in der letzten Lektion besprochen, lässt sich der Datensatz mit ggplot() visualisieren, z.B.:

ggplot(baumkataster, aes(x = Kronendurchmesser)) +
  geom_histogram()

Eine neue Messreihe lässt sich z.B. so errechnen:

alter <- 2020 - baumkataster$Pflanzjahr
head(alter)
## [1] 100 100 100 100 100 100

Der Befehl mutate() funktioniert sehr ähnlich, gibt aber den veränderten Datensatz zurück:

mutate(baumkataster, alter = 2020 - Pflanzjahr)
## # A tibble: 118,403 × 8
##    Gattung/Art/Deutscher …¹ Baumn…² Objekt Pflan…³ Krone…⁴ HOCHW…⁵ RECHT…⁶ alter
##    <chr>                    <chr>   <chr>    <dbl>   <dbl>   <dbl>   <dbl> <dbl>
##  1 Platanus x hispanica, G… 1       Acker…    1920       8  5.55e6 473366.   100
##  2 Platanus x hispanica, G… 2       Acker…    1920       8  5.55e6 473363.   100
##  3 Platanus x hispanica, G… 3       Acker…    1920       8  5.55e6 473360.   100
##  4 Platanus x hispanica, G… 4       Acker…    1920       8  5.55e6 473357.   100
##  5 Platanus x hispanica, G… 5       Acker…    1920       8  5.55e6 473354.   100
##  6 Platanus x hispanica, G… 6       Acker…    1920       8  5.55e6 473351.   100
##  7 Platanus x hispanica, G… 7       Acker…    1920       8  5.55e6 473348.   100
##  8 Platanus x hispanica, G… 8       Acker…    1920       8  5.55e6 473345.   100
##  9 Platanus x hispanica, G… 9       Acker…    1920       8  5.55e6 473343.   100
## 10 Platanus x hispanica, G… 10      Acker…    1920       8  5.55e6 473340.   100
## # … with 118,393 more rows, and abbreviated variable names
## #   ¹​`Gattung/Art/Deutscher Name`, ²​Baumnummer, ³​Pflanzjahr,
## #   ⁴​Kronendurchmesser, ⁵​HOCHWERT, ⁶​RECHTSWERT

Derselbe Befehl mit dem Pipe-Operator:

baumkataster %>%
  mutate(alter = 2020 - Pflanzjahr)
## # A tibble: 118,403 × 8
##    Gattung/Art/Deutscher …¹ Baumn…² Objekt Pflan…³ Krone…⁴ HOCHW…⁵ RECHT…⁶ alter
##    <chr>                    <chr>   <chr>    <dbl>   <dbl>   <dbl>   <dbl> <dbl>
##  1 Platanus x hispanica, G… 1       Acker…    1920       8  5.55e6 473366.   100
##  2 Platanus x hispanica, G… 2       Acker…    1920       8  5.55e6 473363.   100
##  3 Platanus x hispanica, G… 3       Acker…    1920       8  5.55e6 473360.   100
##  4 Platanus x hispanica, G… 4       Acker…    1920       8  5.55e6 473357.   100
##  5 Platanus x hispanica, G… 5       Acker…    1920       8  5.55e6 473354.   100
##  6 Platanus x hispanica, G… 6       Acker…    1920       8  5.55e6 473351.   100
##  7 Platanus x hispanica, G… 7       Acker…    1920       8  5.55e6 473348.   100
##  8 Platanus x hispanica, G… 8       Acker…    1920       8  5.55e6 473345.   100
##  9 Platanus x hispanica, G… 9       Acker…    1920       8  5.55e6 473343.   100
## 10 Platanus x hispanica, G… 10      Acker…    1920       8  5.55e6 473340.   100
## # … with 118,393 more rows, and abbreviated variable names
## #   ¹​`Gattung/Art/Deutscher Name`, ²​Baumnummer, ³​Pflanzjahr,
## #   ⁴​Kronendurchmesser, ⁵​HOCHWERT, ⁶​RECHTSWERT

So lassen sich auch hier verschiedene Befehle verknüpfen. filter() beschränkt den Datensatz auf Merkmalsträger, die den Kriterien entsprechen:

baumkataster %>%
  mutate(alter = 2020 - Pflanzjahr) %>%
  filter(alter > 30) ->
  alte_baeume

summary(alte_baeume)
##  Gattung/Art/Deutscher Name  Baumnummer           Objekt         
##  Length:73859               Length:73859       Length:73859      
##  Class :character           Class :character   Class :character  
##  Mode  :character           Mode  :character   Mode  :character  
##                                                                  
##                                                                  
##                                                                  
##    Pflanzjahr   Kronendurchmesser    HOCHWERT         RECHTSWERT    
##  Min.   :1645   Min.   : 2.000    Min.   :5545117   Min.   :463163  
##  1st Qu.:1960   1st Qu.: 6.000    1st Qu.:5550415   1st Qu.:472667  
##  Median :1974   Median : 8.000    Median :5552480   Median :475708  
##  Mean   :1966   Mean   : 8.503    Mean   :5552593   Mean   :475402  
##  3rd Qu.:1980   3rd Qu.:10.000    3rd Qu.:5554589   3rd Qu.:478539  
##  Max.   :1989   Max.   :35.000    Max.   :5563639   Max.   :485360  
##      alter       
##  Min.   : 31.00  
##  1st Qu.: 40.00  
##  Median : 46.00  
##  Mean   : 53.54  
##  3rd Qu.: 60.00  
##  Max.   :375.00

Schließlich ergibt das Streudiagramm von Koordinaten so eine art Karte:

ggplot(alte_baeume) +
    geom_point(size = 0.1, aes(x = RECHTSWERT, y = HOCHWERT))

Diesen Ansatz werden wir in der nächsten Lektion vertiefen.

5.7 Aufgaben

  1. Besuchen Sie https://pleiades.stoa.org/ - worum geht es hier?

  2. Finden Sie den kompletten aktuellen Datensatz für „locations“ als CSV-Datei.

  3. Importieren Sie ihn in R und weisen Sie dem Datensatz den Namen pleiades zu.

pleiades <- read_csv("http://atlantides.org/downloads/pleiades/dumps/pleiades-locations-latest.csv.gz")
  1. Finden Sie geeignete Werte für (einzelne) Längen- und Breitengrade im Datensatz.
pleiades$reprLong %>% # Längengrad
 summary()
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -19.666   8.059  17.036  19.943  30.847 111.078    7256
pleiades$reprLat %>% # Breitengrad
  summary()
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  -29.54   36.21   39.13   39.52   43.67   62.50    7256
  1. Plotten Sie die Koordinaten auf x- und y-Achse mit ggplot(). Was erkennen Sie?
ggplot(pleiades) +
  geom_point(aes(x = reprLong, y = reprLat))

  1. Halbieren Sie die Größe und setzen Sie den Alpha-Wert der Punkte auf 0,2.
ggplot(pleiades) +
  geom_point(aes(x = reprLong, y = reprLat), size = 0.5, alpha = 0.2)

  1. Bringen Sie die Grafik in die Mercator-Projektion.
ggplot(pleiades) +
  geom_point(aes(x = reprLong, y = reprLat), size = 0.5, alpha = 0.2) +
  coord_map("mercator")

  1. Schauen Sie sich diesen Befehl an:
map_data("world") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = long,
                               y = lat,
                               group = group)) +
    coord_quickmap(xlim = c(-8, 40),
                   ylim = c(26, 48))
  1. Versuchen Sie, jede einzelne Zeile nachzuvollziehen, indem Sie die entsprechenden Funktionen recherchieren.

  2. Führen Sie den Befehl aus.

  1. Ändern Sie die Farbe der Flächen in hellgrau.
map_data("world") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = long, y = lat, group = group),
                 fill    = "grey") +
    coord_quickmap(xlim = c(-8, 40),
                   ylim = c(26, 48))

  1. Wählen Sie einen Kartenausschnitt, auf dem Portugal, Ägypten, Irak und Frankreich komplett zu sehen sind.
map_data("world") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = long, y = lat, group = group),
                 fill    = "grey") +
    coord_quickmap(xlim = c(-8, 48),
                   ylim = c(22, 50))

  1. Plotten Sie auf diesem Hintergrund den Datensatz pleiades. Passen Sie dabei die Parameter so an, dass es Ihnen optisch zusagt.
map_data("world") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = long, y = lat, group = group),
                 fill    = "grey") +
    coord_quickmap(xlim = c(-8, 48), 
                   ylim = c(22, 50)) +
    geom_point(data    = pleiades,
               mapping = aes(x = reprLong, y = reprLat),
               color   = "blue",
               size    = 0.3,
               alpha   = 0.5)

  1. Wählen Sie für die Karte die Bonnesche Projektion mit Standardparallele bei 40°N.
map_data("world") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = long, y = lat, group = group),
                 fill    = "grey") +
    coord_map("bonne", 40,
              xlim = c(-8, 48),
              ylim = c(22, 50)) +
    geom_point(data    = pleiades,
               mapping = aes(x = reprLong, y = reprLat),
               color   = "blue",
               size    = 0.3,
               alpha   = 0.5)

  1. Entfernen Sie alle Achsenbeschriftungen.
map_data("world") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = long, y = lat, group = group),
                 fill    = "grey") +
    coord_map("bonne", 40,
              xlim = c(-8, 48),
              ylim = c(22, 50)) +
    geom_point(data    = pleiades,
               mapping = aes(x = reprLong, y = reprLat),
               color   = "blue",
               size    = 0.3,
               alpha   = 0.5) +
    theme(axis.title = element_blank(),
          axis.ticks = element_blank(),
          axis.text  = element_blank())

  1. (Achtung: extrem knifflig!) Bilden Sie diese Grafik nach, die die Orte geordnet nach ältestem Fund darstellt:
pleiades %>%
  select(long = reprLong,
         lat = reprLat,
         timePeriods) %>%
  mutate(oldest = str_extract(timePeriods, "[ACHRL]")) %>%
  mutate(oldest = factor(oldest, levels = c("A", "C", "H", "R", "L"))) ->
  pleiades_periods
 
map_data("world") %>%
  filter(region != "Antarctica") %>%
  ggplot() +
    geom_polygon(mapping = aes(long, lat, group = group),
                 fill    = "darkgreen") +
    coord_map("albers",
              parameters = c(16, 37),
              xlim       = c(-8, 40),
              ylim       = c(26, 48)) +
    geom_point(data    = pleiades_periods,
               mapping = aes(long, lat, color = oldest),
               size    = 0.2,
               alpha   = 0.3) +
    scale_color_brewer("Früheste Epoche",
                       breaks = c("A", "C", "H", "R", "L"),
                       labels = c("1000-550 v.d.Z.",
                                  "550-330 v.d.Z.",
                                  "330 v.d.Z. - 30 n.d.Z.",
                                  "30-300 n.d.Z.",
                                  "300-640 n.d.Z."),
                       palette = "YlOrBr") +
    theme(axis.title       = element_blank(),
          axis.ticks       = element_blank(),
          axis.text        = element_blank(),
          panel.background = element_rect(fill  = "darkblue"),
          panel.grid       = element_line(color = "blue")) +
    guides(colour = guide_legend(override.aes = list(alpha = 1,
                                                     size  = 8,
                                                     shape = 15)))