In Ostluft/rOstluft: Tools for handling air quality data by Ostluft

set.seed(1014)
options(digits = 3)

knitr::opts_chunk$set(
  collapse = TRUE,
  cache = FALSE,
  comment = "#>"
)

library(dplyr)
library(magrittr)

Einleitung

Dieses Dokument soll einem Anwender des Packages rOstluft::rOstluft einen Überblick der enthaltenen Funktionalität bieten und die gängigsten Arbeitsabläufe aufzeigen.

Der Fokus dieses Packages liegt in der Bereitstellung von Daten aus verschiedenen Datenquellen in einem einheitlichen Format zur Analyse der Luftqualität. Ausserdem enthält es Werkzeuge für einige übliche Aufgaben, die während solchen Analysen anfallen. Diese sind Umrechnungen zwischen verschiedenen Mittelungsintervallen, Statistische Methoden mit Berücksichtung der Datenverfügbarkeit, Umwandlung von Volumen- und Massenkonzentrationen und lesen von Daten vorliegend in verschiedenen Formaten.

Installation

Der Quellcode von rOstluft ist auf github gehosted. Die einfachste Variante ist die Installation mit Hilfe des Packages devtools:

#install.packages("devtools")
devtools::install_github("Ostluft/rOstluft")

Zusätzlich muss das Package aws.s3 manuell aus dem cloudyr Repositorium installiert werden, weil die CRANR Version veraltet ist:

install.packages("aws.s3", repos = c("cloudyr" = "http://cloudyr.github.io/drat"))

Ist dies wegen Einschränkungen durch Firewalls oder Proxies nicht möglich. Muss der Quellcode manuell von github heruntergeladen werden (Clone or download > Download as ZIP), entpackt und manuell installiert werden. Allerdings bestehen Abhängigkeiten zu Packages die auf CRAN bereitgestellt werden. Können auch keine CRAN Packages installiert werden, müssen zuerst alle CRAN Abhängkigkeiten und deren Abhängigkeiten installiert werden.

Zusätzlich besteht noch die Github Abhängkigkeit zu rOstluft.data. Dieses Packages muss auf die gleiche Weise zuerst installiert werden mit folgenden Schritten:

download.file("https://github.com/Ostluft/rOstluft/archive/master.zip", "rOstluft-master.zip")
download.file("https://github.com/Ostluft/rOstluft.data/archive/master.zip", "rOstluft.data-master.zip")

install.packages("devtools")
install.packages("aws.s3", repos = c("cloudyr" = "http://cloudyr.github.io/drat"))

deps <- c('dplyr', 'tidyr', 'lubridate', 'R6', 'rappdirs', 'tibble', 'base64url', 'forcats',
          'fs', 'purrr', 'readr', 'stringr', 'stringi', 'sp', 'rgdal', 'rlang', 'magrittr')

for (p in deps) {
  install.packages(p)
}

devtools::install_local("rOstluft.data-master.zip", dependencies = FALSE)
devtools::install_local("rOstluft-master.zip", dependencies = FALSE)

Falls das installieren von rOstluft scheitert, fehlt vermutlich eine Abhängigkeit. Welche das ist, kann der Fehlermeldung entnommen werden.

Nach der Installation kann das Packages verwendet werden:

library(rOstluft)

Abfrage von Daten

Ostluft Amazon AWS S3

Die zentrale Datenablage innerhalb von Ostluft erfolgt auf Amazon AWS S3.Aus Lizenztechnischen Gründen kann das Bucket nicht öffentlich zugänglich gemacht werden. Die Zugangsdaten werden von Jörg Sintermann vergeben. Die Zugangsdaten werden am einfachsten über eine .Renvirion Datei im Verzeichnis des RStudio Projektes oder im HOME Verzeichnis des Users[^3] dem Package zugänglich gemacht. Inhalt der .Renvirion Datei:

AWS_ACCESS_KEY_ID = "XXXXXXXXXXXXXXXXXXXX"
AWS_SECRET_ACCESS_KEY = "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
AWS_DEFAULT_REGION = "eu-central-1"

Weitere Möglichkeiten sind in der Dokumentation von aws.signature zu finden.

Sämtliche Daten die einmal von Amazon S3 geöffnet wurden, werden lokal auf dem Rechner gespeichert. Bei jedem folgenden Zugriff wird nur überprüft, ob die Daten noch identisch sind.

Der Zugriff auf die Daten erfügt dann über ein store Objekt, welches auf folgende Art initialisiert wird:

store <- storage_s3_rds("tutorial", format = format_rolf(), bucket = "rostluft", prefix = "aqmet")

# oder nutze den vordefinierten store (name ist optional, default ist "aqmet")
store <- store_aqmet("tutorial")

Dieses store Objekt verfügt über verschiedene Methoden, zur Abfrage von Daten sind jedoch nur zwei von Bedeutung:

• $get_content(): holt eine Übersicht über alle im store enthaltene Daten
• $get(): Aktuelle Datenabfrage. Welche Daten geholt werden, muss über Funktionsargumente definiert werden

content <- store$get_content()
content

dplyr::sample_n(content, 10)

Eine Datenabfrage holt immer die Daten von einem komplettem Jahr einer Station für einen bestimmten Mittelungszeitraum:

store$get(site = "Zch_Schimmelstrasse", year = 2014, interval = "min30")

Jedes der Argumente kann auch ein Vector sein mit mehreren Werte. Es wird dann Kombination aller Möglichkeiten abgefragt.

sites <- c("Zch_Schimmelstrasse", "Zch_Rosengartenstrasse", "Zch_Stampfenbachstrasse")
data <- store$get(site = sites, year = 2014:2015, interval = "min30")
data

dplyr::sample_n(data, 10)

Die letzte Funktionalität von $get() ist die Filterung der Daten vor der Rückgabe. Mit dem Argument "filter" als dplyr::filter() kompatibler Ausdruck:

store$get(site = sites, year = 2014, interval = "min30", filter = parameter == "PM10")

Lokales Arbeiten mit den AWS S3 Daten

Ein Nachteil des S3 Storage ist, dass bei jeder Datenabfrage eine Internetverbindung zur Überprüfung, ob aktualisierte Daten verfügbar sind, vorhanden sein muss. Selbst wenn die Daten bereits lokal gespeichert sind. Nicht nur wird eine Internetverbindung benötigt, die Überprüfung braucht auch eine kurze Zeit. Hat man alle notwendigen Daten bereits herunter geladen, kann man jedoch mit einem lokalen Store arbeiten. Praktischerweise hat der S3 Store eine Funktion, welche einen lokalen Store zurückgibt.

lokal <- store$get_local_storage()

Dieser verfügt über die gleichen Funktionalität wie der S3 Store:

lokal$get_content()

lokal$get(site = sites, year = 2014, interval = "min30", filter = parameter == "PM10")

Dem aufmerksamen Leser ist vermutlich aufgefallen, dass die Inhaltsübersicht von $get_content() nicht dem lokalen Inhalt entspricht. Es ist immer noch die Übersicht welche Daten in S3 verfügbar wären. Will man wissen welche Daten lokal verfügbar sind, hilft einem die Funktion $list_chunks() weiter.

tibble::glimpse(lokal$list_chunks())

Auch der S3 Store verfügt über die $list_chunks() Funktion.

tibble::glimpse(store$list_chunks())

Im normal Fall ist Nutzung von $get_content() über $list_chunks() zu bevorzugen. Die Rückgabe von $get_content() enthält die gemessen Parameter und die Anzahl der gültigen Punkte. Bei $list_chunks() hingegen sind einige interne Informationen über den Store enthalten.

Um die Vorbereitung um mit einem lokalen Store zu arbeiten zu vereinfachen verfügt der S3 Store über die Funktion $download(...). Die dot Argumente werden als Filter auf die Rückgabe von $list_chunks() angewendet. Ohne Argumente wird der komplette Store heruntergeladen. Folgendes Beispiel lädt samtliche 30 Minutenmittelwerte für die Station Rosengartenstrasse nach dem Jahr 2015 (Achtung == verwenden, es sind Filter Ausdrücke!):

store$download(site == "Zch_Rosengartenstrasse", interval == "min30", year > 2015)

lokal$list_chunks() %>% dplyr::select("site", "year", "interval")

Eigene Daten in einem lokalen Store

Der aqmet S3 Store enthält nur bereinigte Daten von abgeschlossenen Jahren. Werden Daten benötigt, die nicht in S3 vorhanden sind, bietet sich die Nutzung eines seperaten lokalen Stores an:

my_store = storage_local_rds("eigene_daten", format_rolf(), read.only = FALSE)

In diesen kann man nun Daten mit der Store Funktion $put() oder der Hilfsfunktion import_directory() importieren:

examples_path <- system.file("extdata", package = "rOstluft.data")
import_directory(my_store, examples_path, read_airmo_csv, glob = "*Jan.csv")

fn <- fs::path(examples_path, "Zch_Rosengartenstrasse_2010-2014.csv")
data <- read_airmo_csv(fn)
my_store$put(data)

Wie gewohnt kann erhält man mit $get_content() eine Übersicht und die Daten mit $get():

my_store$get_content()

my_store$get(site = "Zch_Stampfenbachstrasse", interval = c("d1", "h1"), year = 2013)  %>%
  dplyr::arrange(.data$starttime, .data$parameter)

Einlesefunktionen von Daten

Um Daten in R einzulesen exisitieren für folgende Quellen bereits Funktionen:

• Ostluft: read_airmo_csv() und read_airmo_dat() für Exporte aus der AIRMO
• MeteoSchweiz: read_smn() und read_smn_multiple() für SwissMetNet Exporte
• ETHZ/IAC: read_ethz_iac() Meteo Stationen auf dem Hönggerberg und dem Hauptgebäude CHN

Existiert für die vorliegenden Daten keine Funktion und man möchte die selber schreiben an dieser Stelle einige Hinweise und Tipps.

Verwende readr Funktionen und definiere die Zeichenkodierung

Ein generelles Problem beim einlesen von Textdateien ist die Zeichenkodierung. Die meisten haben vermutlich "�" schon in Texten gesehen. Alle Daten im aqmet Store und alle Rückgabewerte von read_xxx() Funktionen sind UTF-8 kodiert. Mit dem readr Packages wird bei korrekter Definition der Locale der Text automatisch zu UTF-8 konvertiert. Im Gegensatz zu den R Base Funktionen oderdata.table::fread().

Ausgabeformat rolf

Für ein reibungsloses Zusammenspiel mit den Stores zu garantieren muss das Format rolf eingehalten werden. rolf setzt einen tibble mit folgenden Spalten und Klassen voraus:

| Spalte | Klasse | |-----------------|---------| | "starttime" | POSIXct | | "site" | factor | | "parameter" | factor | | "interval" | factor | | "unit" | factor | | "value" | double |

Selbst wenn in der Datei keine Einheiten enthalten sind, ist es besser die Spalte unit mit NA zu initialisieren. Ein weiterer Trick ist es die Spalten am Ende explizit zu selektionieren. Hier ein Beispiel, wie es in read_smn() gelöst ist:

data <- dplyr::mutate(data,
  stn = forcats::as_factor(.data$stn),
  parameter = forcats::as_factor(.data$parameter),
  interval = forcats::as_factor(interval),
  unit = factor(NA)
)
dplyr::select(data, starttime = "time", site = "stn", "parameter", "interval", "unit", "value")

Mit dem dplyr::select am Ende ist man sicher nur die Spalten in der Rückgabe zu haben die man haben will und es ist einfach noch Spalten umzubenennen und die Reihenfolge zu ändern.

Zusammenfügen von Daten im rolf Format

Die Verwendung von Factoren im rolf Format hat seine Vor- und Nachteile. Der Hauptvorteil liegt bei schnelleren Lese- und Schreiboperation, der grösste Nachteil beim Zusammenfügen der Daten. Jedoch unterstützt seit dplyr v1.0.0 dplyr::bind_rows das Zusammenfügen von Factor Spalten. Diese Funktion unterstützt ausserdem Quasiquotation von rlang[^1]. Somit kann auch eine grosse Liste von tibbles im rolf Format komfortable aneinander gefügt werden:

df1 <- read_airmo_csv(fs::path(examples_path, "Zch_Stampfenbachstrasse_h1_2013_Jan.csv"))
df2 <- read_airmo_csv(fs::path(examples_path, "Zch_Stampfenbachstrasse_d1_2013_Jan.csv"))
df_comb <- dplyr::bind_rows(df1, df2)
df_comb %>% dplyr::arrange(.data$starttime, .data$parameter)

# liste mit 20+ tibbles
df_list <- read_smn_multiple(fs::path(examples_path, "smn_multi.txt"), as_list = TRUE)
length(df_list)
df_comb <- dplyr::bind_rows(!!!df_list)
df_comb %>% dplyr::arrange(.data$starttime, .data$parameter)

dplyr::bind_rows() entfernt jedoch keine Duplikate aus den Daten. Es hängt die Daten einfach aneinander. Möchte man Duplikate entfernen, bzw. Daten aktualisieren mit neueren Daten muss man die Funktion $merge() des Format Objekts benützen. Dieses wurde beim Erzeugen des Store direkt initialisiert und ist als Feld $format verfügbar. Alternativ kann ein seperates Format Objekt erzeugt werden. Die Daten im ersten Datenframe werden über Daten im zweiten Datenframe priorisiert.

# kopiere Daten und setze einige NA -> Resultat sollte NA enthalten
df1 <- df2
df1$value[1:150] <- NA

df_comb <- store$format$merge(df1, df2)
df_comb %>% dplyr::arrange(.data$starttime, .data$parameter)

rolf <- format_rolf()
df_comb <- rolf$merge(df1, df2)
df_comb %>% dplyr::arrange(.data$starttime, .data$parameter)

Die $merge() Funktion unterstützt nur 2 Argumente. Muss eine Liste von tibbles gemerget werden muss man purrr::reduce() benützen. Mit Hilfe des Argument .dir kann man die Priosierung setzen. Mit "forward"" haben die Daten die zuerst in der Liste auftauchen Priorität, mit "backward" die Letzten.

df_comb <- purrr::reduce(df_list, rolf$merge, .dir = "forward")
df_comb %>% dplyr::arrange(.data$starttime, .data$parameter)

Metadaten

Daten aus verschiedenen Quellen werden verschiedene Schemen für die Vergabe von Namen haben. MeteoSchweiz zum Beispiel kodiert im Parameter Namen gleichzeitig das Mittelungsinterval:

| Parameter | Interval | Messgrösse | |-----------|---------------------------|------------------------------------------| | ta1towb0 | Bürgerliche Stundenmittel | Lufttemperatur, Instrument 1 | | ta1towh0 | Stundenmittel | Lufttemperatur, Instrument 1 | | ta1tows0 | Momentanwerte | Lufttemperatur, Instrument 1 | | fk1towz0 | Zehnminutenmittel | Windgeschwindigkeit skalar, Instrument 1 |

Um die Bezeichnungen zu normalisieren kommt es zu einem Zusammenspiel von im Store bereitgestellte Daten und der Funktion meta_apply(). Das Store Objekt verfügt über die Funktion $get_meta(), welche im AWS S3 bereitgestellte Metadaten zurück gibt. Wird $get_meta() ohne Argument aufgerufen erhält man eine Liste mit sämtlichen Metadaten im Store. Kennt man bereits den Namen kann man direkt die entsprechenden Metadaten holen:

meta <- store$get_meta()
names(meta)

meteoschweiz <- store$get_meta("meteoschweiz")
tibble::glimpse(meteoschweiz)

Die meta_apply() benützt die Metadaten als "Lookup Table". Es wird folgende Operation ausgeführt:

meta_apply(data, meta, data_src, data_dest, meta_key, meta_val)

data$data_dest = meta[meta$meta_key == data$data_src]$meta_val

Wobei die jeweiligen Spalten jeweils als Argument übergeben werden. Die Funktion verfügt über verschiende Modi um mit fehlenden Einträgen umzugehen. Im folgenden Beispiel werden bei den MeteoSchweiz Daten die "unit" Spalte basierend auf den "parameter" Spalte ausgefüllt. Da die Parameter noch die MeteoSchweiz Bezeichnungen haben müssen wir als Meta Key die Spalte "parameter_original" und "unit" als Value verwenden. Es werden die Daten aus rOstluft.data verwendet. In diesen fehlt in den Metadaten der MeteoSchweiz Parameter rre150z0 (Niederschlag). Der Default Modus "strict" stoppt das Script:

meteoschweiz <- readRDS(fs::path(examples_path, "meta_smn.rds"))
data <- read_smn(fs::path(examples_path, "smn.txt"), na.rm = FALSE)
data <- dplyr::arrange(data, .data$starttime, .data$parameter)
data

tryCatch({
  meta_apply(data, meteoschweiz, "parameter", "unit", "parameter_original", "unit")
}, error = function(e) {
  sprintf(e$message)
})

Im einfachsten Fall werden Zeilen mit fehlenden Einträge einfach gelöscht. Dies geschieht mit dem Modus "drop":

df <- meta_apply(data, meteoschweiz, "parameter", "unit", "parameter_original", "unit", mode = "drop")
df

Sollte der Parameter für spätere Bearbeitung behalten werden, kommt der Modus "keep" zum Zug:

df <- meta_apply(data, meteoschweiz, "parameter", "unit", "parameter_original", "unit", mode = "keep")
df

Mit Hilfe des Modus "replace" kann eine zusätzliche "Lookup Table" übergeben werden. Werte in dieser Tabelle haben Priorität über Werte in der Meta Tabelle. Im Beispiel wird die fehlende Einheit für rre150z0 bereit gestellt und die Einheit für dkl010z0 von ° zu deg überschrieben.

df <- meta_apply(data, meteoschweiz, "parameter", "unit", "parameter_original", "unit",
                 mode = "replace", replacements = list(rre150z0 = "mm", dkl010z0 = "deg"))
df

Hier ein Beispiel für die komplette Überführung der MeteoSchweiz Daten, inkl. manuelle Umbennung von rre150z0 und Änderung der Einheit für die Windrichtung (WD) von ° zu deg:

df <- meta_apply(data, meteoschweiz, "parameter", "unit", "parameter_original", "unit",
                 mode = "replace", replacements = list(rre150z0 = "mm", dkl010z0 = "deg"))

df <- meta_apply(df, meteoschweiz, "parameter", "parameter", "parameter_original", "parameter",
                 mode = "replace", replacements = list(rre150z0 = "Niederschlag"))

df <- meta_apply(df, meteoschweiz, "site", "site", "site_short", "site")

df

openair kompatibles Wide Format

Das Long Format ist nicht für jede Analyse geeignet. Aus diesem Grund enthält rOstluft die Funktionalität Daten im rolf Format in ein openair kompatibles Wide Format umzuwandeln. Ein Nachteil am Wide Format ist, dass Informationen die zu einer Spalte gehören, wie die Einheit einer Messgrösse, verloren gehen. Dies führt auch dazu, das man nicht die gleiche Messgrösse mit unterschiedlichen Einheiten in den Daten haben kann. Ausser man nennt die Messgrösse um. Die Funktion rolf_to_openair() geht folgende Kompromisse ein: Die Einheiten werden entfernt, aber im Attribute "units" des Datenframes gespeichert. Messgrössen mit der Einheit ppb oder ppm werden ebenfalls entfernt (Dieser Automatismus lässt sich mit dem Argument keep_ppb abschalten). Openair nimmt an, dass alle Messgrössen als Massenkonzentrationen vorliegen. Weitere Konventionen in openair sind das die Zeitinformation in der Spalte "date" als Date oder POSIXct vorliegen, die Windgeschwindigkeit in der Spalte "ws" und die Windrichtung in der Spalte "wd". Sämtliche Character und Factor Spalten dienen als Gruppierungsspalten. Diese werden von openair ignoriert oder gar entfernt, wenn nicht das Argument type benützt wird. Dieses erlaubt die explizite Nutzung einer Spalte als Gruppierungskriterium. Es können mehrere Gruppierungskriterien getrennt durch ein Komma definiert werden. Die openair Funktion openair::cutData() übernimmt die Gruppierung.

data <- read_airmo_csv(fs::path(examples_path, "Zch_Stampfenbachstrasse_2010-2014.csv"))
wide <- rolf_to_openair(data)
tibble::glimpse(wide)

Diese Daten können nun direkt in openair verwendet werden:

openair::windRose(wide)
openair::pollutionRose(wide, "NO2", type="year")

openair_to_rolf() konventiert Daten zurück ins rolf Format. Der Nutzer muss die fehlenden Informationen bereitstellen:

openair_to_rolf(wide, interval = "min30", units = attr(wide, "units"))

rolf_to_openair_single() bietet die Möglichkeit einen bestimmten Parameter rauszupicken:

rolf_to_openair_single(data, "NO2", unit = "µg/m3", keep_interval = TRUE)

Die weitere Optionen für die Funktion rolf_to_openair() sind in der Dokumentation zu finden.

Berechnung von Statistiken

Statistiken und die Verdichtung von Daten ist ein grundlegender Bestandteil einer Datenanalyse. In rOstluft gibt es die Funktion resample() zur Berechnung von Statistiken von einem Interval zu einem anderen. Die Funktion calculate_statstable() erlaubt die Definition von Berechnungen über mehrere verschiedene Intervale, ist aber weniger flexibel.

resample

resample() aggregiert die Daten in folgenden Schritten:

• Aufsplitten der Daten in Serien (site, parameter, interval, unit)
• pad_serie(): Auffüllen von Lücken mit NA (kann übersprungen werden)
• Gruppiere Serie mit neuem Interval mit Hilfe von lubridate::floor_date()
• Berechne definierte Statistiken für den Parameter
• kombiniere die neu berechneten Serien

data <- read_airmo_csv(fs::path(examples_path, "Zch_Stampfenbachstrasse_min30_2013_Jan.csv"))

# Behalte nur Massenkonzentrationen
data <- dplyr::filter(data, !(.data$unit == "ppb" | .data$unit == "ppm"))

statistics <- list(
  "default_statistic" = "drop",
  "O3" = list("mean", "perc95", "n", "min", "max"),
  "RainDur" = "sum",
  "WD" = "wind.direction",
  "WVs" = "wind.speed_scalar",
  "WVv" = "wind.speed_vector"
)

# Monatsmittelwerte
resample(data, statistics, "m1")

# Tagesmittelwerte
resample(data, statistics, "d1") %>%
  dplyr::arrange(.data$starttime)

resample() verfügt über die Funktionalität Datenverfügbarkeit und Lücken mit den Argumenten data_thresh und max_gap zu berücksichtigen:

# Jahresmittelwerte nur mit Daten vom Januar
resample(data, statistics, "y1")

# Jahresmittelwerte 80% Datenverfügbarkeit => all NA
resample(data, statistics, "y1", data_thresh = 0.8)

# Jahresmittelwerte 10 Tage Lücke = 48' 30min Werte => all NA
resample(data, statistics, "y1", max_gap = 480)

default_statistic

In der Liste werden für jeden Parameter die Statistiken definiert, welche angewendet werden sollen. Der Listeneintrag default_statistic wird auf alle nicht eingetragenen Parameter angewendet. Die Statistik drop schliesst einen Parameter aus den Resultaten aus. Dies erlaubt die schnelle Berechnung von Statistiken für ein paar wenige Parametern. Gleichzeitig ist es auch möglich die gleichen Statistiken auf viele Parameter anzuwenden und ein paar wenige Ausnahmen zu definieren:

statistics <- list(
  "default_statistic" = "mean",
  "RainDur" = "sum",
  "StrGlo" = "max",
  "WD" = "drop",
  "WVv" = "drop"
)

resample(data, statistics, "m1")

Automatismen und Einschränkungen

Die Funktion benennt für die meisten Statistiken den Parameter gemäss dem AIRMO Schema[^2] um. Die Statistiken mean, median, sum, sd, percentile benennen den Paremeter nicht um und sind somit nicht kombinierbar.

Das Auffüllen der Lücken kann mittels des Argument skip_padding übersprungen werden. Zusätzlich kann mittels start_date, end_date und drop_last das Füllen kontrolliert werden. Normalerweise nimmt die Funktion das letzte Datum und füllt bis zum letzten neuen Intervall auf. Werden zum Beispiel Tagesmittelwerte berechnet und es sind nur Daten bis zum 22. Dezember 2017 10:00 Uhr vorhanden, wird der letzte Wert in der Serie 2017-12-22 00:00 sein. Wird eine komplette Jahresreihe benötigt, kann entweder das end_date auf 2017-12-31 oder auf 2018-01-01 und drop_last auf TRUE gesetzt werden.

Windmittelung

Ein weiterer Spezialfall sind die Windberechnungen. Für die vektorielle Mittelung müssen die Paremeter zusammengefasst werden. Es müssen für alle drei Statistiken (wind.direction, wind.speed_vector, wind.speed_scalar) die Parameter definiert werden. Ist nur eine der Geschwindigkeiten vorhanden, wird die andere auf Basis der Vorhandenen berechnet. Zusätzlich ist die Windmittelung nicht kombinierbar mit anderen Statistiken.

Statstable

Werden viele verschiedene und mehrstufige Statistiken benötigt ist die Berechnung mit resample() umständlich. Die Funktion calculate_statstable() hingegen bietet diese Funktionalität, hat aber folgende Einschränkungen:

• Datenreihen werden immer auf komplette Jahre aufgefüllt
• max_gap wird nur auf Berechnungen zu Jahresintervallen berücksichtigt und in Tagen definiert
• data_thresh ist für alle Berechnungen identisch.
• Intervall h8gl erwartet Intervall h1 als Basis
• Die Verwendung von default_statistic und _input_ führt schnell zu unerwarteten Ergebnissen

Zur Analyse von Luftqualität sind diese Einschränkungen kein Problem, bzw. sogar erwünscht. Die Funktion calculate_statstable() erwartet wie der Name schon vermuten lässt die Defintion der Statistiken in einer Tabellenform. Die Tabelle hat die vier Spalten "parameter", "statistic", "from", "to". Jede Zeile enthält eine zu berechnende Statistik für einen Parameter vom Basis Interval ("from") zum neuen Interval ("to"). Die Tabelle wird dann gruppiert zuerst mit "from", dann mit "to" und zuletzt mit Parameter. Auf diese Weise wird eine Liste generiert die mit resample() kompatibel ist. Ist keine default_statistic definiert wird automatisch default_statistic = "drop" zur Liste hinzugefügt. In welcher Reihenfolge die verschiedenen "from" Intervalle berechnet werden, kann optional mit Hilfe des Argument order definiert werden. Der Default Wert ist c("input", "h1", "h8gl", "d1", "m1", "y1"). Um die Tabelle in einer kompakten Form zu erstellen kann in jeder Zelle mehrere Werte getrennt durch "," zusammengefasst werden. Hier die Tabelle zur Berechnung der LRV Grenzwerte aus 30min Werten (entspricht 6 Aufrufe von resample()):

lrv_table <- tibble::tribble(
  ~parameter, ~statistic, ~from, ~to,
  "SO2, NO2, PM10", "mean", "input", "y1",
  "SO2, NO2", "perc95", "input", "y1",
  "O3", "perc98", "input", "m1",
  "O3", "mean", "input", "h1",
  "O3", "n>120", "h1", "y1",
  "SO2, NO2, CO, PM10", "mean", "input", "d1",
  "SO2", "n>100", "d1", "y1",
  "NO2", "n>80", "d1", "y1",
  "CO", "n>8", "d1", "y1",
  "PM10", "n>50", "d1", "y1"
)

knitr::kable(lrv_table)

Wird diese Tabelle verwendet resultiert eine Liste mit h1, m1 und y1 Einträgen, welche die berechneten Grössen enthält:

data <- read_airmo_csv(fs::path(examples_path, "Zch_Stampfenbachstrasse_min30_2017.csv"))

# Berechne Massenkonzentrationen aus Volumenkonzentrationen
data <- calculate_mass_concentrations(data)

stats <- calculate_statstable(data, lrv_table)
stats

Nur die m1 und y1 sind für die LRV relevant:

lrv <- dplyr::bind_rows(stats$y1, stats$m1)

knitr::kable(lrv, digits = 0)

_inputs_

Der Parameter _inputs_ ist ein Hilfskonstrukt wie default_statistic. Ein Problem bei mehrstufigen Berechnungen ist, dass unter Umständen bereits Resultate im Input für weitere Berechnungen enthalten sind. In folgenden Bespiel wird zuerst aus den 30min Werten d1 Mittelwerte, Anzahl Datenpunkte, min und max berechnet für d1, m1 und y1. Danach sollte das Tagesmaximum für das Jahr bestimmt werden:

# Ein Parameter reicht für die Veranschaulichung
O3 <- dplyr::filter(data, .data$parameter == "O3")

statstable <- tibble::tribble(
  ~parameter, ~statistic, ~from, ~to,
  "default_statistic", "mean, n, min, max", "input", "d1, m1, y1",
  "default_statistic", "max", "d1", "y1"
)

# Keine Berücksichtung der Datenverfügbarkeit
stats <- calculate_statstable(O3, statstable, data_thresh = 0, max_gap = NULL)
stats$y1

Wird für die Berechnung des maximalen Tagesmittelwert des Jahres _inputs_ statt default_statistic verwendet, werden nur die Statistiken für Parameter berechnet die in den Input Daten enthalten sind.

statstable <- tibble::tribble(
  ~parameter, ~statistic, ~from, ~to,
  "default_statistic", "mean, n, min, max", "input", "d1, m1, y1",
  "_inputs_", "max", "d1", "y1"
)

# Keine Berücksichtung der Datenverfügbarkeit
stats <- calculate_statstable(O3, statstable, data_thresh = 0, max_gap = NULL)
stats$y1

Im Allgemeinen dürfte die explizite Definition der Statistiken für jeden Parameter gegenüber der Verwendung von default_statistic und _inputs_ vorzuziehen. Ansonsten kann das Resultat unerwartete Ergebnisse enthalten. Die Definition der statstable kann auch in einer Textdatei oder Exceldatei erfolgen und dann eingelesen werden, statt im Code definiert werden.

aqmet_stats <- statstable_default()
expanded_stats <- statstable_expand(aqmet_stats)
n_stats <- nrow(expanded_stats)
n_inputs <- nrow(dplyr::filter(expanded_stats, from == "input") %>% dplyr::distinct(parameter))

An dieser Stelle noch die Berechnungen aller Statistiken im aqmet Store. Aus r n_inputs Input Grössen werden r n_stats Statistiken berechnet:

knitr::kable(aqmet_stats)

Umwandlung von Volumen- und Massenkonzentrationen

Die Umwandlung der Konzentrationen erfolgt äquivalent zu den Berechnungen in der AIRMO. Sollen andere Konstanten verwendet werden, ist die Dokumentation von convert_set_R() zu konsultieren. Die Funktionen convert_conc() und convert_conc_multiple() verfügen beide über Varianten wie die umgewandelten Daten zurückgegeben werden. Definiert wird dies mit dem Argument method. Bei "return" werden nur die umgewandelten Daten zurückgegeben, bei "append" werden sie am Ende angehängt und zuletzt bei "replace" werden die ursprünglichen Daten ersetzt.

data <- read_airmo_csv(fs::path(examples_path, "Zch_Stampfenbachstrasse_min30_2013_Jan.csv"))

# Behalte nur Volumenkonzentrationen
data <- dplyr::filter(data, .data$unit == "ppb" | .data$unit == "ppm")

convert_conc(data, "NO", "ppb", "µg/m3", method = "return")

conversions <- tibble::tribble(
  ~parameter, ~from, ~to,
  "CO", "ppm", "mg/m3",
  "NO", "ppb", "µg/m3",
  "O3", "ppb", "µg/m3",
  "NO2", "ppb", "µg/m3",
  "SO2", "ppb", "µg/m3"
)
convert_conc_multiple(data, conversions, method = "return") %>%
  dplyr::arrange(.data$starttime)

convert_conc_multiple(data, conversions, method = "append") %>%
  dplyr::arrange(.data$starttime)

convert_conc_multiple(data, conversions, method = "replace") %>%
  dplyr::arrange(.data$starttime)

Die Funktion calculate_mass_concentrations() ist eine Hilfsfunktion, die aus Volumenkonzentrationen für CO, NO, NO2, NO2, SO2 in einem tibble automatisch die Massenkonzentrationen berechnet. Mit dem Argument keep_ppb können optional die Volumenkonzentration behalten werden. Das Default vorgehen ist die Volumenkonzentrationen zu ersetzen.

Pluck

Die Pluck Funktionen sind Wrapper um dplyr::filter(). Sie ermöglichen dem Nutzer eine komfortable Filterung von Daten im rolf Format. Als Beispiel nehmen wir die Daten aus dem mehrfach Meteoschweiz Export, dieser enthaltet eine Vielzahl von verschiedenen Stationen, Intervalle und Parameter:

fn <- system.file("extdata", "smn_multi.txt", package = "rOstluft.data")
data <- read_smn_multiple(fn) %>% 
  dplyr::arrange(starttime)

pluck_parameter(data, "ta1towb0")

# entpackt automatisch Vectors
pluck_site(data, c("KLO", "UEB")) %>% dplyr::slice(41:46)

# supported splicing mit den rlang operator !!!
intervals = c("h1", "d1")
pluck_interval(data, !!!intervals) %>% dplyr::slice(40:45)

# pluck_year mit einer sequence
pluck_year(data, 2010:2018) %>% dplyr::slice(51:56)

# NA Werte werden gefiltert
pluck_unit(data, "hPa")

# Verkettbar mittels pipes
data %>%
  pluck_site("KLO", "UEB") %>%
  pluck_parameter("gre000z0") %>%
  pluck_year(2010:2018) %>% 
  dplyr::slice(1:6)

Koordinatentransformationen

Für die Darstellung von Elementen mit Koordinaten im Bezugsystem LV95 oder LV03 auf Karten ist es oft notwendig diese nach WSG84 zu transformieren. Die Funktionen rOstluft::transform_LV95_to_WSG84() und rOstluft::transform_WSG84_to_LV95() sind Wrapper um die allgemeine Funktion rOstluft::transform_crs().

meteoschweiz <- readRDS(fs::path(examples_path, "meta_smn.rds"))
meteoschweiz <- dplyr::distinct(meteoschweiz, site, x, y)

transform_LV95_to_WSG84(meteoschweiz) %>% head

transform_crs(meteoschweiz, c(lon = "x", lat = "y"),
          sp::CRS("+init=epsg:2056"), sp::CRS("+init=epsg:4326"),
          append = FALSE) %>% head

NOAA Hysplit Trajektorien

Neben dem Ostluft AWS Store aqmet existiert noch ein Store mit hysplit Trajektorien von 1980 bis zum jeweiligen Vormonat des aktuellen Datums für die Standorte Zürich Kaserne und St. Gallen Blumenbergplatz. Die Trajektorien sind jeweils für 01:00 und 13:00 UTC+1 (=Luftmassen Herkunft) gerechnet.

Beim Format hysplit wird für die Abfrage nur die site und das year benötigt:

hysplit_store <- storage_s3_rds("hysplit_tutorial", format_hysplit(), "rostluft", "hysplit")

# nutze den vordefinierten store (name ist optional, default ist "hysplit")
hysplit_store <- store_hysplit("hysplit_tutorial")

# hole die Trakektorien für Zürich Kaserne von 2018
traj <- hysplit_store$get(site = "ZH_Kaserne_hysplit", year = 2018)

# plotte die Daten mit openair
openair::trajPlot(openair::selectByDate(traj, start = "2018-01-01", end = "2018-01-05"))

Lokale Daten löschen

Werden die Daten lokal nicht mehr benötigt, kann mit der Funktion $destroy() des Store Objektes die Daten gelöscht werden. Zu beachten ist, dass der AWS S3 Store als read only initialisiert wurde. Um ihn löschen zu können muss dies geändert werden. Um ein versehentliches löschen durch Autocomplete zu verhindern, muss der String "DELETE" als Argument übergeben werden

store$read.only = FALSE
store$destroy("DELETE")

my_store$destroy("DELETE")

hysplit_store$read.only = FALSE
hysplit_store$destroy("DELETE")

[^1]: rlang muss nicht als library importiert werden. ... leitet die Expression an rOstluft weiter. Dieses hat den Operator importiert und kann somit die Expression korrekt auswerten.

[^2]: Ausnahme für die Statistiken perc02, perc95, perc98. Die AIRMO spezifiert nicht das Basis Interval => O3_95%, In rOstluft O3_95%_min30.

[^3]: Der Pfad zum Verzeichnis entspricht der Ausgabe von Sys.getenv("HOME")

Ostluft/rOstluft documentation built on Feb. 6, 2024, 1:26 a.m.