#' \code{target_report}: prepare input for target plot
#'
#' @title Target plot
#' @references Janssen et al., 2017. "Guidance Document on Modelling Quality Objectives and Benchmarking. Version 2.1"
#' @param data a data frame including observations and model forecasts
#' @param obs name of the column with observed data
#' @param mod name of the column with forecasts
#' @param point name of the column with station ID
#' @param pollutant one of \code{"NO2"}, \code{"O3"}, \code{"PM10"}, \code{"PM2.5"}
#' @param beta parameter \eqn{\beta} (default is 2)
#' @param t_rep output of \code{target_report}
#' @name target_plot
#'
#' @return \code{target_report} returns a list of 3:
#' \describe{
#' \item{\code{quality_points}}{a data frame with 9 variables for each station \describe{
#' \item{\code{Point}}{station ID}
#' \item{\code{Obs.ave}}{annual mean of observations (paired with forecasts)}
#' \item{\code{Mod.ave}}{annual mean of forecasts (paired with observations)}
#' \item{\code{rmsu}}{\eqn{RMS_U} (see \link{RMS_U_obs})}
#' \item{\code{crmse_norm}}{\eqn{CRMSE/(\beta*RMS_U)}{CRMSE/(\beta RMS_U)}}
#' \item{\code{r}}{R (see \link{R})}
#' \item{\code{crmse_ratio}}{\eqn{abs(NMSD)/sqrt(2*(1-R))}{\frac{|NMSD|}{\sqrt{2(1-R)}}} as in eq.34, p.30 in Janssen et al., 2017}
#' \item{\code{bias_norm}}{\eqn{BIAS/(\beta*RMS_U)}{BIAS/(\beta RMS_U)}}
#' \item{\code{mqi_ts}}{Model Quality Indicator \eqn{RMSE/(\beta*RMS_U)}{RMSE/(\beta RMS_U)} (see \link{MQI_ts_synth})}
#' \item{\code{mqi_year}}{Model Quality Indicator for yearly averages (see \link{MQI_year})}
#' \item{\code{n_valid}}{no. of valid data}
#' }}
#' \item{\code{quality_overall}}{a data frame with 3 overall indicators \describe{
#' \item{\code{mqi_ts_p90}}{90th percentile of \code{mqi_ts}}
#' \item{\code{mqi_year_p90}}{90th percentile of \code{mqi_year}}
#' \item{\code{n_points}}{no. of valid stations}
#' }}
#' \item{\code{parameters}}{a list of 7 parameters \describe{
#' \item{\code{U_RV95r}}{see \link{params_U}}
#' \item{\code{alpha}}{see \link{params_U}}
#' \item{\code{RV}}{see \link{params_U}}
#' \item{\code{Np}}{see \link{params_U}}
#' \item{\code{Nnp}}{see \link{params_U}}
#' \item{\code{pollutant}}{same as in input, one of \code{"NO2"}, \code{"O3"}, \code{"PM10"}, \code{"PM2.5"}}
#' \item{\code{beta}}{same as in input}
#' }}
#' }
#' @export
#' @import dplyr
#'
#' @examples
#' # prepare dataset
#' require(dplyr)
#' Mod <- dMean(mod.data %>% filter(Var=="c_PM10"),
#' value = "Value", time = "Time", point = "Point")
#' Obs <- obs.pm10 %>% mutate(Day=format(Time,"%Y-%m-%d"),
#' Point=ID)
#' Dat <- inner_join(Mod, Obs, by=c("Point", "Day"), suffix = c(".mod", ".obs"))
#'
#' # calculate indicators
#' t_rep <- target_report(Dat, obs = "Value.obs", mod = "Value.mod",
#' point = "Point", pollutant = "PM10")
#'
#' # plot
#' target_plot(t_rep)
#' @rdname target_plot
target_report <- function(data, obs, mod, point,
pollutant=c("NO2", "O3", "PM10", "PM2.5"),
beta=2) {
pollutant <- match.arg(pollutant)
expected <- switch(pollutant,
"NO2"=365*24*0.75,
"PM10"=365*0.75,
"PM2.5"=365*0.75,
"O3"=365*0.75)
data <- as.data.frame(data)
#Mod.ave=Obs.ave=Point=bias_norm=crmse_norm=n_valid=rmsu=NULL
Dat <- data.frame(Obs=data[,obs], Mod=data[,mod], Point=data[,point])
Dat %>% mutate(Obs = Obs+Mod-Mod,
Mod = Mod+Obs-Obs) %>%
filter(!is.na(Obs))-> Dat
Dat %>% group_by(Point) %>%
summarize(Obs.ave = mean(Obs, na.rm=T),
Mod.ave = mean(Mod, na.rm=T),
rmsu = RMS_U_obs(obs = Obs, pollutant = pollutant),
crmse_norm = CRMSE(mod = Mod, obs = Obs)/(beta*rmsu),
r = R(mod = Mod, obs = Obs),
crmse_ratio= abs(NMSD(mod = Mod, obs = Obs))/sqrt(2*(1-r)),
bias_norm = BIAS(mod = Mod, obs = Obs)/(beta*rmsu),
mqi_ts = MQI_ts_synth(mod = Mod, obs = Obs, pollutant = pollutant),
mqi_year = MQI_year(obs = Obs.ave, mod = Mod.ave, pollutant = pollutant),
n_valid = n()) %>%
filter(n_valid>=expected) -> quality_points
quality_overall <- data.frame(mqi_ts_p90 = quantile(quality_points$mqi_ts, probs=0.9, na.rm=T),
mqi_year_p90 = quantile(quality_points$mqi_year, probs=0.9, na.rm=T),
n_points = nrow(quality_points))
out <- list()
out$quality_points <- quality_points
out$quality_overall <- quality_overall
out$parameters <- params_U(pollutant=pollutant)
out$parameters$pollutant <- pollutant
out$parameters$beta <- beta
return(out)
}
#' \code{target_plot}: target plot
#'
#' @import ggplot2
#' @return \code{target_plot} returns a target plot (object of class \code{ggplot} and \code{gg})
#' @rdname target_plot
#' @export
target_plot <- function(t_rep) {
Point=bias_norm=crmse_norm=NULL
ggplot(data = t_rep$quality_points,
aes(x = crmse_norm*sign(crmse_ratio-1), y = bias_norm, shape = Point)) +
geom_point(size=2.5, col="dodgerblue4") +
scale_shape_manual(name='',
values = unlist(lapply(c(LETTERS,letters)[1:length(t_rep$quality_points$Point)], utf8ToInt))) +
stat_function(fun=function(x){+sqrt(1-x^2)}, geom="area", aes(fill="chartreuse3"), n=1001,
alpha=0.4, col=NA, inherit.aes = FALSE)+
stat_function(fun=function(x){-sqrt(1-x^2)}, geom="area", aes(fill="chartreuse3"), n=1001,
alpha=0.4, col=NA, inherit.aes = FALSE)+
scale_fill_manual(values = c('chartreuse3' = 'chartreuse3'),name = '',
labels = expression(RMSE <= beta*RMS[U]))+
stat_function(fun=function(x){+sqrt((1/t_rep$parameters$beta)^2-x^2)}, geom="line",
colour="darkgrey", n=2001, linetype=3, inherit.aes = FALSE)+
stat_function(fun=function(x){-sqrt((1/t_rep$parameters$beta)^2-x^2)}, geom="line",
colour="darkgrey", n=2001, linetype=3, inherit.aes = FALSE) -> p
p +
theme_bw() + theme(legend.position="bottom")+
geom_abline(intercept=0,slope=+1,col="darkgrey")+
geom_abline(intercept=0,slope=-1,col="darkgrey")+
xlim(-2,2) + ylim(-2,2) +
scale_x_continuous(labels=c(2,1,0,1,2))+
xlab(expression(CRMSE[norm]))+ylab(expression(BIAS[norm])) +
coord_fixed() -> p
p +
annotate(geom = "text", x = 0, y = 1.8, label = "BIAS>0" , color = "darkgrey", alpha=0.8)+
annotate(geom = "text", x = 0, y = -1.8, label = "BIAS<0" , color = "darkgrey", alpha=0.8)+
annotate(geom = "text", x = 2, y = 0, label = "St.Dev.", color = "darkgrey", alpha=0.8, adj=1)+
annotate(geom = "text", x = -2, y = 0, label = "R" , color = "darkgrey", alpha=0.8, adj=0)+
annotate(geom = "text", x = -2, y = 1.8, size=6, parse=T,
label = paste0("MQI==",round(t_rep$quality_overall$mqi_ts_p90,2)),
color = ifelse(round(t_rep$quality_overall$mqi_ts_p90,2)>1,
"firebrick3","green4"), adj=0, alpha=0.8)+
annotate(geom = "text", x = -2, y = 1.6, size=4, parse=T,
label = paste0("MQI[year]==",round(t_rep$quality_overall$mqi_year_p90,2)),
color = ifelse(round(t_rep$quality_overall$mqi_year_p90,2)>1,
"firebrick3","green4"), adj=0, alpha=0.8) +
annotate(geom = "text", x = -2, y = -1 , size=3, parse=T,
label = paste0("U['95,r']^{RV}==",t_rep$parameters$U_RV95r),
color = "grey20", adj=0, alpha=0.8) +
annotate(geom = "text", x = -2, y = -1.2, size=3, parse=T,
label = paste0("alpha==",t_rep$parameters$alpha),
color = "grey20", adj=0, alpha=0.8) +
annotate(geom = "text", x = -2, y = -1.4, size=3, parse=T,
label = paste0("RV==",t_rep$parameters$RV),
color = "grey20", adj=0, alpha=0.8) +
annotate(geom = "text", x = -2, y = -1.6, size=3, parse=T,
label = paste0("beta==",t_rep$parameters$beta),
color = "grey20", adj=0, alpha=0.8)-> p
p +
ggtitle(paste("Target plot",t_rep$parameters$pollutant)) -> p
p$layers <- p$layers[c(2:length(p$layers),1)]
return(p)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.