In jobonaf/dartle: Air Quality Model Benchmarking

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")

Datasets

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)

Target plots

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)

Summary reports

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)

References

• Guidance Document on Modelling Quality Objectives and Benchmarking

---

Logo made by Freepik from www.flaticon.com is licensed by CC 3.0 BY

jobonaf/dartle documentation built on May 29, 2019, 4:52 p.m.