From a69169ded18e4b31e1236fe55539aba26fc4ab45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Wed, 23 Jan 2019 11:58:36 +0100 Subject: [PATCH 01/21] Upload New File --- R/plot_forecasts.R | 372 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 372 insertions(+) create mode 100644 R/plot_forecasts.R diff --git a/R/plot_forecasts.R b/R/plot_forecasts.R new file mode 100644 index 00000000..e8166ca4 --- /dev/null +++ b/R/plot_forecasts.R @@ -0,0 +1,372 @@ +######################################################### +# Barcelona Supercomputing Center # +# llledo@bsc.es - September 2018 # +# Functions to plot several probabilistic forecasts # +######################################################### +library(data.table) +library(ggplot2) +library(reshape2) +library(plyr) + +#======================== +# Function to plot several forecasts for an event, +# either initialized at different moments or by different models. +# Probabilities for extreme categories can be added as hatches. +# Ensemble members can be added as jittered points. +# The observed value can be added as a diamond over the pdf. +#======================== +plot_forecasts <- function(fcst_df,tercile_limits,extreme_limits=NULL,obs=NULL,plotfile=NULL,title="Set a title",varname="Varname (units)",fcst_names=NULL,add_ensmemb=c("above","below","no"),colors=c("ggplot","dst","hydro")) { + + #------------------------ + # Color definitions + #------------------------ + colors <- match.arg(colors) + if(colors=="dst") { + colorFill <- rev(c("#FF764D","#b5b5b5","#33BFD1")) + colorHatch <- c("deepskyblue3","indianred3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue","red") + } else if(colors=="hydro") { + colorFill <- rev(c("#41CBC9","#b5b5b5","#FFAB38")) + colorHatch <- c("darkorange1","deepskyblue3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("darkorange3","blue") + } else if(colors=="ggplot") { + gg_color_hue <- function(n) { + hues = seq(15, 375, length = n + 1) + hcl(h = hues, l = 65, c = 100)[1:n] + } + colorFill <- rev(gg_color_hue(3)) + colorHatch <- c("deepskyblue3","indianred1") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue","red") + } else { stop("Unknown color set") } + + #------------------------ + # Check args + #------------------------ + add_ensmemb <- match.arg(add_ensmemb) + if(length(tercile_limits)!=2) stop("Provide two limits for delimiting tercile categories") + if(tercile_limits[1]>=tercile_limits[2]) stop("The provided tercile limits are in the wrong order") + if(!is.null(extreme_limits)) { + if(length(extreme_limits)!=2) stop("Provide two limits for delimiting extreme categories") + if(extreme_limits[1]>=tercile_limits[1] | extreme_limits[2]<=tercile_limits[2]) stop("The provided extreme limits are not consistent with tercile limits") + } + + #------------------------ + # Set proper fcst names + #------------------------ + if(!is.null(fcst_names)) { + colnames(fcst_df) <- factor(fcst_names,levels=fcst_names) + } + + #------------------------ + # Produce a first plot with the pdf for each init in a panel + #------------------------ + melt_df <- melt(fcst_df,variable.name="init") + plot <- ggplot(melt_df,aes(x=value)) + geom_density(alpha=1,na.rm=T) + coord_flip() + facet_wrap(~init,strip.position="top",nrow=1) + xlim(range(c(obs,density(melt_df$value,na.rm=T)$x))) + ggp <- ggplot_build(plot) + + #------------------------ + # Gather the coordinates of the plots + # together with init and corresponding terciles + #------------------------ + tmp.df <- ggp$data[[1]][,c("x","ymin","ymax","PANEL")] + tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] + tmp.df$tercile <- factor(ifelse(tmp.df$xthr) { + lev[i] <- level + last <- x[i] + } + } + level <- level+1 + last <- 1/0 + for(i in seq(length(x),1,-1)) { + if (lev[i]!=0) { next } + if (last-x[i]>thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level+1 + } + return(lev*thr*sqrt(3)/2) +} + + +#================== +# Hatch one polygon +# Based on base polygon function +#================== +polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven=F) { + halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) + cross <- halfplane[-1L] - halfplane[-length(halfplane)] + does.cross <- cross != 0 + if (!any(does.cross)) + return() + x1 <- x[-length(x)][does.cross] + y1 <- y[-length(y)][does.cross] + x2 <- x[-1L][does.cross] + y2 <- y[-1L][does.cross] + t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * (x2 - x1))) + o <- order(t) + tsort <- t[o] + crossings <- cumsum(cross[does.cross][o]) + if (fillOddEven) + crossings <- crossings%%2 + drawline <- crossings != 0 + lx <- x0 + xd * tsort + ly <- y0 + yd * tsort + lx1 <- lx[-length(lx)][drawline] + ly1 <- ly[-length(ly)][drawline] + lx2 <- lx[-1L][drawline] + ly2 <- ly[-1L][drawline] + return(data.frame(x=lx1,y=ly1,xend=lx2,yend=ly2)) +} + +#================== +# Hatch one polygon +# Based on base polygon function +#================== +polygon.fullhatch <- function(x, y, density, angle, width_units, height_units, inches=c(5,1)) { + x <- c(x, x[1L]) + y <- c(y, y[1L]) + angle <- angle%%180 + #usr <- par("usr") + #inches <- par("pin") + upi <- c(width_units, height_units)/inches + if (upi[1L] < 0) + angle <- 180 - angle + if (upi[2L] < 0) + angle <- 180 - angle + upi <- abs(upi) + xd <- cos(angle/180 * pi) * upi[1L] + yd <- sin(angle/180 * pi) * upi[2L] + hatch_ls <- list() + i <- 1 + if (angle < 45 || angle > 135) { + if (angle < 45) { + first.x <- max(x) + last.x <- min(x) + } + else { + first.x <- min(x) + last.x <- max(x) + } + y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) + x0 <- 0 + y0 <- floor((min(y) - first.x * yd/xd)/y.shift) * y.shift + y.end <- max(y) - last.x * yd/xd + while (y0 < y.end) { + hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i+1 + y0 <- y0 + y.shift + } + } + else { + if (angle < 90) { + first.y <- max(y) + last.y <- min(y) + } + else { + first.y <- min(y) + last.y <- max(y) + } + x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) + x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift + y0 <- 0 + x.end <- max(x) - last.y * xd/yd + while (x0 < x.end) { + hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i+1 + x0 <- x0 + x.shift + } + } + return(do.call("rbind",hatch_ls)) +} -- GitLab From e3222571cd012503a6a7cb35de93fee24d3176f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Thu, 24 Jan 2019 16:46:34 +0100 Subject: [PATCH 02/21] Update plot_forecasts.R to match code style --- R/plot_forecasts.R | 640 +++++++++++++++++++++++---------------------- 1 file changed, 329 insertions(+), 311 deletions(-) diff --git a/R/plot_forecasts.R b/R/plot_forecasts.R index e8166ca4..bd5bf197 100644 --- a/R/plot_forecasts.R +++ b/R/plot_forecasts.R @@ -15,268 +15,287 @@ library(plyr) # Ensemble members can be added as jittered points. # The observed value can be added as a diamond over the pdf. #======================== -plot_forecasts <- function(fcst_df,tercile_limits,extreme_limits=NULL,obs=NULL,plotfile=NULL,title="Set a title",varname="Varname (units)",fcst_names=NULL,add_ensmemb=c("above","below","no"),colors=c("ggplot","dst","hydro")) { - - #------------------------ - # Color definitions - #------------------------ - colors <- match.arg(colors) - if(colors=="dst") { - colorFill <- rev(c("#FF764D","#b5b5b5","#33BFD1")) - colorHatch <- c("deepskyblue3","indianred3") - colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("blue","red") - } else if(colors=="hydro") { - colorFill <- rev(c("#41CBC9","#b5b5b5","#FFAB38")) - colorHatch <- c("darkorange1","deepskyblue3") - colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("darkorange3","blue") - } else if(colors=="ggplot") { - gg_color_hue <- function(n) { - hues = seq(15, 375, length = n + 1) - hcl(h = hues, l = 65, c = 100)[1:n] - } - colorFill <- rev(gg_color_hue(3)) - colorHatch <- c("deepskyblue3","indianred1") - colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("blue","red") - } else { stop("Unknown color set") } +PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", varname="Varname (units)", fcst_names=NULL, add_ensmemb=c("above", "below", "no"), colors=c("ggplot", "s2s4e", "hydro")) { + #------------------------ + # Color definitions + #------------------------ + colors <- match.arg(colors) + if (colors == "s2s4e") { + colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1")) + colorHatch <- c("deepskyblue3", "indianred3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue", "red") + } else if (colors == "hydro") { + colorFill <- rev(c("#41CBC9", "#b5b5b5", "#FFAB38")) + colorHatch <- c("darkorange1", "deepskyblue3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("darkorange3", "blue") + } else if (colors == "ggplot") { + gg_color_hue <- function(n) { + hues=seq(15, 375, length=n+1) + hcl(h=hues, l=65, c=100)[1:n] + } + colorFill <- rev(gg_color_hue(3)) + colorHatch <- c("deepskyblue3", "indianred1") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue", "red") + } else { + stop("Unknown color set") + } - #------------------------ - # Check args - #------------------------ - add_ensmemb <- match.arg(add_ensmemb) - if(length(tercile_limits)!=2) stop("Provide two limits for delimiting tercile categories") - if(tercile_limits[1]>=tercile_limits[2]) stop("The provided tercile limits are in the wrong order") - if(!is.null(extreme_limits)) { - if(length(extreme_limits)!=2) stop("Provide two limits for delimiting extreme categories") - if(extreme_limits[1]>=tercile_limits[1] | extreme_limits[2]<=tercile_limits[2]) stop("The provided extreme limits are not consistent with tercile limits") - } + #------------------------ + # Check input arguments + #------------------------ + add_ensmemb <- match.arg(add_ensmemb) + if (length(tercile_limits) != 2) { + stop("Provide two limits for delimiting tercile categories") + } + if (tercile_limits[1] >= tercile_limits[2]) { + stop("The provided tercile limits are in the wrong order") + } + if (!is.null(extreme_limits)) { + if (length(extreme_limits) != 2) { + stop("Provide two limits for delimiting extreme categories") + } + if (extreme_limits[1] >= tercile_limits[1] | extreme_limits[2] <= tercile_limits[2]) { + stop("The provided extreme limits are not consistent with tercile limits") + } + } - #------------------------ - # Set proper fcst names - #------------------------ - if(!is.null(fcst_names)) { - colnames(fcst_df) <- factor(fcst_names,levels=fcst_names) - } + #------------------------ + # Set proper fcst names + #------------------------ + if (!is.null(fcst_names)) { + colnames(fcst_df) <- factor(fcst_names, levels=fcst_names) + } - #------------------------ - # Produce a first plot with the pdf for each init in a panel - #------------------------ - melt_df <- melt(fcst_df,variable.name="init") - plot <- ggplot(melt_df,aes(x=value)) + geom_density(alpha=1,na.rm=T) + coord_flip() + facet_wrap(~init,strip.position="top",nrow=1) + xlim(range(c(obs,density(melt_df$value,na.rm=T)$x))) - ggp <- ggplot_build(plot) + #------------------------ + # Produce a first plot with the pdf for each init in a panel + #------------------------ + melt_df <- melt(fcst_df, variable.name="init") + plot <- ggplot(melt_df, aes(x=value)) + geom_density(alpha=1, na.rm=T) + coord_flip() + facet_wrap(~init, strip.position="top", nrow=1) + xlim(range(c(obs, density(melt_df$value, na.rm=T)$x))) + ggp <- ggplot_build(plot) - #------------------------ - # Gather the coordinates of the plots - # together with init and corresponding terciles - #------------------------ - tmp.df <- ggp$data[[1]][,c("x","ymin","ymax","PANEL")] - tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] - tmp.df$tercile <- factor(ifelse(tmp.df$xthr) { - lev[i] <- level - last <- x[i] - } - } - level <- level+1 - last <- 1/0 - for(i in seq(length(x),1,-1)) { - if (lev[i]!=0) { next } - if (last-x[i]>thr) { - lev[i] <- level - last <- x[i] - } - } - level <- level+1 - } - return(lev*thr*sqrt(3)/2) + lev <- x * 0 + level <- 1 + while (any(lev == 0)) { + last <- -1/0 + for(i in 1:length(x)) { + if (lev[i] != 0) { + next + } + if (x[i] - last > thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level + 1 + last <- 1 / 0 + for (i in seq(length(x), 1, -1)) { + if (lev[i] != 0) { + next + } + if (last - x[i] > thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level + 1 + } + return(lev * thr * sqrt(3) / 2) } @@ -284,89 +303,88 @@ jitter_ensmemb <- function(x,thr=0.1) { # Hatch one polygon # Based on base polygon function #================== -polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven=F) { - halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) - cross <- halfplane[-1L] - halfplane[-length(halfplane)] - does.cross <- cross != 0 - if (!any(does.cross)) - return() - x1 <- x[-length(x)][does.cross] - y1 <- y[-length(y)][does.cross] - x2 <- x[-1L][does.cross] - y2 <- y[-1L][does.cross] - t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * (x2 - x1))) - o <- order(t) - tsort <- t[o] - crossings <- cumsum(cross[does.cross][o]) - if (fillOddEven) - crossings <- crossings%%2 - drawline <- crossings != 0 - lx <- x0 + xd * tsort - ly <- y0 + yd * tsort - lx1 <- lx[-length(lx)][drawline] - ly1 <- ly[-length(ly)][drawline] - lx2 <- lx[-1L][drawline] - ly2 <- ly[-1L][drawline] - return(data.frame(x=lx1,y=ly1,xend=lx2,yend=ly2)) +.polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven=F) { + halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) + cross <- halfplane[-1L] - halfplane[-length(halfplane)] + does.cross <- cross != 0 + if (!any(does.cross)) { + return() + } + x1 <- x[-length(x)][does.cross] + y1 <- y[-length(y)][does.cross] + x2 <- x[-1L][does.cross] + y2 <- y[-1L][does.cross] + t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * (x2 - x1))) + o <- order(t) + tsort <- t[o] + crossings <- cumsum(cross[does.cross][o]) + if (fillOddEven) { + crossings <- crossings%%2 + } + drawline <- crossings != 0 + lx <- x0 + xd * tsort + ly <- y0 + yd * tsort + lx1 <- lx[-length(lx)][drawline] + ly1 <- ly[-length(ly)][drawline] + lx2 <- lx[-1L][drawline] + ly2 <- ly[-1L][drawline] + return(data.frame(x=lx1,y=ly1,xend=lx2,yend=ly2)) } #================== # Hatch one polygon # Based on base polygon function #================== -polygon.fullhatch <- function(x, y, density, angle, width_units, height_units, inches=c(5,1)) { - x <- c(x, x[1L]) - y <- c(y, y[1L]) - angle <- angle%%180 - #usr <- par("usr") - #inches <- par("pin") - upi <- c(width_units, height_units)/inches - if (upi[1L] < 0) - angle <- 180 - angle - if (upi[2L] < 0) - angle <- 180 - angle - upi <- abs(upi) - xd <- cos(angle/180 * pi) * upi[1L] - yd <- sin(angle/180 * pi) * upi[2L] - hatch_ls <- list() - i <- 1 - if (angle < 45 || angle > 135) { - if (angle < 45) { - first.x <- max(x) - last.x <- min(x) - } - else { - first.x <- min(x) - last.x <- max(x) - } - y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) - x0 <- 0 - y0 <- floor((min(y) - first.x * yd/xd)/y.shift) * y.shift - y.end <- max(y) - last.x * yd/xd - while (y0 < y.end) { - hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) - i <- i+1 - y0 <- y0 + y.shift - } +.polygon.fullhatch <- function(x, y, density, angle, width_units, height_units, inches=c(5, 1)) { + x <- c(x, x[1L]) + y <- c(y, y[1L]) + angle <- angle%%180 + upi <- c(width_units, height_units) / inches + if (upi[1L] < 0) { + angle <- 180 - angle + } + if (upi[2L] < 0) { + angle <- 180 - angle + } + upi <- abs(upi) + xd <- cos(angle/180 * pi) * upi[1L] + yd <- sin(angle/180 * pi) * upi[2L] + hatch_ls <- list() + i <- 1 + if (angle < 45 || angle > 135) { + if (angle < 45) { + first.x <- max(x) + last.x <- min(x) + } else { + first.x <- min(x) + last.x <- max(x) + } + y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) + x0 <- 0 + y0 <- floor((min(y) - first.x * yd/xd) / y.shift) * y.shift + y.end <- max(y) - last.x * yd/xd + while (y0 < y.end) { + hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i+1 + y0 <- y0 + y.shift + } + } else { + if (angle < 90) { + first.y <- max(y) + last.y <- min(y) + } else { + first.y <- min(y) + last.y <- max(y) } - else { - if (angle < 90) { - first.y <- max(y) - last.y <- min(y) - } - else { - first.y <- min(y) - last.y <- max(y) - } - x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) - x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift - y0 <- 0 - x.end <- max(x) - last.y * xd/yd - while (x0 < x.end) { - hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) - i <- i+1 - x0 <- x0 + x.shift - } + x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) + x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift + y0 <- 0 + x.end <- max(x) - last.y * xd/yd + while (x0 < x.end) { + hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i+1 + x0 <- x0 + x.shift } - return(do.call("rbind",hatch_ls)) + } + return(do.call("rbind", hatch_ls)) } -- GitLab From b8632e8af570b4c4a97657a663c8c42ca1d92e75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Thu, 24 Jan 2019 16:49:16 +0100 Subject: [PATCH 03/21] Update plot_forecasts.R to call hidden functions correctly --- R/plot_forecasts.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/plot_forecasts.R b/R/plot_forecasts.R index bd5bf197..7e4550dc 100644 --- a/R/plot_forecasts.R +++ b/R/plot_forecasts.R @@ -101,7 +101,7 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NU tmp.df$extremes <- factor(ifelse(tmp.df$x Date: Thu, 24 Jan 2019 17:17:34 +0100 Subject: [PATCH 04/21] More style changes --- R/plot_forecasts.R | 71 +++++++++++++++++++++++++++++++++------------- 1 file changed, 52 insertions(+), 19 deletions(-) diff --git a/R/plot_forecasts.R b/R/plot_forecasts.R index 7e4550dc..810cd57b 100644 --- a/R/plot_forecasts.R +++ b/R/plot_forecasts.R @@ -15,9 +15,13 @@ library(plyr) # Ensemble members can be added as jittered points. # The observed value can be added as a diamond over the pdf. #======================== -PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", varname="Varname (units)", fcst_names=NULL, add_ensmemb=c("above", "below", "no"), colors=c("ggplot", "s2s4e", "hydro")) { +PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, + obs=NULL, plotfile=NULL, title="Set a title", + varname="Varname (units)", fcst_names=NULL, + add_ensmemb=c("above", "below", "no"), + colors=c("ggplot", "s2s4e", "hydro")) { #------------------------ - # Color definitions + # Define color sets #------------------------ colors <- match.arg(colors) if (colors == "s2s4e") { @@ -60,7 +64,8 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NU if (length(extreme_limits) != 2) { stop("Provide two limits for delimiting extreme categories") } - if (extreme_limits[1] >= tercile_limits[1] | extreme_limits[2] <= tercile_limits[2]) { + if (extreme_limits[1] >= tercile_limits[1] | + extreme_limits[2] <= tercile_limits[2]) { stop("The provided extreme limits are not consistent with tercile limits") } } @@ -76,7 +81,11 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NU # Produce a first plot with the pdf for each init in a panel #------------------------ melt_df <- melt(fcst_df, variable.name="init") - plot <- ggplot(melt_df, aes(x=value)) + geom_density(alpha=1, na.rm=T) + coord_flip() + facet_wrap(~init, strip.position="top", nrow=1) + xlim(range(c(obs, density(melt_df$value, na.rm=T)$x))) + plot <- ggplot(melt_df, aes(x=value)) + + geom_density(alpha=1, na.rm=T) + + coord_flip() + + facet_wrap(~init, strip.position="top", nrow=1) + + xlim(range(c(obs, density(melt_df$value, na.rm=T)$x))) ggp <- ggplot_build(plot) #------------------------ @@ -85,7 +94,9 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NU #------------------------ tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] - tmp.df$tercile <- factor(ifelse(tmp.df$x Date: Thu, 24 Jan 2019 19:51:12 +0100 Subject: [PATCH 05/21] more code style changes --- R/plot_forecasts.R | 58 +++++++++++++++++++++++++++++----------------- 1 file changed, 37 insertions(+), 21 deletions(-) diff --git a/R/plot_forecasts.R b/R/plot_forecasts.R index 810cd57b..7c2b404c 100644 --- a/R/plot_forecasts.R +++ b/R/plot_forecasts.R @@ -214,15 +214,16 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, if (!is.null(obs)) { plot <- plot + # this adds the obs diamond - geom_point(data=obs_xy, aes(x=x, y=ymax, size="Observation"), shape=23, - color="black", fill=colorObs) + geom_point(data=obs_xy, aes(x=x, y=ymax, size="Observation"), + shape=23, color="black", fill=colorObs) } #------------------------ # Compute probability for each tercile and identify MLT #------------------------ tmp.dt <- data.table(tmp.df) - pct <- tmp.dt[, .(pct=integrate(approxfun(x, ymax), lower=min(x), upper=max(x))$value), by=.(init, tercile)] + pct <- tmp.dt[, .(pct=integrate(approxfun(x, ymax), lower=min(x), upper=max(x))$value), + by=.(init, tercile)] tot <- pct[, .(tot=sum(pct)), by=init] pct <- merge(pct, tot, by="init") pct$pct <- round(100 * pct$pct / pct$tot, 0) @@ -233,14 +234,16 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, # Compute probability for extremes #------------------------ if (!is.null(extreme_limits)) { - pct2 <- tmp.dt[, .(pct=integrate(approxfun(x, ymax), lower=min(x), upper=max(x))$value), by=.(init, extremes)] + pct2 <- tmp.dt[, .(pct=integrate(approxfun(x, ymax), lower=min(x), upper=max(x))$value), + by=.(init, extremes)] tot2 <- pct2[, .(tot=sum(pct)), by=init] pct2 <- merge(pct2, tot2, by="init") pct2$pct <- round(100 * pct2$pct / pct2$tot, 0) lab_pos2 <- c(extreme_limits[1], NA, extreme_limits[2]) pct2 <- merge(pct2, max.df, by=c("init", "extremes")) # include potentially missing groups - pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), levels=c("Below P10", "Normal", "Above P90"))), ] + pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), + levels=c("Below P10", "Normal", "Above P90"))), ] } #------------------------ @@ -254,15 +257,20 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, vjust <- -0.5 } plot <- plot + - geom_text(data=pct, aes(x=lab_pos[as.integer(tercile)], y=labpos, label=paste0(pct, "%"), hjust=as.integer(tercile)*1.5-2.5), vjust=vjust, angle=-90, size=3.2) + - geom_text(data=pct[MLT==T, ], aes(x=lab_pos[as.integer(tercile)], y=labpos, label="*", hjust=as.integer(tercile)*3.5-5.0), vjust=0.1, angle=-90, size=7, color="black") + geom_text(data=pct, aes(x=lab_pos[as.integer(tercile)], y=labpos, label=paste0(pct, "%"), + hjust=as.integer(tercile) * 1.5 - 2.5), vjust=vjust, angle=-90, size=3.2) + + geom_text(data=pct[MLT==T, ], aes(x=lab_pos[as.integer(tercile)], y=labpos, + label="*", hjust=as.integer(tercile) * 3.5 - 5.0), vjust=0.1, angle=-90, + size=7, color="black") #------------------------ # Add probability labels for extremes #------------------------ if (!is.null(extreme_limits)) { plot <- plot + - geom_text(data=pct2[extremes!="Normal", ], aes(x=lab_pos2[as.integer(extremes)], y=0.9*y, label=paste0(pct, "%"), hjust=as.integer(extremes)*1.5-2.5), vjust=-0.5, angle=-90, size=3.2, color=rep(colorLab, dim(fcst_df)[2])) + geom_text(data=pct2[extremes!="Normal", ], aes(x=lab_pos2[as.integer(extremes)], + y=0.9*y, label=paste0(pct, "%"), hjust=as.integer(extremes) * 1.5 - 2.5), + vjust=-0.5, angle=-90, size=3.2, color=rep(colorLab, dim(fcst_df)[2])) } #------------------------ @@ -270,13 +278,21 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, #------------------------ plot <- plot + theme_minimal() + - scale_fill_manual(name="Probability of\nterciles", breaks=c("Above normal", "Normal", "Below normal"), values=colorFill, drop=F) + + scale_fill_manual(name="Probability of\nterciles", + breaks=c("Above normal", "Normal", "Below normal"), + values=colorFill, drop=F) + scale_color_manual(name="Probability of\nextremes", values=colorHatch) + scale_shape_manual(name="Ensemble\nmembers", values=c(21)) + scale_size_manual(name="Observation", values=c(3)) + labs(x=varname, y="Probabilty density\n(total area=1)", title=title) + - theme(axis.text.x=element_blank(), panel.grid.minor.x=element_blank(), legend.key.size=unit(0.3, "in"), panel.border=element_rect(fill=NA, color="black"), strip.background=element_rect(colour="black", fill="gray80"), panel.spacing=unit(0.2, "in"), panel.grid.major.x=element_line(color="grey93")) + - guides(fill=guide_legend(order=1), color=guide_legend(order=2, reverse=T), shape=guide_legend(order=3, label=F), size=guide_legend(order=4, label=F)) + theme(axis.text.x=element_blank(), panel.grid.minor.x=element_blank(), + legend.key.size=unit(0.3, "in"), + panel.border=element_rect(fill=NA, color="black"), + strip.background=element_rect(colour="black", fill="gray80"), + panel.spacing=unit(0.2, "in"), + panel.grid.major.x=element_line(color="grey93")) + + guides(fill=guide_legend(order=1), color=guide_legend(order=2, reverse=T), + shape=guide_legend(order=3, label=F), size=guide_legend(order=4, label=F)) #------------------------ # Save to plotfile if needed, and return plot @@ -306,7 +322,7 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, level <- 1 while (any(lev == 0)) { last <- -1/0 - for(i in 1:length(x)) { + for (i in 1:length(x)) { if (lev[i] != 0) { next } @@ -361,7 +377,7 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, ly1 <- ly[-length(ly)][drawline] lx2 <- lx[-1L][drawline] ly2 <- ly[-1L][drawline] - return(data.frame(x=lx1,y=ly1,xend=lx2,yend=ly2)) + return(data.frame(x=lx1, y=ly1, xend=lx2, yend=ly2)) } #================== @@ -380,8 +396,8 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, angle <- 180 - angle } upi <- abs(upi) - xd <- cos(angle/180 * pi) * upi[1L] - yd <- sin(angle/180 * pi) * upi[2L] + xd <- cos(angle / 180 * pi) * upi[1L] + yd <- sin(angle / 180 * pi) * upi[2L] hatch_ls <- list() i <- 1 if (angle < 45 || angle > 135) { @@ -392,13 +408,13 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, first.x <- min(x) last.x <- max(x) } - y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) + y.shift <- upi[2L] / density / abs(cos(angle / 180 * pi)) x0 <- 0 y0 <- floor((min(y) - first.x * yd/xd) / y.shift) * y.shift y.end <- max(y) - last.x * yd/xd while (y0 < y.end) { hatch_ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) - i <- i+1 + i <- i + 1 y0 <- y0 + y.shift } } else { @@ -409,13 +425,13 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, first.y <- min(y) last.y <- max(y) } - x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) - x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift + x.shift <- upi[1L] / density / abs(sin(angle / 180 * pi)) + x0 <- floor((min(x) - first.y * xd / yd) / x.shift) * x.shift y0 <- 0 - x.end <- max(x) - last.y * xd/yd + x.end <- max(x) - last.y * xd / yd while (x0 < x.end) { hatch_ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) - i <- i+1 + i <- i + 1 x0 <- x0 + x.shift } } -- GitLab From e9e034fab7a45db49d201e9ad601ab58f8a0cf60 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Thu, 24 Jan 2019 19:55:20 +0100 Subject: [PATCH 06/21] Rename file --- R/{plot_forecasts.R => plotForecastPDF.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{plot_forecasts.R => plotForecastPDF.R} (100%) diff --git a/R/plot_forecasts.R b/R/plotForecastPDF.R similarity index 100% rename from R/plot_forecasts.R rename to R/plotForecastPDF.R -- GitLab From 0e2d88842d8941b6e143e1d48f912f0c46734a21 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Thu, 24 Jan 2019 19:56:05 +0100 Subject: [PATCH 07/21] Rename file --- R/{plotForecastPDF.R => PlotForecastPDF.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{plotForecastPDF.R => PlotForecastPDF.R} (100%) diff --git a/R/plotForecastPDF.R b/R/PlotForecastPDF.R similarity index 100% rename from R/plotForecastPDF.R rename to R/PlotForecastPDF.R -- GitLab From 64b9ac9bdfe76445fae2f18b9da3ecdd39e2cd73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Thu, 24 Jan 2019 20:40:43 +0100 Subject: [PATCH 08/21] Added the documentation (TBC) --- R/PlotForecastPDF.R | 43 +++++++++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 7c2b404c..77f6bedd 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -1,20 +1,31 @@ -######################################################### -# Barcelona Supercomputing Center # -# llledo@bsc.es - September 2018 # -# Functions to plot several probabilistic forecasts # -######################################################### -library(data.table) -library(ggplot2) -library(reshape2) -library(plyr) -#======================== -# Function to plot several forecasts for an event, -# either initialized at different moments or by different models. -# Probabilities for extreme categories can be added as hatches. -# Ensemble members can be added as jittered points. -# The observed value can be added as a diamond over the pdf. -#======================== +#'Plot one or multiple ensemble forecast pdfs for the same event +#' +#'@author Llorenç Lledó \email{llledo@bsc.es} +#'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. +#' +#'@param fcst_df a dataframe blabla TODO +#'@param tercile_limits an array with P33 and P66 values that define the tercile categories. +#'@param extreme_limits (optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown). +#'@param obs (optional) a number with the observed value. (Default: observation is not shown). +#'@param plotfile (optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved). +#'@param title a string with the plot title. +#'@param varname a string with the variable name and units. +#'@param fcst_names (optional) an array of strings with the titles of each individual forecast. +#'@param add_ensmemb either to add the ensemble members (above or below the pdf) or not. +#'@param colors a selection of predefined color sets: ggplot for blue/green/red, s2s4e for , or hydro for a suitable scale for precipitation/inflows. +#' +#'@return a ggplot object containing the plot. +#' +#'@import data.table +#'@import ggplot2 +#'@import reshape2 +#'@import plyr +#' +#'@examples +#' Blabla +#' +#'@export PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", varname="Varname (units)", fcst_names=NULL, -- GitLab From 7bc0cb346df4b48dfa8461e645a24793a75d3aca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 25 Jan 2019 13:11:12 +0100 Subject: [PATCH 09/21] Enhancing documentation and code style. --- R/PlotForecastPDF.R | 78 ++++++++++++++++++++++----------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 77f6bedd..8d643ae8 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -4,16 +4,16 @@ #'@author Llorenç Lledó \email{llledo@bsc.es} #'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. #' -#'@param fcst_df a dataframe blabla TODO -#'@param tercile_limits an array with P33 and P66 values that define the tercile categories. -#'@param extreme_limits (optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown). +#'@param fcst a dataframe or array containing all the ensember members for each frecast. If \code{"fcst"} is an array, it should have two labelled dimensions, and one of them should be \code{"members"}. If \code{"fcsts"} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. +#'@param tercile.limits an array with P33 and P66 values that define the tercile categories. +#'@param extreme.limits (optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown). #'@param obs (optional) a number with the observed value. (Default: observation is not shown). #'@param plotfile (optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved). #'@param title a string with the plot title. -#'@param varname a string with the variable name and units. -#'@param fcst_names (optional) an array of strings with the titles of each individual forecast. -#'@param add_ensmemb either to add the ensemble members (above or below the pdf) or not. -#'@param colors a selection of predefined color sets: ggplot for blue/green/red, s2s4e for , or hydro for a suitable scale for precipitation/inflows. +#'@param var.name a string with the variable name and units. +#'@param fcst.names (optional) an array of strings with the titles of each individual forecast. +#'@param add.ensmemb either to add the ensemble members (above or below the pdf) or not. +#'@param color.set a selection of predefined color sets: use \code{"ggplot"} (default) for blue/green/red, \code{"s2s4e"} for blue/grey/orange, or \code{"hydro"} for yellow/gray/blue (suitable for precipitation and inflows). #' #'@return a ggplot object containing the plot. #' @@ -26,28 +26,28 @@ #' Blabla #' #'@export -PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, +PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", - varname="Varname (units)", fcst_names=NULL, - add_ensmemb=c("above", "below", "no"), - colors=c("ggplot", "s2s4e", "hydro")) { + var.name="Varname (units)", fcst.names=NULL, + add.ensmemb=c("above", "below", "no"), + color.set=c("ggplot", "s2s4e", "hydro")) { #------------------------ # Define color sets #------------------------ - colors <- match.arg(colors) - if (colors == "s2s4e") { + color.set <- match.arg(color.set) + if (color.set == "s2s4e") { colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1")) colorHatch <- c("deepskyblue3", "indianred3") colorMember <- c("#ffff7f") colorObs <- "purple" colorLab <- c("blue", "red") - } else if (colors == "hydro") { + } else if (color.set == "hydro") { colorFill <- rev(c("#41CBC9", "#b5b5b5", "#FFAB38")) colorHatch <- c("darkorange1", "deepskyblue3") colorMember <- c("#ffff7f") colorObs <- "purple" colorLab <- c("darkorange3", "blue") - } else if (colors == "ggplot") { + } else if (color.set == "ggplot") { gg_color_hue <- function(n) { hues=seq(15, 375, length=n+1) hcl(h=hues, l=65, c=100)[1:n] @@ -64,19 +64,19 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, #------------------------ # Check input arguments #------------------------ - add_ensmemb <- match.arg(add_ensmemb) - if (length(tercile_limits) != 2) { + add.ensmemb <- match.arg(add.ensmemb) + if (length(tercile.limits) != 2) { stop("Provide two limits for delimiting tercile categories") } - if (tercile_limits[1] >= tercile_limits[2]) { + if (tercile.limits[1] >= tercile.limits[2]) { stop("The provided tercile limits are in the wrong order") } - if (!is.null(extreme_limits)) { - if (length(extreme_limits) != 2) { + if (!is.null(extreme.limits)) { + if (length(extreme.limits) != 2) { stop("Provide two limits for delimiting extreme categories") } - if (extreme_limits[1] >= tercile_limits[1] | - extreme_limits[2] <= tercile_limits[2]) { + if (extreme.limits[1] >= tercile.limits[1] | + extreme.limits[2] <= tercile.limits[2]) { stop("The provided extreme limits are not consistent with tercile limits") } } @@ -84,8 +84,8 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, #------------------------ # Set proper fcst names #------------------------ - if (!is.null(fcst_names)) { - colnames(fcst_df) <- factor(fcst_names, levels=fcst_names) + if (!is.null(fcst.names)) { + colnames(fcst_df) <- factor(fcst.names, levels=fcst.names) } #------------------------ @@ -105,8 +105,8 @@ PlotForecastPDF <- function(fcst_df, tercile_limits, extreme_limits=NULL, #------------------------ tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] - tmp.df$tercile <- factor(ifelse(tmp.df$x Date: Fri, 25 Jan 2019 13:48:14 +0100 Subject: [PATCH 10/21] Adapt input to read arrays. Enhancing code style. --- R/PlotForecastPDF.R | 103 ++++++++++++++++++++++++++------------------ 1 file changed, 61 insertions(+), 42 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 8d643ae8..6ef00d0a 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -48,17 +48,17 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, colorObs <- "purple" colorLab <- c("darkorange3", "blue") } else if (color.set == "ggplot") { - gg_color_hue <- function(n) { + gg.color.hue <- function(n) { hues=seq(15, 375, length=n+1) hcl(h=hues, l=65, c=100)[1:n] } - colorFill <- rev(gg_color_hue(3)) + colorFill <- rev(gg.color.hue(3)) colorHatch <- c("deepskyblue3", "indianred1") colorMember <- c("#ffff7f") colorObs <- "purple" colorLab <- c("blue", "red") } else { - stop("Unknown color set") + stop("Parameter 'color.set' should be one of ggplot/s2s4e/hydro") } #------------------------ @@ -66,14 +66,14 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ add.ensmemb <- match.arg(add.ensmemb) if (length(tercile.limits) != 2) { - stop("Provide two limits for delimiting tercile categories") + stop("Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") } if (tercile.limits[1] >= tercile.limits[2]) { stop("The provided tercile limits are in the wrong order") } if (!is.null(extreme.limits)) { if (length(extreme.limits) != 2) { - stop("Provide two limits for delimiting extreme categories") + stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") } if (extreme.limits[1] >= tercile.limits[1] | extreme.limits[2] <= tercile.limits[2]) { @@ -81,22 +81,41 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, } } + #------------------------ + # Check fcst type and convert to data.frame if needed + #------------------------ + if (is.array(fcst) { + if (!"members" %in% names(dim(fcst1)) | length(dim(fcst)) != 2) { + stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be 'members'") + } + dim.members <- which(names(dim(fcst)) == "members") + if (dim.members == 1) { + fcst.df <- data.frame(fcst) + } else { + fcst.df <- data.frame(t(fcst)) + } + } else if (is.data.frame(fcst)) { + fcst.df <- fcst + } else { + stop("Parameter 'fcst' should be an array or a data.frame") + } + #------------------------ # Set proper fcst names #------------------------ if (!is.null(fcst.names)) { - colnames(fcst_df) <- factor(fcst.names, levels=fcst.names) + colnames(fcst.df) <- factor(fcst.names, levels=fcst.names) } #------------------------ # Produce a first plot with the pdf for each init in a panel #------------------------ - melt_df <- melt(fcst_df, variable.name="init") - plot <- ggplot(melt_df, aes(x=value)) + + melt.df <- melt(fcst.df, variable.name="init") + plot <- ggplot(melt.df, aes(x=value)) + geom_density(alpha=1, na.rm=T) + coord_flip() + facet_wrap(~init, strip.position="top", nrow=1) + - xlim(range(c(obs, density(melt_df$value, na.rm=T)$x))) + xlim(range(c(obs, density(melt.df$value, na.rm=T)$x))) ggp <- ggplot_build(plot) #------------------------ @@ -112,9 +131,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ # Get the height and width of a panel #------------------------ - pan_width <- diff(range(tmp.df$x)) - pan_height <- max(tmp.df$ymax) - magicratio <- 9 * pan_height / pan_width + pan.width <- diff(range(tmp.df$x)) + pan.height <- max(tmp.df$ymax) + magic.ratio <- 9 * pan.height / pan.width #------------------------ # Compute hatch coordinates for extremes @@ -128,7 +147,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, tmp.df2 <- data.frame(x=c(x$x, max(x$x), min(x$x)), y=c(x$ymax, 0, 0)) # compute the hatches for this polygon hatches <- .polygon.fullhatch(tmp.df2$x, tmp.df2$y, angle=60, density=10, - width_units=pan_width, height_units=pan_height) + width_units=pan.width, height_units=pan.height) # add bottom segment end1 <- data.frame(x=x$x[1], y=x$ymax[1], xend=x$x[1], yend=0) # add top segment @@ -144,7 +163,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Compute max y for each extreme category max.ls <- dlply(tmp.df, .(init,extremes), function(x) { - data.frame(y=min(0.85 * pan_height, max(x$ymax))) + data.frame(y=min(0.85 * pan.height, max(x$ymax))) }) attr <- attr(max.ls, "split_labels") for (i in 1:length(max.ls)) { @@ -156,10 +175,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ # Compute jitter space for ensemble members #------------------------ - jitter_df <- melt(data.frame(dlply(melt_df, .(init), function(x) { - .jitter_ensmemb(sort(x$value), pan_width / 100) }), check.names=F), + jitter.df <- melt(data.frame(dlply(melt.df, .(init), function(x) { + .jitter.ensmemb(sort(x$value), pan.width / 100) }), check.names=F), value.name="yjitter", variable.name="init") - jitter_df$x <- melt(data.frame(dlply(melt_df, .(init), function(x) { + jitter.df$x <- melt(data.frame(dlply(melt.df, .(init), function(x) { sort(x$value)}) ), value.name="x")$x #------------------------ @@ -168,10 +187,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ if (!is.null(obs)) { tmp.dt <- data.table(tmp.df, key=c("init", "x")) - obs_dt <- data.table(init=factor(colnames(fcst_df), levels=colnames(fcst_df)), - value=rep(obs, dim(fcst_df)[2])) - setkey(obs_dt, init, value) - obs_xy <- tmp.dt[obs_dt, roll="nearest"] + obs.dt <- data.table(init=factor(colnames(fcst.df), levels=colnames(fcst.df)), + value=rep(obs, dim(fcst.df)[2])) + setkey(obs.dt, init, value) + obs.xy <- tmp.dt[obs.dt, roll="nearest"] } #------------------------ @@ -199,7 +218,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ if (!is.null(obs)) { plot <- plot + - geom_vline(data=obs_dt, aes(xintercept=value), linetype="dashed", color=colorObs) + geom_vline(data=obs.dt, aes(xintercept=value), linetype="dashed", color=colorObs) } #------------------------ @@ -208,15 +227,15 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, if (add.ensmemb == "below") { plot <- plot + # this adds a grey box for ensmembers - geom_rect(aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=-pan_height/10), + geom_rect(aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=-pan.height/10), fill="gray95", color="black", width=0.2) + # this adds the ensemble members - geom_point(data=jitter_df, color="black", fill=colorMember, alpha=1, - aes(x=x, y=-pan_height/10-magicratio*yjitter, shape="Ensemble members")) + geom_point(data=jitter.df, color="black", fill=colorMember, alpha=1, + aes(x=x, y= -pan.height / 10 - magic.ratio * yjitter, shape="Ensemble members")) } else if (add.ensmemb == "above") { plot <- plot + - geom_point(data=jitter_df, color="black", fill=colorMember, alpha=1, - aes(x=x, y=0.7*magicratio*yjitter, shape="Ensemble members")) + geom_point(data=jitter.df, color="black", fill=colorMember, alpha=1, + aes(x=x, y=0.7 * magic.ratio * yjitter, shape="Ensemble members")) } #------------------------ @@ -225,7 +244,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, if (!is.null(obs)) { plot <- plot + # this adds the obs diamond - geom_point(data=obs_xy, aes(x=x, y=ymax, size="Observation"), + geom_point(data=obs.xy, aes(x=x, y=ymax, size="Observation"), shape=23, color="black", fill=colorObs) } @@ -239,7 +258,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, pct <- merge(pct, tot, by="init") pct$pct <- round(100 * pct$pct / pct$tot, 0) pct$MLT <- pct[, .(MLT=pct==max(pct)), by=init]$MLT - lab_pos <- c(tercile.limits[1], mean(tercile.limits), tercile.limits[2]) + lab.pos <- c(tercile.limits[1], mean(tercile.limits), tercile.limits[2]) #------------------------ # Compute probability for extremes @@ -250,7 +269,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, tot2 <- pct2[, .(tot=sum(pct)), by=init] pct2 <- merge(pct2, tot2, by="init") pct2$pct <- round(100 * pct2$pct / pct2$tot, 0) - lab_pos2 <- c(extreme.limits[1], NA, extreme.limits[2]) + lab.pos2 <- c(extreme.limits[1], NA, extreme.limits[2]) pct2 <- merge(pct2, max.df, by=c("init", "extremes")) # include potentially missing groups pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), @@ -261,16 +280,16 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Add probability labels for terciles #------------------------ if (add.ensmemb == "above") { - labpos <- -0.2 * pan_height + labpos <- -0.2 * pan.height vjust <- 0 } else { labpos <- 0 vjust <- -0.5 } plot <- plot + - geom_text(data=pct, aes(x=lab_pos[as.integer(tercile)], y=labpos, label=paste0(pct, "%"), + geom_text(data=pct, aes(x=lab.pos[as.integer(tercile)], y=labpos, label=paste0(pct, "%"), hjust=as.integer(tercile) * 1.5 - 2.5), vjust=vjust, angle=-90, size=3.2) + - geom_text(data=pct[MLT==T, ], aes(x=lab_pos[as.integer(tercile)], y=labpos, + geom_text(data=pct[MLT==T, ], aes(x=lab.pos[as.integer(tercile)], y=labpos, label="*", hjust=as.integer(tercile) * 3.5 - 5.0), vjust=0.1, angle=-90, size=7, color="black") @@ -279,9 +298,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ if (!is.null(extreme.limits)) { plot <- plot + - geom_text(data=pct2[extremes!="Normal", ], aes(x=lab_pos2[as.integer(extremes)], + geom_text(data=pct2[extremes!="Normal", ], aes(x=lab.pos2[as.integer(extremes)], y=0.9*y, label=paste0(pct, "%"), hjust=as.integer(extremes) * 1.5 - 2.5), - vjust=-0.5, angle=-90, size=3.2, color=rep(colorLab, dim(fcst_df)[2])) + vjust=-0.5, angle=-90, size=3.2, color=rep(colorLab, dim(fcst.df)[2])) } #------------------------ @@ -319,7 +338,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # A function to distribute ensemble members so not to overlap # Requires a sorted array of values. #================== -.jitter_ensmemb <- function(x, thr=0.1) { +.jitter.ensmemb <- function(x, thr=0.1) { # Idea: start with first level. Loop all points, # and if distance to last point in the level is more than a threshold, # include the point to the level. @@ -395,11 +414,11 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Hatch one polygon # Based on base polygon function #================== -.polygon.fullhatch <- function(x, y, density, angle, width_units, height_units, inches=c(5, 1)) { +.polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches=c(5, 1)) { x <- c(x, x[1L]) y <- c(y, y[1L]) angle <- angle%%180 - upi <- c(width_units, height_units) / inches + upi <- c(width.units, height.units) / inches if (upi[1L] < 0) { angle <- 180 - angle } @@ -409,7 +428,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, upi <- abs(upi) xd <- cos(angle / 180 * pi) * upi[1L] yd <- sin(angle / 180 * pi) * upi[2L] - hatch_ls <- list() + hatch.ls <- list() i <- 1 if (angle < 45 || angle > 135) { if (angle < 45) { @@ -424,7 +443,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, y0 <- floor((min(y) - first.x * yd/xd) / y.shift) * y.shift y.end <- max(y) - last.x * yd/xd while (y0 < y.end) { - hatch_ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) + hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) i <- i + 1 y0 <- y0 + y.shift } @@ -441,10 +460,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, y0 <- 0 x.end <- max(x) - last.y * xd / yd while (x0 < x.end) { - hatch_ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) + hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) i <- i + 1 x0 <- x0 + x.shift } } - return(do.call("rbind", hatch_ls)) + return(do.call("rbind", hatch.ls)) } -- GitLab From 33ce5019f9005baa4d67fa7c6f640678c76a1179 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 25 Jan 2019 14:35:25 +0100 Subject: [PATCH 11/21] Included examples, corrected bugs --- R/PlotForecastPDF.R | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 6ef00d0a..55569cdf 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -23,7 +23,11 @@ #'@import plyr #' #'@examples -#' Blabla +#'fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) +#'PlotForecastPDF(fcsts,c(-1,1)) +#'fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) +#'PlotForecastPDF(fcsts2,c(-1,1)) +#' #' #'@export PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, @@ -84,7 +88,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ # Check fcst type and convert to data.frame if needed #------------------------ - if (is.array(fcst) { + if (is.array(fcst)) { if (!"members" %in% names(dim(fcst1)) | length(dim(fcst)) != 2) { stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be 'members'") } @@ -123,9 +127,15 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # together with init and corresponding terciles #------------------------ tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] - tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] - tmp.df$tercile <- factor(ifelse(tmp.df$x Date: Fri, 25 Jan 2019 14:52:34 +0100 Subject: [PATCH 12/21] Correct small bugs --- R/PlotForecastPDF.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 55569cdf..592e8325 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -26,7 +26,7 @@ #'fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) #'PlotForecastPDF(fcsts,c(-1,1)) #'fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) -#'PlotForecastPDF(fcsts2,c(-1,1)) +#'PlotForecastPDF(fcsts2,c(-0.66,0.66)) #' #' #'@export @@ -89,7 +89,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Check fcst type and convert to data.frame if needed #------------------------ if (is.array(fcst)) { - if (!"members" %in% names(dim(fcst1)) | length(dim(fcst)) != 2) { + if (!"members" %in% names(dim(fcst)) | length(dim(fcst)) != 2) { stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be 'members'") } dim.members <- which(names(dim(fcst)) == "members") @@ -157,7 +157,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, tmp.df2 <- data.frame(x=c(x$x, max(x$x), min(x$x)), y=c(x$ymax, 0, 0)) # compute the hatches for this polygon hatches <- .polygon.fullhatch(tmp.df2$x, tmp.df2$y, angle=60, density=10, - width_units=pan.width, height_units=pan.height) + width.units=pan.width, height.units=pan.height) # add bottom segment end1 <- data.frame(x=x$x[1], y=x$ymax[1], xend=x$x[1], yend=0) # add top segment -- GitLab From 580ff962ac06c2bbaed1f4ae992e575821789833 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 25 Jan 2019 15:02:44 +0100 Subject: [PATCH 13/21] Expanded examples, added extra check for fcst.names --- R/PlotForecastPDF.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 592e8325..7dc092c8 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -26,7 +26,7 @@ #'fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) #'PlotForecastPDF(fcsts,c(-1,1)) #'fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) -#'PlotForecastPDF(fcsts2,c(-0.66,0.66)) +#'PlotForecastPDF(fcsts2,c(-0.66,0.66),extreme.limits=c(-1.2,1.2),fcst.names=paste0("random fcst ",1:5),obs=0.7) #' #' #'@export @@ -108,6 +108,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Set proper fcst names #------------------------ if (!is.null(fcst.names)) { + if (length(fcst.names) != dim(fcst.df)[2]) { + stop("Parameter 'fcst.names' should be an array with as many names as distinct forecasts provided") + } colnames(fcst.df) <- factor(fcst.names, levels=fcst.names) } -- GitLab From a57697e77a9e69ae7ba1d7aaf9f49b3ebc5cbd5f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 25 Jan 2019 15:05:50 +0100 Subject: [PATCH 14/21] Cosmetic change --- R/PlotForecastPDF.R | 1 - 1 file changed, 1 deletion(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 7dc092c8..c78b1436 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -1,4 +1,3 @@ - #'Plot one or multiple ensemble forecast pdfs for the same event #' #'@author Llorenç Lledó \email{llledo@bsc.es} -- GitLab From 41aa28f75939b645dc27462a52f2e326c63cac16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 8 Feb 2019 16:33:12 +0100 Subject: [PATCH 15/21] Remove messages from melt --- R/PlotForecastPDF.R | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index c78b1436..3e946b3c 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -116,7 +116,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ # Produce a first plot with the pdf for each init in a panel #------------------------ - melt.df <- melt(fcst.df, variable.name="init") + melt.df <- melt(fcst.df, variable.name="init", id.vars=NULL) plot <- ggplot(melt.df, aes(x=value)) + geom_density(alpha=1, na.rm=T) + coord_flip() + @@ -187,11 +187,13 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ # Compute jitter space for ensemble members #------------------------ - jitter.df <- melt(data.frame(dlply(melt.df, .(init), function(x) { - .jitter.ensmemb(sort(x$value), pan.width / 100) }), check.names=F), - value.name="yjitter", variable.name="init") - jitter.df$x <- melt(data.frame(dlply(melt.df, .(init), function(x) { - sort(x$value)}) ), value.name="x")$x + if (add.ensmemb != "no") { + jitter.df <- melt(data.frame(dlply(melt.df, .(init), function(x) { + .jitter.ensmemb(sort(x$value), pan.width / 100) }), check.names=F), + value.name="yjitter", variable.name="init", id.vars=NULL) + jitter.df$x <- melt(data.frame(dlply(melt.df, .(init), function(x) { + sort(x$value)}) ), value.name="x", id.vars=NULL)$x + } #------------------------ # Get y coordinates for observed x values, -- GitLab From f909c99d24326a7b6e61a793426c99933cc77aea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Wed, 6 Mar 2019 17:34:41 +0100 Subject: [PATCH 16/21] Fix for plotting ensmembers when there are NAs in a fcst. --- R/PlotForecastPDF.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 3e946b3c..cca7c916 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -189,10 +189,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ if (add.ensmemb != "no") { jitter.df <- melt(data.frame(dlply(melt.df, .(init), function(x) { - .jitter.ensmemb(sort(x$value), pan.width / 100) }), check.names=F), + .jitter.ensmemb(sort(x$value,na.last=T), pan.width / 100) }), check.names=F), value.name="yjitter", variable.name="init", id.vars=NULL) jitter.df$x <- melt(data.frame(dlply(melt.df, .(init), function(x) { - sort(x$value)}) ), value.name="x", id.vars=NULL)$x + sort(x$value,na.last=T)}) ), value.name="x", id.vars=NULL)$x } #------------------------ @@ -358,11 +358,12 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # include the point to the level. # Otherwise keep the point for another round. # Do one round in each direction to avoid uggly patterns. - if (is.unsorted(x)) { + if (is.unsorted(x,na.rm=T)) { stop("Provide a sorted array!") } lev <- x * 0 + lev[is.na(lev)] <- -1 level <- 1 while (any(lev == 0)) { last <- -1/0 @@ -388,6 +389,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, } level <- level + 1 } + lev[lev==-1] <- NA return(lev * thr * sqrt(3) / 2) } -- GitLab From 0f3103efb449315ac94f756aeaade58313da89b2 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 15 Mar 2019 11:28:54 +0100 Subject: [PATCH 17/21] Unittest, Imports added and .Rbuildignore fixed --- .Rbuildignore | 13 +++-- DESCRIPTION | 7 ++- NAMESPACE | 8 ++++ R/PlotForecastPDF.R | 28 ++++------- man/PlotForecastPDF.Rd | 47 +++++++++++++++++++ {tests => tets}/testthat.R | 0 .../testthat/test-AnoMultiMetric.R | 0 tets/testthat/test-PlotForecastPDF.R | 39 +++++++++++++++ 8 files changed, 115 insertions(+), 27 deletions(-) create mode 100644 man/PlotForecastPDF.Rd rename {tests => tets}/testthat.R (100%) rename {tests => tets}/testthat/test-AnoMultiMetric.R (100%) create mode 100644 tets/testthat/test-PlotForecastPDF.R diff --git a/.Rbuildignore b/.Rbuildignore index 4b988093..cdf20b53 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,7 +1,6 @@ -.git -.gitignore -.tar.gz -.pdf -./.nc -.data/ -.RData +.*\.git$ +.*\.gitignore$ +.*\.tar.gz$ +.*\.pdf$ +./.nc$ +.*\.RData$ diff --git a/DESCRIPTION b/DESCRIPTION index 362b4f8c..e787c9d3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,6 +10,7 @@ Authors@R: c( person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")), person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"), person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "ctb"), + person("Llorenc", "Lledo", , "llledo@bsc.es", role = "aut"), person("Other authors", role = "aut"), person("Other contributors", role = "ctb")) Description: Improved global climate model calibration and regionalization @@ -21,7 +22,11 @@ Depends: Imports: s2dverification, abind, - stats + stats, + ggplot2, + plyr, + data.table, + reshape2 Suggests: testthat License: LGPL-3 diff --git a/NAMESPACE b/NAMESPACE index 576f3a5d..03a3b388 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,5 +2,13 @@ export(AnoMultiMetric) export(MultivarRMSE) +export(PlotForecastPDF) +import(ggplot2) import(s2dverification) import(stats) +importFrom(data.table,CJ) +importFrom(data.table,data.table) +importFrom(data.table,setkey) +importFrom(plyr,.) +importFrom(plyr,dlply) +importFrom(reshape2,melt) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index cca7c916..ec157958 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -1,6 +1,6 @@ #'Plot one or multiple ensemble forecast pdfs for the same event #' -#'@author Llorenç Lledó \email{llledo@bsc.es} +#'@author Llorenc Lledo \email{llledo@bsc.es} #'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. #' #'@param fcst a dataframe or array containing all the ensember members for each frecast. If \code{"fcst"} is an array, it should have two labelled dimensions, and one of them should be \code{"members"}. If \code{"fcsts"} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. @@ -16,24 +16,27 @@ #' #'@return a ggplot object containing the plot. #' -#'@import data.table +#'@importFrom data.table data.table +#'@importFrom data.table CJ +#'@importFrom data.table setkey #'@import ggplot2 -#'@import reshape2 -#'@import plyr +#'@importFrom reshape2 melt +#'@importFrom plyr . +#'@importFrom plyr dlply #' #'@examples #'fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) #'PlotForecastPDF(fcsts,c(-1,1)) #'fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) #'PlotForecastPDF(fcsts2,c(-0.66,0.66),extreme.limits=c(-1.2,1.2),fcst.names=paste0("random fcst ",1:5),obs=0.7) -#' -#' #'@export PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", var.name="Varname (units)", fcst.names=NULL, add.ensmemb=c("above", "below", "no"), color.set=c("ggplot", "s2s4e", "hydro")) { + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- + y <- xend <- yend <- yjitter <- MLT <- NULL #------------------------ # Define color sets #------------------------ @@ -348,10 +351,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, return(plot) } -#================== -# A function to distribute ensemble members so not to overlap -# Requires a sorted array of values. -#================== .jitter.ensmemb <- function(x, thr=0.1) { # Idea: start with first level. Loop all points, # and if distance to last point in the level is more than a threshold, @@ -393,11 +392,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, return(lev * thr * sqrt(3) / 2) } - -#================== -# Hatch one polygon -# Based on base polygon function -#================== .polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven=F) { halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) cross <- halfplane[-1L] - halfplane[-length(halfplane)] @@ -426,10 +420,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, return(data.frame(x=lx1, y=ly1, xend=lx2, yend=ly2)) } -#================== -# Hatch one polygon -# Based on base polygon function -#================== .polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches=c(5, 1)) { x <- c(x, x[1L]) y <- c(y, y[1L]) diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd new file mode 100644 index 00000000..6d034108 --- /dev/null +++ b/man/PlotForecastPDF.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotForecastPDF.R +\name{PlotForecastPDF} +\alias{PlotForecastPDF} +\title{Plot one or multiple ensemble forecast pdfs for the same event} +\usage{ +PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, + plotfile = NULL, title = "Set a title", var.name = "Varname (units)", + fcst.names = NULL, add.ensmemb = c("above", "below", "no"), + color.set = c("ggplot", "s2s4e", "hydro")) +} +\arguments{ +\item{fcst}{a dataframe or array containing all the ensember members for each frecast. If \code{"fcst"} is an array, it should have two labelled dimensions, and one of them should be \code{"members"}. If \code{"fcsts"} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members.} + +\item{tercile.limits}{an array with P33 and P66 values that define the tercile categories.} + +\item{extreme.limits}{(optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown).} + +\item{obs}{(optional) a number with the observed value. (Default: observation is not shown).} + +\item{plotfile}{(optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved).} + +\item{title}{a string with the plot title.} + +\item{var.name}{a string with the variable name and units.} + +\item{fcst.names}{(optional) an array of strings with the titles of each individual forecast.} + +\item{add.ensmemb}{either to add the ensemble members (above or below the pdf) or not.} + +\item{color.set}{a selection of predefined color sets: use \code{"ggplot"} (default) for blue/green/red, \code{"s2s4e"} for blue/grey/orange, or \code{"hydro"} for yellow/gray/blue (suitable for precipitation and inflows).} +} +\value{ +a ggplot object containing the plot. +} +\description{ +This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. +} +\examples{ +fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) +PlotForecastPDF(fcsts,c(-1,1)) +fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) +PlotForecastPDF(fcsts2,c(-0.66,0.66),extreme.limits=c(-1.2,1.2),fcst.names=paste0("random fcst ",1:5),obs=0.7) +} +\author{ +Llorenc Lledo \email{llledo@bsc.es} +} diff --git a/tests/testthat.R b/tets/testthat.R similarity index 100% rename from tests/testthat.R rename to tets/testthat.R diff --git a/tests/testthat/test-AnoMultiMetric.R b/tets/testthat/test-AnoMultiMetric.R similarity index 100% rename from tests/testthat/test-AnoMultiMetric.R rename to tets/testthat/test-AnoMultiMetric.R diff --git a/tets/testthat/test-PlotForecastPDF.R b/tets/testthat/test-PlotForecastPDF.R new file mode 100644 index 00000000..c55c79e2 --- /dev/null +++ b/tets/testthat/test-PlotForecastPDF.R @@ -0,0 +1,39 @@ +context("Generic tests") +test_that("some checks", { + fcsts2 <- array(1:100, dim = c(members = 20, fcst = 5)) + expect_equal( + length(PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), + fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7)), + 10) + expect_equal( + names(PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), + fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7)), + c("data", "layers", "scales", "mapping", "theme", "coordinates", "facet", + "plot_env", "labels", "guides")) + expect_equal( + length(PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), + fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7)[[5]]), + 61) +}) + +test_that("Sanity checks", { + expect_error( + PlotForecastPDF(fcst, tercile.limits), + "object 'tercile.limits' not found") + expect_error( + PlotForecastPDF(fcst, tercile.limits = c(0.25, 0.55)), + "object 'fcst' not found") + expect_error( + PlotForecastPDF(fcst, tercile.limits = 10), + "Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") + expect_error( + PlotForecastPDF(fcst, tercile.limits = c(10, 20)), + "object 'fcst' not found") + fcsts2 <- array(rnorm(100),dim = c(members = 20, fcst = 5)) + expect_error( + PlotForecastPDF(fcst = fcsts2, tercile.limits), + "object 'tercile.limits' not found") + expect_error( + PlotForecastPDF(fcst = fcsts2, tercile.limits = c(-0.5, 0.5), extreme.limits = NA), + "Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") +}) \ No newline at end of file -- GitLab From 414edc1d6b379b7fedeed44c7f51d6391f9b9dae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Mon, 18 Mar 2019 14:16:28 +0100 Subject: [PATCH 18/21] Corrected many style changes and corrected documentation according to reviewer comments. Author name was corrected. --- R/PlotForecastPDF.R | 240 ++++++++++++++++++++++---------------------- 1 file changed, 120 insertions(+), 120 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index ec157958..12cacb96 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -1,7 +1,7 @@ #'Plot one or multiple ensemble forecast pdfs for the same event #' -#'@author Llorenc Lledo \email{llledo@bsc.es} -#'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. +#'@author Llorenç Lledó \email{llledo@bsc.es} +#'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. #' #'@param fcst a dataframe or array containing all the ensember members for each frecast. If \code{"fcst"} is an array, it should have two labelled dimensions, and one of them should be \code{"members"}. If \code{"fcsts"} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. #'@param tercile.limits an array with P33 and P66 values that define the tercile categories. @@ -11,7 +11,7 @@ #'@param title a string with the plot title. #'@param var.name a string with the variable name and units. #'@param fcst.names (optional) an array of strings with the titles of each individual forecast. -#'@param add.ensmemb either to add the ensemble members (above or below the pdf) or not. +#'@param add.ensmemb either to add the ensemble members \code{"above"} (default) or \code{"below"} the pdf, or not (\code{"no"}). #'@param color.set a selection of predefined color sets: use \code{"ggplot"} (default) for blue/green/red, \code{"s2s4e"} for blue/grey/orange, or \code{"hydro"} for yellow/gray/blue (suitable for precipitation and inflows). #' #'@return a ggplot object containing the plot. @@ -30,12 +30,12 @@ #'fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) #'PlotForecastPDF(fcsts2,c(-0.66,0.66),extreme.limits=c(-1.2,1.2),fcst.names=paste0("random fcst ",1:5),obs=0.7) #'@export -PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, +PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", - var.name="Varname (units)", fcst.names=NULL, - add.ensmemb=c("above", "below", "no"), + var.name="Varname (units)", fcst.names=NULL, + add.ensmemb=c("above", "below", "no"), color.set=c("ggplot", "s2s4e", "hydro")) { - value <- init <- extremes <- x <- ymin <- ymax <- tercile <- + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- y <- xend <- yend <- yjitter <- MLT <- NULL #------------------------ # Define color sets @@ -46,24 +46,24 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, colorHatch <- c("deepskyblue3", "indianred3") colorMember <- c("#ffff7f") colorObs <- "purple" - colorLab <- c("blue", "red") + colorLab <- c("blue", "red") } else if (color.set == "hydro") { colorFill <- rev(c("#41CBC9", "#b5b5b5", "#FFAB38")) colorHatch <- c("darkorange1", "deepskyblue3") colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("darkorange3", "blue") + colorObs <- "purple" + colorLab <- c("darkorange3", "blue") } else if (color.set == "ggplot") { - gg.color.hue <- function(n) { - hues=seq(15, 375, length=n+1) - hcl(h=hues, l=65, c=100)[1:n] + ggColorHue <- function(n) { + hues <- seq(15, 375, length = n + 1) + hcl(h = hues, l = 65, c = 100)[1:n] } - colorFill <- rev(gg.color.hue(3)) + colorFill <- rev(ggColorHue(3)) colorHatch <- c("deepskyblue3", "indianred1") colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("blue", "red") - } else { + colorObs <- "purple" + colorLab <- c("blue", "red") + } else { stop("Parameter 'color.set' should be one of ggplot/s2s4e/hydro") } @@ -81,7 +81,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, if (length(extreme.limits) != 2) { stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") } - if (extreme.limits[1] >= tercile.limits[1] | + if (extreme.limits[1] >= tercile.limits[1] | extreme.limits[2] <= tercile.limits[2]) { stop("The provided extreme limits are not consistent with tercile limits") } @@ -90,12 +90,12 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, #------------------------ # Check fcst type and convert to data.frame if needed #------------------------ - if (is.array(fcst)) { + if (is.array(fcst)) { if (!"members" %in% names(dim(fcst)) | length(dim(fcst)) != 2) { stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be 'members'") } dim.members <- which(names(dim(fcst)) == "members") - if (dim.members == 1) { + if (dim.members == 1) { fcst.df <- data.frame(fcst) } else { fcst.df <- data.frame(t(fcst)) @@ -110,21 +110,21 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Set proper fcst names #------------------------ if (!is.null(fcst.names)) { - if (length(fcst.names) != dim(fcst.df)[2]) { + if (length(fcst.names) != dim(fcst.df)[2]) { stop("Parameter 'fcst.names' should be an array with as many names as distinct forecasts provided") - } - colnames(fcst.df) <- factor(fcst.names, levels=fcst.names) + } + colnames(fcst.df) <- factor(fcst.names, levels = fcst.names) } #------------------------ # Produce a first plot with the pdf for each init in a panel #------------------------ - melt.df <- melt(fcst.df, variable.name="init", id.vars=NULL) - plot <- ggplot(melt.df, aes(x=value)) + - geom_density(alpha=1, na.rm=T) + - coord_flip() + - facet_wrap(~init, strip.position="top", nrow=1) + - xlim(range(c(obs, density(melt.df$value, na.rm=T)$x))) + melt.df <- melt(fcst.df, variable.name = "init", id.vars = NULL) + plot <- ggplot(melt.df, aes(x = value)) + + geom_density(alpha = 1, na.rm = T) + + coord_flip() + + facet_wrap(~init, strip.position = "top", nrow = 1) + + xlim(range(c(obs, density(melt.df$value, na.rm = T)$x))) ggp <- ggplot_build(plot) #------------------------ @@ -139,9 +139,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, } else { stop("Cannot find PANELS in ggp object") } - tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[1], "Below normal", - ifelse(tmp.df$x < tercile.limits[2], "Normal", "Above normal")), - levels=c("Below normal", "Normal", "Above normal")) + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[1], "Below normal", + ifelse(tmp.df$x < tercile.limits[2], "Normal", "Above normal")), + levels = c("Below normal", "Normal", "Above normal")) #------------------------ # Get the height and width of a panel @@ -154,35 +154,35 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, # Compute hatch coordinates for extremes #------------------------ if (!is.null(extreme.limits)) { - tmp.df$extremes <- factor(ifelse(tmp.df$x thr) { @@ -388,7 +388,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, } level <- level + 1 } - lev[lev==-1] <- NA + lev[lev == -1] <- NA return(lev * thr * sqrt(3) / 2) } @@ -403,12 +403,12 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, y1 <- y[-length(y)][does.cross] x2 <- x[-1L][does.cross] y2 <- y[-1L][does.cross] - t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * (x2 - x1))) + t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1)) / (xd * (y2 - y1) - yd * (x2 - x1))) o <- order(t) tsort <- t[o] crossings <- cumsum(cross[does.cross][o]) if (fillOddEven) { - crossings <- crossings%%2 + crossings <- crossings %% 2 } drawline <- crossings != 0 lx <- x0 + xd * tsort @@ -417,13 +417,13 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, ly1 <- ly[-length(ly)][drawline] lx2 <- lx[-1L][drawline] ly2 <- ly[-1L][drawline] - return(data.frame(x=lx1, y=ly1, xend=lx2, yend=ly2)) + return(data.frame(x = lx1, y = ly1, xend = lx2, yend = ly2)) } .polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches=c(5, 1)) { x <- c(x, x[1L]) y <- c(y, y[1L]) - angle <- angle%%180 + angle <- angle %% 180 upi <- c(width.units, height.units) / inches if (upi[1L] < 0) { angle <- 180 - angle @@ -446,8 +446,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, } y.shift <- upi[2L] / density / abs(cos(angle / 180 * pi)) x0 <- 0 - y0 <- floor((min(y) - first.x * yd/xd) / y.shift) * y.shift - y.end <- max(y) - last.x * yd/xd + y0 <- floor((min(y) - first.x * yd / xd) / y.shift) * y.shift + y.end <- max(y) - last.x * yd / xd while (y0 < y.end) { hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) i <- i + 1 -- GitLab From 72a17b55a90f0f65b695dc54812d918c474cca7b Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 19 Mar 2019 17:31:25 +0100 Subject: [PATCH 19/21] Merging branches correction --- .Rbuildignore | 1 + DESCRIPTION | 1 - NAMESPACE | 5 ++--- R/PlotForecastPDF.R | 8 +++++--- man/PlotForecastPDF.Rd | 14 ++++++++------ tests/testthat.R | 5 +++++ 6 files changed, 21 insertions(+), 13 deletions(-) create mode 100644 tests/testthat.R diff --git a/.Rbuildignore b/.Rbuildignore index cdf20b53..43690b70 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,4 @@ .*\.pdf$ ./.nc$ .*\.RData$ +.*\.gitlab-ci.yml$ diff --git a/DESCRIPTION b/DESCRIPTION index 5c834d86..493b057a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,7 +14,6 @@ Authors@R: c( person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")), person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"), person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut")) - Description: Improved global climate model calibration and regionalization techniques as well as better forecast verification methods need to be developed for Mediterranean region to extract as much climate information as possible from diff --git a/NAMESPACE b/NAMESPACE index 5a0ac42f..f74095a4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,10 +1,9 @@ # Generated by roxygen2: do not edit by hand - -export(PlotForecastPDF) -import(ggplot2) export(CST_MultiMetric) export(CST_MultivarRMSE) +export(PlotForecastPDF) +import(ggplot2) import(s2dverification) import(stats) importFrom(data.table,CJ) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 12cacb96..e4dc6cfa 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -25,10 +25,12 @@ #'@importFrom plyr dlply #' #'@examples -#'fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) +#'fcsts <- data.frame(fcst1 = rnorm(50), fcst2 = rnorm(50, 0.5, 1.2), +#' fcst3 = rnorm(50, -0.5, 0.9)) #'PlotForecastPDF(fcsts,c(-1,1)) -#'fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) -#'PlotForecastPDF(fcsts2,c(-0.66,0.66),extreme.limits=c(-1.2,1.2),fcst.names=paste0("random fcst ",1:5),obs=0.7) +#'fcsts2 <- array(rnorm(100), dim = c(members = 20, fcst = 5)) +#'PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), +#' fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7) #'@export PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, obs=NULL, plotfile=NULL, title="Set a title", diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd index 6d034108..89fecb49 100644 --- a/man/PlotForecastPDF.Rd +++ b/man/PlotForecastPDF.Rd @@ -26,7 +26,7 @@ PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, \item{fcst.names}{(optional) an array of strings with the titles of each individual forecast.} -\item{add.ensmemb}{either to add the ensemble members (above or below the pdf) or not.} +\item{add.ensmemb}{either to add the ensemble members \code{"above"} (default) or \code{"below"} the pdf, or not (\code{"no"}).} \item{color.set}{a selection of predefined color sets: use \code{"ggplot"} (default) for blue/green/red, \code{"s2s4e"} for blue/grey/orange, or \code{"hydro"} for yellow/gray/blue (suitable for precipitation and inflows).} } @@ -34,14 +34,16 @@ PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, a ggplot object containing the plot. } \description{ -This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. +This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. } \examples{ -fcsts <- data.frame(fcst1=rnorm(50),fcst2=rnorm(50,0.5,1.2),fcst3=rnorm(50,-0.5,0.9)) +fcsts <- data.frame(fcst1 = rnorm(50), fcst2 = rnorm(50, 0.5, 1.2), + fcst3 = rnorm(50, -0.5, 0.9)) PlotForecastPDF(fcsts,c(-1,1)) -fcsts2 <- array(rnorm(100),dim=c(members=20,fcst=5)) -PlotForecastPDF(fcsts2,c(-0.66,0.66),extreme.limits=c(-1.2,1.2),fcst.names=paste0("random fcst ",1:5),obs=0.7) +fcsts2 <- array(rnorm(100), dim = c(members = 20, fcst = 5)) +PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), + fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7) } \author{ -Llorenc Lledo \email{llledo@bsc.es} +Llorenç Lledó \email{llledo@bsc.es} } diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 00000000..c069d753 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,5 @@ +library(testthat) +library(CSTools) + +test_check("CSTools") + -- GitLab From d09fb7103058d154134fbda789cc909e4f23a112 Mon Sep 17 00:00:00 2001 From: nperez Date: Tue, 19 Mar 2019 18:10:46 +0100 Subject: [PATCH 20/21] correction of automatic check failed Format failue --- R/PlotForecastPDF.R | 812 +++++++++++++++++++++----------------------- 1 file changed, 386 insertions(+), 426 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index e4dc6cfa..83b4e239 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -3,7 +3,7 @@ #'@author Llorenç Lledó \email{llledo@bsc.es} #'@description This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond. #' -#'@param fcst a dataframe or array containing all the ensember members for each frecast. If \code{"fcst"} is an array, it should have two labelled dimensions, and one of them should be \code{"members"}. If \code{"fcsts"} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. +#'@param fcst a dataframe or array containing all the ensember members for each frecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members. #'@param tercile.limits an array with P33 and P66 values that define the tercile categories. #'@param extreme.limits (optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown). #'@param obs (optional) a number with the observed value. (Default: observation is not shown). @@ -11,8 +11,8 @@ #'@param title a string with the plot title. #'@param var.name a string with the variable name and units. #'@param fcst.names (optional) an array of strings with the titles of each individual forecast. -#'@param add.ensmemb either to add the ensemble members \code{"above"} (default) or \code{"below"} the pdf, or not (\code{"no"}). -#'@param color.set a selection of predefined color sets: use \code{"ggplot"} (default) for blue/green/red, \code{"s2s4e"} for blue/grey/orange, or \code{"hydro"} for yellow/gray/blue (suitable for precipitation and inflows). +#'@param add.ensmemb either to add the ensemble members \code{'above'} (default) or \code{'below'} the pdf, or not (\code{'no'}). +#'@param color.set a selection of predefined color sets: use \code{'ggplot'} (default) for blue/green/red, \code{'s2s4e'} for blue/grey/orange, or \code{'hydro'} for yellow/gray/blue (suitable for precipitation and inflows). #' #'@return a ggplot object containing the plot. #' @@ -30,448 +30,408 @@ #'PlotForecastPDF(fcsts,c(-1,1)) #'fcsts2 <- array(rnorm(100), dim = c(members = 20, fcst = 5)) #'PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), -#' fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7) +#' fcst.names = paste0('random fcst ', 1 : 5), obs = 0.7) #'@export -PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits=NULL, - obs=NULL, plotfile=NULL, title="Set a title", - var.name="Varname (units)", fcst.names=NULL, - add.ensmemb=c("above", "below", "no"), - color.set=c("ggplot", "s2s4e", "hydro")) { - value <- init <- extremes <- x <- ymin <- ymax <- tercile <- - y <- xend <- yend <- yjitter <- MLT <- NULL - #------------------------ - # Define color sets - #------------------------ - color.set <- match.arg(color.set) - if (color.set == "s2s4e") { - colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1")) - colorHatch <- c("deepskyblue3", "indianred3") - colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("blue", "red") - } else if (color.set == "hydro") { - colorFill <- rev(c("#41CBC9", "#b5b5b5", "#FFAB38")) - colorHatch <- c("darkorange1", "deepskyblue3") - colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("darkorange3", "blue") - } else if (color.set == "ggplot") { - ggColorHue <- function(n) { - hues <- seq(15, 375, length = n + 1) - hcl(h = hues, l = 65, c = 100)[1:n] +PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, + plotfile = NULL, title = "Set a title", var.name = "Varname (units)", fcst.names = NULL, + add.ensmemb = c("above", "below", "no"), color.set = c("ggplot", "s2s4e", "hydro")) { + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- y <- xend <- yend <- yjitter <- MLT <- NULL + #------------------------ + # Define color sets + #------------------------ + color.set <- match.arg(color.set) + if (color.set == "s2s4e") { + colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1")) + colorHatch <- c("deepskyblue3", "indianred3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue", "red") + } else if (color.set == "hydro") { + colorFill <- rev(c("#41CBC9", "#b5b5b5", "#FFAB38")) + colorHatch <- c("darkorange1", "deepskyblue3") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("darkorange3", "blue") + } else if (color.set == "ggplot") { + ggColorHue <- function(n) { + hues <- seq(15, 375, length = n + 1) + hcl(h = hues, l = 65, c = 100)[1:n] + } + colorFill <- rev(ggColorHue(3)) + colorHatch <- c("deepskyblue3", "indianred1") + colorMember <- c("#ffff7f") + colorObs <- "purple" + colorLab <- c("blue", "red") + } else { + stop("Parameter 'color.set' should be one of ggplot/s2s4e/hydro") + } + #------------------------ + # Check input arguments + #------------------------ + add.ensmemb <- match.arg(add.ensmemb) + if (length(tercile.limits) != 2) { + stop("Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") } - colorFill <- rev(ggColorHue(3)) - colorHatch <- c("deepskyblue3", "indianred1") - colorMember <- c("#ffff7f") - colorObs <- "purple" - colorLab <- c("blue", "red") - } else { - stop("Parameter 'color.set' should be one of ggplot/s2s4e/hydro") - } - - #------------------------ - # Check input arguments - #------------------------ - add.ensmemb <- match.arg(add.ensmemb) - if (length(tercile.limits) != 2) { - stop("Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") - } - if (tercile.limits[1] >= tercile.limits[2]) { - stop("The provided tercile limits are in the wrong order") - } - if (!is.null(extreme.limits)) { - if (length(extreme.limits) != 2) { - stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") + if (tercile.limits[1] >= tercile.limits[2]) { + stop("The provided tercile limits are in the wrong order") } - if (extreme.limits[1] >= tercile.limits[1] | - extreme.limits[2] <= tercile.limits[2]) { - stop("The provided extreme limits are not consistent with tercile limits") + if (!is.null(extreme.limits)) { + if (length(extreme.limits) != 2) { + stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") + } + if (extreme.limits[1] >= tercile.limits[1] | extreme.limits[2] <= tercile.limits[2]) { + stop("The provided extreme limits are not consistent with tercile limits") + } + } + #------------------------ + # Check fcst type and convert to data.frame if needed + #------------------------ + if (is.array(fcst)) { + if (!"members" %in% names(dim(fcst)) | length(dim(fcst)) != 2) { + stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be 'members'") + } + dim.members <- which(names(dim(fcst)) == "members") + if (dim.members == 1) { + fcst.df <- data.frame(fcst) + } else { + fcst.df <- data.frame(t(fcst)) + } + } else if (is.data.frame(fcst)) { + fcst.df <- fcst + } else { + stop("Parameter 'fcst' should be an array or a data.frame") } - } - - #------------------------ - # Check fcst type and convert to data.frame if needed - #------------------------ - if (is.array(fcst)) { - if (!"members" %in% names(dim(fcst)) | length(dim(fcst)) != 2) { - stop("Parameter 'fcst' should be a two-dimensional array with labelled dimensions and one of them should be 'members'") + #------------------------ + # Set proper fcst names + #------------------------ + if (!is.null(fcst.names)) { + if (length(fcst.names) != dim(fcst.df)[2]) { + stop("Parameter 'fcst.names' should be an array with as many names as distinct forecasts provided") + } + colnames(fcst.df) <- factor(fcst.names, levels = fcst.names) } - dim.members <- which(names(dim(fcst)) == "members") - if (dim.members == 1) { - fcst.df <- data.frame(fcst) + #------------------------ + # Produce a first plot with the pdf for each init in a panel + #------------------------ + melt.df <- melt(fcst.df, variable.name = "init", id.vars = NULL) + plot <- ggplot(melt.df, aes(x = value)) + geom_density(alpha = 1, na.rm = T) + + coord_flip() + facet_wrap(~init, strip.position = "top", nrow = 1) + xlim(range(c(obs, + density(melt.df$value, na.rm = T)$x))) + ggp <- ggplot_build(plot) + #------------------------ + # Gather the coordinates of the plots together with init and corresponding + # terciles + #------------------------ + tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] + if (!is.null(ggp$layout$layout)) { + tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] + } else if (!is.null(ggp$layout$panel_layout)) { + tmp.df$init <- ggp$layout$panel_layout$init[as.numeric(tmp.df$PANEL)] } else { - fcst.df <- data.frame(t(fcst)) + stop("Cannot find PANELS in ggp object") + } + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[1], "Below normal", + ifelse(tmp.df$x < tercile.limits[2], "Normal", "Above normal")), levels = c("Below normal", + "Normal", "Above normal")) + #------------------------ + # Get the height and width of a panel + #------------------------ + pan.width <- diff(range(tmp.df$x)) + pan.height <- max(tmp.df$ymax) + magic.ratio <- 9 * pan.height/pan.width + #------------------------ + # Compute hatch coordinates for extremes + #------------------------ + if (!is.null(extreme.limits)) { + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[1], "Below P10", + ifelse(tmp.df$x < extreme.limits[2], "Normal", "Above P90")), levels = c("Below P10", + "Normal", "Above P90")) + hatch.ls <- dlply(tmp.df, .(init, extremes), function(x) { + # close the polygon + tmp.df2 <- data.frame(x = c(x$x, max(x$x), min(x$x)), y = c(x$ymax, 0, + 0)) + # compute the hatches for this polygon + hatches <- .polygon.fullhatch(tmp.df2$x, tmp.df2$y, angle = 60, density = 10, + width.units = pan.width, height.units = pan.height) + # add bottom segment + end1 <- data.frame(x = x$x[1], y = x$ymax[1], xend = x$x[1], yend = 0) + # add top segment + end2 <- data.frame(x = x$x[length(x$x)], y = x$ymax[length(x$x)], xend = x$x[length(x$x)], + yend = 0) + return(rbind(hatches, end1, end2)) + }) + attr <- attr(hatch.ls, "split_labels") + for (i in 1:length(hatch.ls)) { + hatch.ls[[i]] <- cbind(hatch.ls[[i]], attr[i, ]) + } + hatch.df <- do.call("rbind", hatch.ls) + # Compute max y for each extreme category + max.ls <- dlply(tmp.df, .(init, extremes), function(x) { + data.frame(y = min(0.85 * pan.height, max(x$ymax))) + }) + attr <- attr(max.ls, "split_labels") + for (i in 1:length(max.ls)) { + max.ls[[i]] <- cbind(max.ls[[i]], attr[i, ]) + } + max.df <- do.call("rbind", max.ls) + } + #------------------------ + # Compute jitter space for ensemble members + #------------------------ + if (add.ensmemb != "no") { + jitter.df <- melt(data.frame(dlply(melt.df, .(init), function(x) { + .jitter.ensmemb(sort(x$value, na.last = T), pan.width/100) + }), check.names = F), value.name = "yjitter", variable.name = "init", id.vars = NULL) + jitter.df$x <- melt(data.frame(dlply(melt.df, .(init), function(x) { + sort(x$value, na.last = T) + })), value.name = "x", id.vars = NULL)$x + } + #------------------------ + # Get y coordinates for observed x values, using a cool data.table feature: merge + # to nearest value + #------------------------ + if (!is.null(obs)) { + tmp.dt <- data.table(tmp.df, key = c("init", "x")) + obs.dt <- data.table(init = factor(colnames(fcst.df), levels = colnames(fcst.df)), + value = rep(obs, dim(fcst.df)[2])) + setkey(obs.dt, init, value) + obs.xy <- tmp.dt[obs.dt, roll = "nearest"] + } + #------------------------ + # Fill each pdf with different colors for the terciles + #------------------------ + plot <- plot + geom_ribbon(data = tmp.df, aes(x = x, ymin = ymin, ymax = ymax, + fill = tercile), alpha = 0.7) + #------------------------ + # Add hatches for extremes + #------------------------ + if (!is.null(extreme.limits)) { + if (nrow(hatch.df[hatch.df$extremes != "Normal", ]) == 0) { + warning("The provided extreme categories are outside the plot bounds. The extremes will not be drawn.") + extreme.limits <- NULL + } else { + plot <- plot + geom_segment(data = hatch.df[hatch.df$extremes != "Normal", + ], aes(x = x, y = y, xend = xend, yend = yend, color = extremes)) + } } - } else if (is.data.frame(fcst)) { - fcst.df <- fcst - } else { - stop("Parameter 'fcst' should be an array or a data.frame") - } - - #------------------------ - # Set proper fcst names - #------------------------ - if (!is.null(fcst.names)) { - if (length(fcst.names) != dim(fcst.df)[2]) { - stop("Parameter 'fcst.names' should be an array with as many names as distinct forecasts provided") + #------------------------ + # Add obs line + #------------------------ + if (!is.null(obs)) { + plot <- plot + geom_vline(data = obs.dt, aes(xintercept = value), linetype = "dashed", + color = colorObs) } - colnames(fcst.df) <- factor(fcst.names, levels = fcst.names) - } - - #------------------------ - # Produce a first plot with the pdf for each init in a panel - #------------------------ - melt.df <- melt(fcst.df, variable.name = "init", id.vars = NULL) - plot <- ggplot(melt.df, aes(x = value)) + - geom_density(alpha = 1, na.rm = T) + - coord_flip() + - facet_wrap(~init, strip.position = "top", nrow = 1) + - xlim(range(c(obs, density(melt.df$value, na.rm = T)$x))) - ggp <- ggplot_build(plot) - - #------------------------ - # Gather the coordinates of the plots - # together with init and corresponding terciles - #------------------------ - tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] - if (!is.null(ggp$layout$layout)) { - tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] - } else if (!is.null(ggp$layout$panel_layout)) { - tmp.df$init <- ggp$layout$panel_layout$init[as.numeric(tmp.df$PANEL)] - } else { - stop("Cannot find PANELS in ggp object") - } - tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[1], "Below normal", - ifelse(tmp.df$x < tercile.limits[2], "Normal", "Above normal")), - levels = c("Below normal", "Normal", "Above normal")) - - #------------------------ - # Get the height and width of a panel - #------------------------ - pan.width <- diff(range(tmp.df$x)) - pan.height <- max(tmp.df$ymax) - magic.ratio <- 9 * pan.height / pan.width - - #------------------------ - # Compute hatch coordinates for extremes - #------------------------ - if (!is.null(extreme.limits)) { - tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[1], "Below P10", - ifelse(tmp.df$x < extreme.limits[2], "Normal", "Above P90")), - levels = c("Below P10", "Normal", "Above P90")) - hatch.ls <- dlply(tmp.df, .(init, extremes), function(x) { - # close the polygon - tmp.df2 <- data.frame(x = c(x$x, max(x$x), min(x$x)), y = c(x$ymax, 0, 0)) - # compute the hatches for this polygon - hatches <- .polygon.fullhatch(tmp.df2$x, tmp.df2$y, angle = 60, density = 10, - width.units = pan.width, height.units = pan.height) - # add bottom segment - end1 <- data.frame(x = x$x[1], y = x$ymax[1], xend = x$x[1], yend = 0) - # add top segment - end2 <- data.frame(x = x$x[length(x$x)], y = x$ymax[length(x$x)], - xend = x$x[length(x$x)], yend = 0) - return(rbind(hatches, end1, end2)) - }) - attr <- attr(hatch.ls, "split_labels") - for (i in 1:length(hatch.ls)) { - hatch.ls[[i]] <- cbind(hatch.ls[[i]], attr[i, ]) + #------------------------ + # Add ensemble members + #------------------------ + if (add.ensmemb == "below") { + plot <- plot + # this adds a grey box for ensmembers + geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = -pan.height/10), + fill = "gray95", color = "black", width = 0.2) + # this adds the ensemble members + geom_point(data = jitter.df, color = "black", fill = colorMember, alpha = 1, + aes(x = x, y = -pan.height/10 - magic.ratio * yjitter, shape = "Ensemble members")) + } else if (add.ensmemb == "above") { + plot <- plot + geom_point(data = jitter.df, color = "black", fill = colorMember, + alpha = 1, aes(x = x, y = 0.7 * magic.ratio * yjitter, shape = "Ensemble members")) } - hatch.df <- do.call("rbind", hatch.ls) - - # Compute max y for each extreme category - max.ls <- dlply(tmp.df, .(init, extremes), function(x) { - data.frame(y = min(0.85 * pan.height, max(x$ymax))) - }) - attr <- attr(max.ls, "split_labels") - for (i in 1:length(max.ls)) { - max.ls[[i]] <- cbind(max.ls[[i]], attr[i, ]) + #------------------------ + # Add obs diamond + #------------------------ + if (!is.null(obs)) { + plot <- plot + # this adds the obs diamond + geom_point(data = obs.xy, aes(x = x, y = ymax, size = "Observation"), shape = 23, + color = "black", fill = colorObs) } - max.df <- do.call("rbind", max.ls) - } - - #------------------------ - # Compute jitter space for ensemble members - #------------------------ - if (add.ensmemb != "no") { - jitter.df <- melt(data.frame(dlply(melt.df, .(init), function(x) { - .jitter.ensmemb(sort(x$value, na.last = T), pan.width / 100) }), check.names = F), - value.name = "yjitter", variable.name = "init", id.vars = NULL) - jitter.df$x <- melt(data.frame(dlply(melt.df, .(init), function(x) { - sort(x$value, na.last = T)}) ), value.name = "x", id.vars = NULL)$x - } - - #------------------------ - # Get y coordinates for observed x values, - # using a cool data.table feature: merge to nearest value - #------------------------ - if (!is.null(obs)) { - tmp.dt <- data.table(tmp.df, key = c("init", "x")) - obs.dt <- data.table(init = factor(colnames(fcst.df), levels = colnames(fcst.df)), - value = rep(obs, dim(fcst.df)[2])) - setkey(obs.dt, init, value) - obs.xy <- tmp.dt[obs.dt, roll = "nearest"] - } - - #------------------------ - # Fill each pdf with different colors for the terciles - #------------------------ - plot <- plot + - geom_ribbon(data = tmp.df, aes(x = x, ymin = ymin, ymax = ymax, fill = tercile), alpha = 0.7) - - #------------------------ - # Add hatches for extremes - #------------------------ - if (!is.null(extreme.limits)) { - if (nrow(hatch.df[hatch.df$extremes != "Normal", ]) == 0) { - warning("The provided extreme categories are outside the plot bounds. The extremes will not be drawn.") - extreme.limits <- NULL + #------------------------ + # Compute probability for each tercile and identify MLT + #------------------------ + tmp.dt <- data.table(tmp.df) + pct <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), + by = .(init, tercile)] + tot <- pct[, .(tot = sum(pct)), by = init] + pct <- merge(pct, tot, by = "init") + pct$pct <- round(100 * pct$pct/pct$tot, 0) + pct$MLT <- pct[, .(MLT = pct == max(pct)), by = init]$MLT + lab.pos <- c(tercile.limits[1], mean(tercile.limits), tercile.limits[2]) + #------------------------ + # Compute probability for extremes + #------------------------ + if (!is.null(extreme.limits)) { + pct2 <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), + by = .(init, extremes)] + tot2 <- pct2[, .(tot = sum(pct)), by = init] + pct2 <- merge(pct2, tot2, by = "init") + pct2$pct <- round(100 * pct2$pct/pct2$tot, 0) + lab.pos2 <- c(extreme.limits[1], NA, extreme.limits[2]) + pct2 <- merge(pct2, max.df, by = c("init", "extremes")) + # include potentially missing groups + pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), + levels = c("Below P10", "Normal", "Above P90"))), ] + } + #------------------------ + # Add probability labels for terciles + #------------------------ + if (add.ensmemb == "above") { + labpos <- -0.2 * pan.height + vjust <- 0 } else { - plot <- plot + - geom_segment(data = hatch.df[hatch.df$extremes != "Normal", ], aes(x = x, y = y, - xend = xend, yend = yend, color = extremes)) + labpos <- 0 + vjust <- -0.5 + } + plot <- plot + geom_text(data = pct, aes(x = lab.pos[as.integer(tercile)], y = labpos, + label = paste0(pct, "%"), hjust = as.integer(tercile) * 1.5 - 2.5), vjust = vjust, + angle = -90, size = 3.2) + geom_text(data = pct[MLT == T, ], aes(x = lab.pos[as.integer(tercile)], + y = labpos, label = "*", hjust = as.integer(tercile) * 3.5 - 5), vjust = 0.1, + angle = -90, size = 7, color = "black") + #------------------------ + # Add probability labels for extremes + #------------------------ + if (!is.null(extreme.limits)) { + plot <- plot + geom_text(data = pct2[extremes != "Normal", ], aes(x = lab.pos2[as.integer(extremes)], + y = 0.9 * y, label = paste0(pct, "%"), hjust = as.integer(extremes) * + 1.5 - 2.5), vjust = -0.5, angle = -90, size = 3.2, color = rep(colorLab, + dim(fcst.df)[2])) } - } - - #------------------------ - # Add obs line - #------------------------ - if (!is.null(obs)) { - plot <- plot + - geom_vline(data = obs.dt, aes(xintercept = value), linetype = "dashed", color = colorObs) - } - - #------------------------ - # Add ensemble members - #------------------------ - if (add.ensmemb == "below") { - plot <- plot + - # this adds a grey box for ensmembers - geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = -pan.height / 10), - fill = "gray95", color = "black", width = 0.2) + - # this adds the ensemble members - geom_point(data = jitter.df, color = "black", fill = colorMember, alpha = 1, - aes(x = x, y = -pan.height / 10 - magic.ratio * yjitter, shape = "Ensemble members")) - } else if (add.ensmemb == "above") { - plot <- plot + - geom_point(data = jitter.df, color = "black", fill = colorMember, alpha = 1, - aes(x = x, y = 0.7 * magic.ratio * yjitter, shape = "Ensemble members")) - } - - #------------------------ - # Add obs diamond - #------------------------ - if (!is.null(obs)) { - plot <- plot + - # this adds the obs diamond - geom_point(data = obs.xy, aes(x = x, y = ymax, size = "Observation"), - shape = 23, color = "black", fill = colorObs) - } - - #------------------------ - # Compute probability for each tercile and identify MLT - #------------------------ - tmp.dt <- data.table(tmp.df) - pct <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), - by = .(init, tercile)] - tot <- pct[, .(tot = sum(pct)), by = init] - pct <- merge(pct, tot, by = "init") - pct$pct <- round(100 * pct$pct / pct$tot, 0) - pct$MLT <- pct[, .(MLT = pct == max(pct)), by = init]$MLT - lab.pos <- c(tercile.limits[1], mean(tercile.limits), tercile.limits[2]) - - #------------------------ - # Compute probability for extremes - #------------------------ - if (!is.null(extreme.limits)) { - pct2 <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), - by = .(init, extremes)] - tot2 <- pct2[, .(tot = sum(pct)), by = init] - pct2 <- merge(pct2, tot2, by = "init") - pct2$pct <- round(100 * pct2$pct / pct2$tot, 0) - lab.pos2 <- c(extreme.limits[1], NA, extreme.limits[2]) - pct2 <- merge(pct2, max.df, by = c("init", "extremes")) - # include potentially missing groups - pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), - levels = c("Below P10", "Normal", "Above P90"))), ] - } - - #------------------------ - # Add probability labels for terciles - #------------------------ - if (add.ensmemb == "above") { - labpos <- -0.2 * pan.height - vjust <- 0 - } else { - labpos <- 0 - vjust <- -0.5 - } - plot <- plot + - geom_text(data = pct, aes(x = lab.pos[as.integer(tercile)], y = labpos, label = paste0(pct, "%"), - hjust = as.integer(tercile) * 1.5 - 2.5), vjust = vjust, angle = -90, size = 3.2) + - geom_text(data = pct[MLT == T, ], aes(x = lab.pos[as.integer(tercile)], y = labpos, - label = "*", hjust = as.integer(tercile) * 3.5 - 5.0), vjust = 0.1, angle = -90, - size = 7, color = "black") - - #------------------------ - # Add probability labels for extremes - #------------------------ - if (!is.null(extreme.limits)) { - plot <- plot + - geom_text(data = pct2[extremes != "Normal", ], aes(x = lab.pos2[as.integer(extremes)], - y = 0.9 * y, label = paste0(pct, "%"), hjust = as.integer(extremes) * 1.5 - 2.5), - vjust = -0.5, angle = -90, size = 3.2, color = rep(colorLab, dim(fcst.df)[2])) - } - - #------------------------ - # Finish all theme and legend details - #------------------------ - plot <- plot + - theme_minimal() + - scale_fill_manual(name = "Probability of\nterciles", - breaks = c("Above normal", "Normal", "Below normal"), - values = colorFill, drop = F) + - scale_color_manual(name = "Probability of\nextremes", values = colorHatch) + - scale_shape_manual(name = "Ensemble\nmembers", values = c(21)) + - scale_size_manual(name = "Observation", values = c(3)) + - labs(x = var.name, y = "Probability density\n(total area=1)", title = title) + - theme(axis.text.x = element_blank(), panel.grid.minor.x = element_blank(), - legend.key.size = unit(0.3, "in"), - panel.border = element_rect(fill = NA, color = "black"), - strip.background = element_rect(colour = "black", fill = "gray80"), - panel.spacing = unit(0.2, "in"), - panel.grid.major.x = element_line(color = "grey93")) + - guides(fill = guide_legend(order = 1), color = guide_legend(order = 2, reverse = T), - shape = guide_legend(order = 3, label = F), size = guide_legend(order = 4, label = F)) - - #------------------------ - # Save to plotfile if needed, and return plot - #------------------------ - if (!is.null(plotfile)) { - ggsave(plotfile, plot) - } - - return(plot) + #------------------------ + # Finish all theme and legend details + #------------------------ + plot <- plot + theme_minimal() + scale_fill_manual(name = "Probability of\nterciles", + breaks = c("Above normal", "Normal", "Below normal"), values = colorFill, + drop = F) + scale_color_manual(name = "Probability of\nextremes", values = colorHatch) + + scale_shape_manual(name = "Ensemble\nmembers", values = c(21)) + scale_size_manual(name = "Observation", + values = c(3)) + labs(x = var.name, y = "Probability density\n(total area=1)", + title = title) + theme(axis.text.x = element_blank(), panel.grid.minor.x = element_blank(), + legend.key.size = unit(0.3, "in"), panel.border = element_rect(fill = NA, + color = "black"), strip.background = element_rect(colour = "black", fill = "gray80"), + panel.spacing = unit(0.2, "in"), panel.grid.major.x = element_line(color = "grey93")) + + guides(fill = guide_legend(order = 1), color = guide_legend(order = 2, reverse = T), + shape = guide_legend(order = 3, label = F), size = guide_legend(order = 4, + label = F)) + #------------------------ + # Save to plotfile if needed, and return plot + #------------------------ + if (!is.null(plotfile)) { + ggsave(plotfile, plot) + } + return(plot) } - -.jitter.ensmemb <- function(x, thr=0.1) { - # Idea: start with first level. Loop all points, - # and if distance to last point in the level is more than a threshold, - # include the point to the level. - # Otherwise keep the point for another round. - # Do one round in each direction to avoid uggly patterns. - if (is.unsorted(x, na.rm = T)) { - stop("Provide a sorted array!") - } - - lev <- x * 0 - lev[is.na(lev)] <- -1 - level <- 1 - while (any(lev == 0)) { - last <- -1 / 0 - for (i in 1:length(x)) { - if (lev[i] != 0) { - next - } - if (x[i] - last > thr) { - lev[i] <- level - last <- x[i] - } +.jitter.ensmemb <- function(x, thr = 0.1) { + # Idea: start with first level. Loop all points, and if distance to last point in + # the level is more than a threshold, include the point to the level. Otherwise + # keep the point for another round. Do one round in each direction to avoid + # uggly patterns. + if (is.unsorted(x, na.rm = T)) { + stop("Provide a sorted array!") } - level <- level + 1 - last <- 1 / 0 - for (i in seq(length(x), 1, -1)) { - if (lev[i] != 0) { - next - } - if (last - x[i] > thr) { - lev[i] <- level - last <- x[i] - } + lev <- x * 0 + lev[is.na(lev)] <- -1 + level <- 1 + while (any(lev == 0)) { + last <- -1/0 + for (i in 1:length(x)) { + if (lev[i] != 0) { + next + } + if (x[i] - last > thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level + 1 + last <- 1/0 + for (i in seq(length(x), 1, -1)) { + if (lev[i] != 0) { + next + } + if (last - x[i] > thr) { + lev[i] <- level + last <- x[i] + } + } + level <- level + 1 } - level <- level + 1 - } - lev[lev == -1] <- NA - return(lev * thr * sqrt(3) / 2) + lev[lev == -1] <- NA + return(lev * thr * sqrt(3)/2) } - -.polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven=F) { - halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) - cross <- halfplane[-1L] - halfplane[-length(halfplane)] - does.cross <- cross != 0 - if (!any(does.cross)) { - return() - } - x1 <- x[-length(x)][does.cross] - y1 <- y[-length(y)][does.cross] - x2 <- x[-1L][does.cross] - y2 <- y[-1L][does.cross] - t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1)) / (xd * (y2 - y1) - yd * (x2 - x1))) - o <- order(t) - tsort <- t[o] - crossings <- cumsum(cross[does.cross][o]) - if (fillOddEven) { - crossings <- crossings %% 2 - } - drawline <- crossings != 0 - lx <- x0 + xd * tsort - ly <- y0 + yd * tsort - lx1 <- lx[-length(lx)][drawline] - ly1 <- ly[-length(ly)][drawline] - lx2 <- lx[-1L][drawline] - ly2 <- ly[-1L][drawline] - return(data.frame(x = lx1, y = ly1, xend = lx2, yend = ly2)) +.polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven = F) { + halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) + cross <- halfplane[-1L] - halfplane[-length(halfplane)] + does.cross <- cross != 0 + if (!any(does.cross)) { + return() + } + x1 <- x[-length(x)][does.cross] + y1 <- y[-length(y)][does.cross] + x2 <- x[-1L][does.cross] + y2 <- y[-1L][does.cross] + t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * + (x2 - x1))) + o <- order(t) + tsort <- t[o] + crossings <- cumsum(cross[does.cross][o]) + if (fillOddEven) { + crossings <- crossings%%2 + } + drawline <- crossings != 0 + lx <- x0 + xd * tsort + ly <- y0 + yd * tsort + lx1 <- lx[-length(lx)][drawline] + ly1 <- ly[-length(ly)][drawline] + lx2 <- lx[-1L][drawline] + ly2 <- ly[-1L][drawline] + return(data.frame(x = lx1, y = ly1, xend = lx2, yend = ly2)) } - -.polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches=c(5, 1)) { - x <- c(x, x[1L]) - y <- c(y, y[1L]) - angle <- angle %% 180 - upi <- c(width.units, height.units) / inches - if (upi[1L] < 0) { - angle <- 180 - angle - } - if (upi[2L] < 0) { - angle <- 180 - angle - } - upi <- abs(upi) - xd <- cos(angle / 180 * pi) * upi[1L] - yd <- sin(angle / 180 * pi) * upi[2L] - hatch.ls <- list() - i <- 1 - if (angle < 45 || angle > 135) { - if (angle < 45) { - first.x <- max(x) - last.x <- min(x) - } else { - first.x <- min(x) - last.x <- max(x) +.polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches = c(5, + 1)) { + x <- c(x, x[1L]) + y <- c(y, y[1L]) + angle <- angle%%180 + upi <- c(width.units, height.units)/inches + if (upi[1L] < 0) { + angle <- 180 - angle } - y.shift <- upi[2L] / density / abs(cos(angle / 180 * pi)) - x0 <- 0 - y0 <- floor((min(y) - first.x * yd / xd) / y.shift) * y.shift - y.end <- max(y) - last.x * yd / xd - while (y0 < y.end) { - hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) - i <- i + 1 - y0 <- y0 + y.shift + if (upi[2L] < 0) { + angle <- 180 - angle } - } else { - if (angle < 90) { - first.y <- max(y) - last.y <- min(y) + upi <- abs(upi) + xd <- cos(angle/180 * pi) * upi[1L] + yd <- sin(angle/180 * pi) * upi[2L] + hatch.ls <- list() + i <- 1 + if (angle < 45 || angle > 135) { + if (angle < 45) { + first.x <- max(x) + last.x <- min(x) + } else { + first.x <- min(x) + last.x <- max(x) + } + y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) + x0 <- 0 + y0 <- floor((min(y) - first.x * yd/xd)/y.shift) * y.shift + y.end <- max(y) - last.x * yd/xd + while (y0 < y.end) { + hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i + 1 + y0 <- y0 + y.shift + } } else { - first.y <- min(y) - last.y <- max(y) - } - x.shift <- upi[1L] / density / abs(sin(angle / 180 * pi)) - x0 <- floor((min(x) - first.y * xd / yd) / x.shift) * x.shift - y0 <- 0 - x.end <- max(x) - last.y * xd / yd - while (x0 < x.end) { - hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) - i <- i + 1 - x0 <- x0 + x.shift + if (angle < 90) { + first.y <- max(y) + last.y <- min(y) + } else { + first.y <- min(y) + last.y <- max(y) + } + x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) + x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift + y0 <- 0 + x.end <- max(x) - last.y * xd/yd + while (x0 < x.end) { + hatch.ls[[i]] <- .polygon.onehatch(x, y, x0, y0, xd, yd) + i <- i + 1 + x0 <- x0 + x.shift + } } - } - return(do.call("rbind", hatch.ls)) + return(do.call("rbind", hatch.ls)) } -- GitLab From 1fd06ebd51a3c69a26357a2538f1ffb6fbae450d Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 20 Mar 2019 11:32:13 +0100 Subject: [PATCH 21/21] Automatic documentation updated --- man/PlotForecastPDF.Rd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd index 89fecb49..bd273aaf 100644 --- a/man/PlotForecastPDF.Rd +++ b/man/PlotForecastPDF.Rd @@ -10,7 +10,7 @@ PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, color.set = c("ggplot", "s2s4e", "hydro")) } \arguments{ -\item{fcst}{a dataframe or array containing all the ensember members for each frecast. If \code{"fcst"} is an array, it should have two labelled dimensions, and one of them should be \code{"members"}. If \code{"fcsts"} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members.} +\item{fcst}{a dataframe or array containing all the ensember members for each frecast. If \code{'fcst'} is an array, it should have two labelled dimensions, and one of them should be \code{'members'}. If \code{'fcsts'} is a data.frame, each column shoul be a separate forecast, with the rows beeing the different ensemble members.} \item{tercile.limits}{an array with P33 and P66 values that define the tercile categories.} @@ -26,9 +26,9 @@ PlotForecastPDF(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, \item{fcst.names}{(optional) an array of strings with the titles of each individual forecast.} -\item{add.ensmemb}{either to add the ensemble members \code{"above"} (default) or \code{"below"} the pdf, or not (\code{"no"}).} +\item{add.ensmemb}{either to add the ensemble members \code{'above'} (default) or \code{'below'} the pdf, or not (\code{'no'}).} -\item{color.set}{a selection of predefined color sets: use \code{"ggplot"} (default) for blue/green/red, \code{"s2s4e"} for blue/grey/orange, or \code{"hydro"} for yellow/gray/blue (suitable for precipitation and inflows).} +\item{color.set}{a selection of predefined color sets: use \code{'ggplot'} (default) for blue/green/red, \code{'s2s4e'} for blue/grey/orange, or \code{'hydro'} for yellow/gray/blue (suitable for precipitation and inflows).} } \value{ a ggplot object containing the plot. @@ -42,7 +42,7 @@ fcsts <- data.frame(fcst1 = rnorm(50), fcst2 = rnorm(50, 0.5, 1.2), PlotForecastPDF(fcsts,c(-1,1)) fcsts2 <- array(rnorm(100), dim = c(members = 20, fcst = 5)) PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), - fcst.names = paste0("random fcst ", 1 : 5), obs = 0.7) + fcst.names = paste0('random fcst ', 1 : 5), obs = 0.7) } \author{ Llorenç Lledó \email{llledo@bsc.es} -- GitLab