The R package dartle is a toolkit of functions for air quality model benchmarking, inspired by the DELTA tool (JRC-IES) and the work of the FAIRMODE WP1.
For dartle to work, you need R
(version 2.10 or higher) and the packages stats
, dplyr
, RcppRoll
, data.table
and ggplot2
. To install dartle downloading it from GitHub, use the package devtools
, as follows:
require("devtools") devtools::install_github("jobonaf/dartle")
Once installed, you can load it
library("dartle")
The package includes some datasets:
|Item |Title | |:-------------------|:-----------------------------| |mod.data (mod_data) |Forecasted concentrations | |obs.no2 (obs_data) |Observed NO2 concentrations | |obs.o3 (obs_data) |Observed ozone concentrations | |obs.pm10 (obs_data) |Observed PM10 concentrations | |obs.pm25 (obs_data) |Observed PM2.5 concentrations |
Let's have a look to the forecasted data (ignore EmissionTime
, Intercept
and Slope
)
library(knitr) kable(head(mod.data))
Let's see where and when have been observed the highest concentrations. Of PM10
require(dplyr) kable(head(obs.pm10 %>% arrange(desc(Value))))
PM2.5
kable(head(obs.pm25 %>% arrange(desc(Value))))
NO2
kable(head(obs.no2 %>% arrange(desc(Value))))
and ozone
kable(head(obs.o3 %>% arrange(desc(Value))))
Note that for PM10, NO2 and ozone only background stations are provided, while obs.pm25
includes also traffic and industrial stations. All the data are from the Environmental Agency of Friuli Venezia Giulia region (Italy).
unique(obs.pm25$StationType)
To produce a target plot, we must arrange forecasted and observed data in the same data.frame. First, we extract PM10 forecasts with filter
and calculate the daily averages with dMean
Mod <- dMean(mod.data %>% filter(Var=="c_PM10"), value = "Value", time = "Time", point = "Point")
Then, we adjust the format of the observed dataset
Obs <- obs.pm10 %>% mutate(Day=format(Time,"%Y-%m-%d"), Point=ID)
So, observed and forecasted datasets can be meld in a single data.frame (Dat.pm10
)
Dat.pm10 <- inner_join(Mod, Obs, by=c("Point", "Day"), suffix = c(".mod", ".obs"))
Now Dat
can be passed to the target_report
function, to calculate some quality indicators
t_rep <- target_report(Dat.pm10, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "PM10")
Finally, t_rep
(the output of the target_report
function) is ready to be passed to the target_plot
function
target_plot(t_rep)
The same processing for PM2.5:
Mod <- dMean(mod.data %>% filter(Var=="c_PM25"), value = "Value", time = "Time", point = "Point") Obs <- obs.pm25 %>% mutate(Day=format(Time,"%Y-%m-%d"), Point=ID) Dat.pm25 <- inner_join(Mod, Obs, by=c("Point", "Day"), suffix = c(".mod", ".obs")) t_rep <- target_report(Dat.pm25, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "PM2.5")
And the plot
target_plot(t_rep)
For NO2, the target plot is based on hourly data, therefore no daily average is performed
Mod <- mod.data %>% filter(Var=="c_NO2") %>% mutate() Obs <- obs.no2 %>% rename(Point=ID) Dat.no2 <- inner_join(Mod, Obs, by=c("Point", "Time"), suffix = c(".mod", ".obs")) t_rep <- target_report(Dat.no2, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "NO2")
Note that the function target_report
needs the argument pollutant
because the parameters used to assess the measurements uncertainties are pollutant-dependent.
target_plot(t_rep)
For ozone, we have to calculate daily maxima of the 8-hours running mean.
Mod <- dMaxAvg8h(mod.data %>% filter(Var=="c_O3"), value = "Value", time = "Time", point = "Point") Obs <- obs.o3 %>% rename(Point=ID)
In this case, both observations and forecasts are provided as hourly averages in the datasets, so we need the function dMaxAvg8h
again.
Obs <- dMaxAvg8h(Obs, value = "Value", time = "Time", point = "Point") Dat.o3 <- inner_join(Mod, Obs, by=c("Point", "Day"), suffix = c(".mod", ".obs")) t_rep <- target_report(Dat.o3, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "O3")
target_plot(t_rep)
s_rep <- summary_report(Dat.pm10, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "PM10") summary_plot(s_rep)
s_rep <- summary_report(Dat.pm25, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "PM2.5") summary_plot(s_rep)
s_rep <- summary_report(Dat.no2, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "NO2") summary_plot(s_rep)
s_rep <- summary_report(Dat.o3, obs = "Value.obs", mod = "Value.mod", point = "Point", pollutant = "O3") summary_plot(s_rep)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.