Veröffentlichungsdatum

17. Oktober 2023

Aufgabe 1: Zähldaten

1a)

Die Projektstruktur steht. Nun können die Daten eingelesen und die nötigen Datentypen definiert werden.

Lädt die Daten zuerst von Moodle herunter.

Hinweise:

  • Siehe [Einführung] für den Standort der Zähler 211 und 502.

  • Alle für die Fallstudie Profil S benötigten Daten könnt ihr unter folgendem Link herunterladen.

  1. Zähldaten zu eurem Standort (211_sihlwaldstrasse_2017_2022.csv, 502_sihluferweg_2016_2022.csv)
  2. Meteodaten (order_105742_data.txt)

Die Zähldaten des WPZ wurden vorgängig bereinigt (z.B. wurden Stundenwerte entfernt, an denen am Zähler Wartungsarbeiten stattgefunden haben). Das macht es für uns einfach, denn wir können die Daten ohne vorgängige Bereinigung einlesen. Behaltet aber im Hinterkopf, dass die Datenaufbereitung, die Datenbereinigung mit viel Aufwand verbunden ist.

  • Lest die Zählaten ein, speichert ihn unter der Variable depo und sichtet den Datensatz (z.B. str(), head(), view() usw.).
Code 
# Speicherort sowie Dateiname anpassen
depo <- read_delim("./HIER RELATIVEN DATEIPFAD EINGEBEN", "HIER SEPERATOR EINGEBEN")

Hinweis: Im Stundenformat zeigen die Werte bei 11:00 die Zähldaten zwischen 11:00 bis 12:00 Uhr.

Code 
# lese die Daten ein
depo <- read_delim("datasets/fallstudie_s/WPZ/211_sihlwaldstrasse_2017_2022.csv", ";")

# erstes Sichten und anpassen der Datentypen
str(depo)

1b)

  • Nun muss das Datum als solches definiert werden. Welches Format hat es (im Code: format = “HIER DATUMSFORMAT”)?
Code 
depo <- depo |>
  mutate(
    Datetime = as.POSIXct(DatumUhrzeit, format = "HIER STEHT DAS DATUMSFORMAT", tz = "CET"),
    # nun schreiben wir uns das Datum in eine seperate Spalte
    Datum = as.Date(Datetime)
  )
Code 
# hier der code mit dem richtigen Format
depo <- depo |>
  mutate(
    Datetime = as.POSIXct(DatumUhrzeit, format = "%d.%m.%Y %H:%M", tz = "CET"),
    Datum = as.Date(Datetime)
  )

1c)

Ihr könnt selbst wählen, ob ihr Fussgänger:innen oder Fahrräder untersuchen wollt (je nachdem ob sie in eurem Datensatz vorhanden sind).

  • Entfernt die überflüssigen Spalten aus dem Datensatz. Ich schlage vor, dass ihr dafür den Befehl dplyr::select() verwendet.
Code 
# In dieser Auswertung werden nur Personen zu Fuss betrachtet!
# it select werden spalten ausgewaehlt oder eben fallengelassen
depo <- depo |>
  dplyr::select(-c(Velo_IN, Velo_OUT, Zeit, DatumUhrzeit))

1d)

  • Berechnen des Totals (IN + OUT), da dieses in den Daten nicht vorhanden ist (wiederum mit piping).

Tipp: Wenn man R sagt: “addiere mir Spalte x mit Spalte y”, dann macht R das für alle Zeilen in diesen zwei Spalten. Wenn man nun noch sagt: “speichere mir das Ergebnis dieser Addition in einer neuen Spalte namens Total”, dann hat man die Aufgabe bereits gelöst. Arbeitet mit mutate()).

Hinweis: Ihr habt das auch schon in Kapitel [Einführung und Installation] gemacht.

  • Entfernt nun alle NA-Werte mit na.omit().
Code 
# Berechnen des Totals, da dieses in den Daten nicht vorhanden ist
depo <- depo |>
  mutate(Total = Fuss_IN + Fuss_OUT)

# Entferne die NA's in dem df.
depo <- na.omit(depo)

Aufgabe 2: Meteodaten

2a)

  • Lest die Meteodaten ein und speichert sie unter meteo.
Code 
# Einlesen
meteo <- read_delim("datasets/fallstudie_s/WPZ/order_105742_data.txt", ";")

2b)

  • Auch hier müssen die Datentypen manuell gesetzt werden.

Tipp: Das Datum wird als Integer erkannt. Zuerst muss es in Text umgewandelt werden aus dem dann das eigentliche Datum herausgelesen werden kann. Das ist mühsam - darum hier der Code.

Code 
meteo <- mutate(meteo, time = as.Date(as.character(time), "%Y%m%d"))

Die Zeitangaben sind in UTC: 00:40 UTC = 02:40 Sommerzeit = 01:40 Winterzeit, Beispiel: 13 = beinhaltet Messperiode von 12:01 bis 13:00

–> Da wir mit Tageshöchstwerten oder -summen rechnen, können wir zum Glück ignorieren, dass das nicht mit den Daten der Zählstellen übereinstimmt. Learning: es ist zentral immer die Metadaten zu checken.

Hinweis Was ist eigentlich Niederschlag:

Link Meteo Schweiz

  • Werden den anderen Spalten die richtigen Typen zugewiesen? Falls nicht, ändert die Datentypen.

  • Nun schneiden wir den Datensatz auf die Untersuchungsdauer zu.

  • Dann müssen auch hier alle nicht verfügbare Werte (NA’s) herausgefiltert werden.

  • Prüft nun, wie die Struktur des data.frame (df) aussieht und ob alle NA Werte entfernt wurden (sum(is.na(df$Variable))). Stimmen alle Datentypen?

Code 
# Die eigentlichen Messwerte sind alle nummerisch
meteo <- meteo |>
    mutate(
        tre200nx = as.numeric(tre200nx),
        tre200jx = as.numeric(tre200jx),
        rre150n0 = as.numeric(rre150n0),
        rre150j0 = as.numeric(rre150j0),
        sremaxdv = as.numeric(sremaxdv)
    ) |>
    filter(time >= depo_start, time <= depo_end) # schneide dann auf Untersuchungsdauer

# Was ist eigentlich Niederschlag:
# https://www.meteoschweiz.admin.ch/home/wetter/wetterbegriffe/niederschlag.html

# Filtere Werte mit NA
meteo <- na.omit(meteo)
# Pruefe ob alles funktioniert hat
str(meteo)
sum(is.na(meteo)) # zeigt die Anzahl NA's im data.frame an

Aufgabe 3: Datenvorverarbeitung (Mutationen)

3a)

Jetzt fügen wir viele Convinience Variabeln hinzu. Wir brauchen:

  1. Wochentag; der Befehl dazu ist wday(). Danach als Faktor speichern.
  2. Werktag oder Wochenende als Faktor. Der Code dazu könnte so aussehen:
Code 
  ...|>
  mutate(Wochenende = ifelse(Wochentag %in% c(6,7), "Wochenende", "Werktag")) |>
  mutate(Wochenende = as.factor(Wochenende)) |>
  ...

je als Faktor: 3. Kalenderwoche: isoweek() 4. Monat: month() 5. Jahr: year()

Code 
depo <- depo |>
  # wday sortiert die Wochentage automatisch in der richtigen Reihenfolge
  mutate(
    Wochentag = wday(Datetime, week_start = 1),
    Wochentag = factor(Wochentag),
    # Werktag oder Wochenende hinzufuegen
    Wochenende = ifelse(Wochentag %in% c(6, 7), "Wochenende", "Werktag"),
    Wochenende = as.factor(Wochenende),
    # Kalenderwoche hinzufuegen
    KW = isoweek(Datetime),
    KW = factor(KW),
    # monat und Jahr
    Monat = month(Datetime),
    Monat = factor(Monat),
    Jahr = year(Datetime),
    Jahr = factor(Jahr))

Dies machen wir auch mit dem “meteo” Datensatz.

  • Wiederum bitte Wochentag, Werktag oder Wochenende, Kalenderwoche, Monat und Jahr. Ebenfalls alles als Faktor speichern.
Code 
# Wir gruppieren die Meteodaten noch nach Kalenderwoche und Werktag / Wochenende
# Dafür brauchen wir zuerst diese als Convenience Variablen
meteo <- meteo |>
  # wday sortiert die Wochentage automatisch in der richtigen Reihenfolge
  mutate(
    Wochentag = wday(time, week_start = 1),
    Wochentag = factor(Wochentag),
    # Werktag oder Wochenende hinzufuegen
    Wochenende = ifelse(Wochentag %in% c(6, 7), "Wochenende", "Werktag"),
    Wochenende = as.factor(Wochenende),
    # Kalenderwoche hinzufuegen
    KW = isoweek(time),
    KW = factor(KW),
    # monat und Jahr
    Monat = month(time),
    Monat = factor(Monat),
    Jahr = year(time),
    Jahr = factor(Jahr))
  • Später werden wir nicht nur Analysen pro Tag machen, sondern auch zusammengefasst nach Woche. Dafür müssen wir nun den meteo-Datensaz gruppieren und den mean berechnen. Hier der Code dazu, wie das aussehen könnte:
Code 
meteo_day <- meteo |>
  group_by(Jahr, Monat, KW, Wochenende) |>
  summarise(
    tre200nx = mean(tre200nx),
    tre200jx = mean(tre200jx),
    rre150n0 = mean(rre150n0),
    rre150j0 = mean(rre150j0),
    sremaxdv= mean(sremaxdv))

Wieder zurück zum Depo-Datensazt.

Ich mache den folgenden Punkt nachgelagert, da zu viele Operationen in einem Schritt auch schon mal etwas durcheinander erzeugen können.

Phase Covid (Code untenstehend). Wir definieren 5 Phasen:

  • von Anfang Untersuchungsperiode bis 1 Jahr vor Lockdown 1

  • Lockdown 1

  • zwischen den Lockdowns

  • Lockdown 2

  • Ende 2. Lockdown bis Ende Untersuchungsperiode

  • Wir packen alle Phasen (normal, die beiden Lockdowns und Covid aber ohne Lockdown) in eine Spalte –> long-format ist schöner (und praktischer für das plotten) als wide-format.

  • Später im multivariaten Modell werden die Levels der Variablen (z.B. bei der Phase Covid: Pre, Normal, Lockdown 1 und 2, Covid) per “default” alphabetisch geordnet und die Effektstärken der einzelnen Levels gegenüber dem ersten Level gerechnet. Das macht wenig Sinn, den die Levels sind nicht alphabetisch, sondern gemäss der Liste oben. Das passen wir ebenfalls an.

  • Hier der Code dazu:

Code 
depo <- depo |>
    mutate(Phase = case_when(
        Datetime < lock_1_start ~ "Pre",
        Datetime >= lock_1_start & Datetime <= lock_1_end ~ "Lockdown_1",
        Datetime > lock_1_end & Datetime < lock_2_start ~ "inter",
        Datetime >= lock_2_start & Datetime <= lock_2_end ~ "Lockdown_2",
        Datetime > lock_2_end ~ "Post"
    ))

# hat das gepklappt?!
unique(depo$Phase)

depo <- depo |>
    # mit factor() koennen die levels direkt einfach selbst definiert werden.
    # wichtig: speizfizieren, dass aus R base, ansonsten kommt es zu einem
    # mix-up mit anderen packages
    mutate(Phase = base::factor(Phase, levels = c("Pre", "Lockdown_1", "Inter", "Lockdown_2", "Post")))

str(depo)

Neben dem Lockdown können auch die Schulferien einen Einfluss auf die Besuchszahlen haben. Wir haben die Schulferien bereits als .csv eingelesen. Allerdings können wir die Schulferien nicht mit der case_when()-Funktion zuweisen, da diese mit dieser Funktion alle Vektoren im Datensatz “schulferien” verglichen werden, und nicht elementweise für jede Zeile im “depo”-Datensatz. Dies führt dazu, dass die Bedingungen nur einmal überprüft werden und dann auf den gesamten Vektor angewendet werden, anstatt Zeile für Zeile.

  • Weil dies etwas kompliziert ist, hier eine Funktion zur Zuweisung der Ferien, welche ihr kopieren könnt:
Code 
# schreibe nun eine Funktion zur zuweisung Ferien. WENN groesser als start UND kleiner als
# ende, DANN schreibe ein 1
for (i in 1:nrow(schulferien)) {
  depo$Ferien[depo$Datum >= schulferien[i, "Start"] & depo$Datum <= schulferien[i, "Ende"]] <- 1
}
depo$Ferien[is.na(depo$Ferien)] <- 0

# als faktor speichern
depo$Ferien <- factor(depo$Ferien)

3b)

  • Nun soll noch die volle Stunde als Integer im Datensatz stehen. Macht das mit dem Befehl hour()
Code 
# Fuer einige Auswertungen muss auf die Stunden als nummerischer Wert zurueckgegriffen werden
depo$Stunde <- hour(depo$Datetime)
# hour gibt uns den integer
typeof(depo$Stunde)

3c)

Die Daten wurden durch den WPZ kalibriert (Kommastellen).

  • Rundet sie auf 0 Nachkommastellen (Ganzzahl; unser Modell kann nicht mit Kommazahlen in der ahbängigen Variable umgehen). Der Befeht lautet round()

  • Definiert sie sicherheitshalber als Integer

  • Macht das für IN, OUT und Total.

Code 
depo$Total <- round(depo$Total, digits = 0)
depo$Fuss_IN <- round(depo$Fuss_IN, digits = 0)
depo$Fuss_OUT <- round(depo$Fuss_OUT, digits = 0)

3d) Tageszeit

Wir setzen den Fokus unserer Untersuchung auf die Veränderung der Besuchszahlen in der Abend- und Morgendämmerung sowie der Nacht. Dafür müssen wir diese tageszeitliche Einteilung der Daten erst machen. Da dies über den Umfang dieser Fallstudie geht, liefere ich euch hier den Code dazu.

Die wichtigsten Punkte:

  • Die Tageslänge wurde für den Standort Zürich (Zeitzone CET) mit dem Package “suncalc” berechnet. Dabei wurden Sommer- und Winterzeit berücksichtigt.

  • Die Einteilung der Tageszeit beruht auf dem Start und dem Ende der astronomischen Dämmerung sowie der Golden Hour. Der Morgen und der Abend wurden nach dieser Definition berechnet und um je eine Stunde Richtung Tag verlängert.

  • Untenstehenden Code könnt ihr einfach kopieren.

Hinweis: damit case_when() funktioniert, müsst ihr dplyr Version als 1.1.1 oder neuer haben.

Code 
# Einteilung Standort Zuerich
Latitude <- 47.38598
Longitude <- 8.50806

# Start und das Ende der Sommerzeit:
# https://www.schulferien.org/schweiz/zeit/zeitumstellung/


# Welche Zeitzone haben wir eigentlich?
# Switzerland uses Central European Time (CET) during the winter as standard time,
# which is one hour ahead of Coordinated Universal Time (UTC+01:00), and
# Central European Summer Time (CEST) during the summer as daylight saving time,
# which is two hours ahead of Coordinated Universal Time (UTC+02:00).
# https://en.wikipedia.org/wiki/Time_in_Switzerland

# Was sind Astronomische Dämmerung und Golden Hour ueberhaupt?
# https://sunrisesunset.de/sonne/schweiz/zurich-kreis-1-city/
# https://www.rdocumentation.org/packages/suncalc/versions/0.5.0/topics/getSunlightTimes

# Wir arbeiten mit folgenden Variablen:
# "nightEnd" : night ends (morning astronomical twilight starts)
# "goldenHourEnd" : morning golden hour (soft light, best time for photography) ends
# "goldenHour" : evening golden hour starts
# "night" : night starts (dark enough for astronomical observations)

lumidata <-
    getSunlightTimes(
        date = seq.Date(depo_start, depo_end, by = 1),
        keep = c("nightEnd", "goldenHourEnd", "goldenHour", "night"),
        lat = Latitude,
        lon = Longitude,
        tz = "CET"
    ) |>
    as_tibble()

# jetzt haben wir alle noetigen Angaben zu Sonnenaufgang, Tageslaenge usw.
# diese Angaben koennen wir nun mit unseren Zaehldaten verbinden:
depo <- depo |>
    left_join(lumidata, by = c(Datum = "date"))

depo <- depo |>
    mutate(Tageszeit = case_when(
        Datetime >= nightEnd & Datetime <= goldenHourEnd ~ "Morgen",
        Datetime > goldenHourEnd & Datetime < goldenHour ~ "Tag",
        Datetime >= goldenHour & Datetime <= night ~ "Abend",
        .default = "Nacht"
    )) |>
    mutate(Tageszeit = factor(Tageszeit, levels = c("Morgen", "Tag", "Abend", "Nacht"), ordered = TRUE))

# behalte die relevanten Var
depo <- depo |> dplyr::select(-nightEnd, -goldenHourEnd, -goldenHour, -night, -lat, -lon)

# Plotte zum pruefn ob das funktioniert hat
ggplot(depo, aes(y = Datetime, color = Tageszeit, x = Stunde)) +
    geom_jitter() +
    scale_color_manual(values = mycolors)

sum(is.na(depo))

# bei mir hat der Zusatz der Tageszeit noch zu einigen NA-Wertren gefueht.
# Diese loesche ich einfach:
depo <- na.omit(depo)
# hat das funktioniert?
sum(is.na(depo))

Aufgabe 4: Aggregierung der Stundendaten

4a)

Unsere Daten liegen im Stundenformat vor. Für einige Auswertungen müssen wir aber auf ganze Tage zurückgreifen können.

  • Die Stundendaten müssen zu ganzen Tagen aggregiert werden. Macht das wiederum einer Pipe. Bezieht folgende Gruppierungen (group_by()) mit ein: Datum, Wochentag, Wochenende, KW, Monat, Jahr, Phase. Summiert die Zählmengen separat (Total, IN, OUT) auf und speichert das Resultat unter depo_d.

Tipp: Wenn man die Convinience Variablen als grouping variable einspeisst, dann werden sie in das neue df übernommen und müssen nicht nochmals hinzugefügt werden

Code 
depo_d <- depo |> 
  group_by(VARIABLE1, VARIABLE2, ...) |>   # Gruppieren nach den Variablen
  summarise(Total = sum(Fuss_IN + Fuss_OUT),# Berechnen der gewünschten Werte
            Fuss_IN = sum(Fuss_IN),
            ...
Code 
# hier werden also pro Nutzergruppe und Richtung die Stundenwerte pro Tag aufsummiert
depo_d <- depo |>
    group_by(Datum, Wochentag, Wochenende, KW, Monat, Jahr, Phase) |>
    summarise(
        Total = sum(Fuss_IN + Fuss_OUT),
        Fuss_IN = sum(Fuss_IN),
        Fuss_OUT = sum(Fuss_OUT)
    )
# Wenn man die Convinience Variablen als grouping variable einspeisst, dann werden sie in
# das neue df uebernommen und muessen nicht nochmals hinzugefuegt werden
# pruefe das df
head(depo_d)
  • Erstellt nun einen Datensatz depo_daytime, in welchem ihr gruppiet nach:
  1. Jahr
  2. Monat
  3. Kalenderwoche
  4. Phase
  5. Ferien
  6. Wochenende oder Werktag
  7. Tageszeit
Code 
# nun gruppieren wir nicht nur nach Tag sondern auch noch nach Tageszeit
depo_daytime <- depo |>
  group_by(Jahr, Monat, KW, Phase, Ferien, Wochenende, Tageszeit) |>
  summarise(
    Total = sum(Fuss_IN + Fuss_OUT),
    Fuss_IN = sum(Fuss_IN),
    Fuss_OUT = sum(Fuss_OUT))
  • Weiter benötigen wir für für die Berechnung der Verteilung der Besuchenden über den Tag die durchschnittliche Besucheranzahl pro Stunde (mean), unterteilt nach Tageszeit und Phase (Gruppierungen Tageszeit, Phase). Speichert das unter “mean_phase_d”.
Code 
mean_phase_d <- depo_daytime |>
  group_by(Phase, Tageszeit) |>
  summarise(
    Total = mean(Total),
    IN = mean(Fuss_IN),
    OUT = mean(Fuss_OUT))

4b)

  • Aggregiere die Stundenwerte nach dem Monat (Gruppierungen Monat, Jahr) und speichert das neue df unter depo_m.

Tipp: Braucht wiederum group_by() und summarise(). Nun brauchen wir nur noch das Total, keine Richtungstrennung mehr.

Code 
depo_m <- depo |>
    group_by(Jahr, Monat) |>
    summarise(Total = sum(Total))
  • Fügt den neu erstellten df eine Spalte mit Jahr + Monat hinzu. Hier der fertige Code dazu:
Code 
depo_m <- depo_m |>
    mutate(
        Ym = paste(Jahr, Monat), # und mache eine neue Spalte, in der Jahr und
        Ym = lubridate::ym(Ym)
    ) # formatiere als Datum
  • Wiederholt diesen Schritt, diesmal aber mit der Gruppierung “Tageszeit” neben “Jahr” und “Monat” und speichert das Resultat unter “depo_m_daytime
Code 
# Gruppiere die Werte nach Monat und TAGESZEIT
depo_m_daytime <- depo |>
    group_by(Jahr, Monat, Tageszeit) |>
    summarise(Total = sum(Total))
# sortiere das df aufsteigend (nur das es sicher stimmt)

depo_m_daytime <- depo_m_daytime |>
    mutate(
        Ym = paste(Jahr, Monat), # und mache eine neue Spalte, in der Jahr und
        Ym = lubridate::ym(Ym)
    ) # formatiere als Datum

4c)

Macht euch mit den Daten vertraut. Plottet sie, seht euch die df’s an, versteht, was sie representieren.

Z.B. sind folgende Befehle und Plots wichtig:

  • str()

  • summarize()

  • head()

  • Scatterplot, x = Datum, y = Anzahl pro Zeiteinheit

  • Histrogram

  • usw.

Hinweis: Geht noch nicht zu weit mit euren Plots. Die Idee ist, dass man sich einen Überblick über die Daten verschafft und noch keine “analysierenden” Plots erstellt.

Nachdem nun alle Daten vorbereitet sind folgt im nächsten Schritt die deskriptive Analyse.