From 55851826e9f7618e229c6a77a91f092d17ebd67f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Mon, 28 Oct 2019 12:03:09 +0100 Subject: [PATCH 01/21] Typo in description. Corrected teerciles definition. --- R/PlotForecastPDF.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index cf531063..1bb74b63 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -3,8 +3,8 @@ #'@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. +#'@param fcst a dataframe or array containing all the ensember members for each forecast. 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 for each panel. Should have 2*npanels dimensions. #'@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). -- GitLab From 50132946d20d00e5fab3c29b6f055a722117e80c Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Wed, 27 Nov 2019 13:01:36 +0100 Subject: [PATCH 02/21] extreme.limits and tercile.limits for length(tercile.limits) > 1 --- R/PlotForecastPDF.R | 140 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 120 insertions(+), 20 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 1bb74b63..e2d801c5 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -70,21 +70,47 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # 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") + + #add.ensmemb <- match.arg(add.ensmemb) + if(!is.array(tercile.limits) & !is.array(extreme.limits)){ + stop("tercile.limits or extreme.limits are not arrays") } - if (tercile.limits[1] >= tercile.limits[2]) { - stop("The provided tercile limits are in the wrong order") + + for( j in 1:length(tercile.limits)){ + + if (length(tercile.limits[[j]]) != 2) { + stop("Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") } - 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") - } + if (tercile.limits[[j]][1] >= tercile.limits[[j]][2]) { + stop("The provided tercile limits are in the wrong order") } + if (!is.null(extreme.limits[[j]])) { + if (length(extreme.limits[[j]]) != 2) { + stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") + } + if (extreme.limits[[j]][1] >= tercile.limits[[j]][1] | extreme.limits[[j]][2] <= + tercile.limits[[j]][2]) { + stop("The provided extreme limits are not consistent with tercile limits") + } + } + + } + +# 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 (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 #------------------------ @@ -124,6 +150,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # 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)] @@ -132,9 +159,22 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } 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")) +tercile__limits <- abind(tercile.limits,along=0) #aqui . per _ + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile__limits[tmp.df$PANEL,1], + "Below normal", + ifelse(tmp.df$x < tercile__limits[tmp.df$PANEL,2],"Normal", "Above normal")), + levels = c("Below normal","Normal", "Above normal")) +# if (!is.null(ggp$layout$layout)) { +# tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] +# 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 #------------------------ @@ -144,10 +184,13 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # 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")) +extreme__limits <- abind(extreme.limits,along=0) #aqui . per _ + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,1], + "Below P10", + ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,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, @@ -177,6 +220,39 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } max.df <- do.call("rbind", max.ls) } +# 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 #------------------------ @@ -247,6 +323,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # 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)] @@ -254,22 +331,45 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N 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]][1], mean(tercile.limits[[1]]), tercile.limits[[1]][2]) + + # 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]) + lab.pos2 <- c(extreme.limits[[1]][1], NA, extreme.limits[[1]][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"))), ] } + +# 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 #------------------------ -- GitLab From c7a4f7d5b050426f8b876b8e294219217375d741 Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Tue, 3 Dec 2019 14:43:30 +0100 Subject: [PATCH 03/21] try different data type for terciles and extremes --- R/PlotForecastPDF.R | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index e2d801c5..f615a2b2 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -45,7 +45,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Define color sets #------------------------ - color.set <- match.arg(color.set) + # color.set <- match.arg(color.set) if (color.set == "s2s4e") { colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1")) colorHatch <- c("deepskyblue3", "indianred3") @@ -71,11 +71,21 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Check input arguments #------------------------ + + # tercile.limits <- as.array(InsertDim(tercile.limits,2,1)) + # extreme.limits <- as.array(InsertDim(extreme.limits,2,1)) + # tercile.limits <- as.data.table(tercile.limits) + # extreme.limits <- as.data.table(extreme.limits) + #add.ensmemb <- match.arg(add.ensmemb) - if(!is.array(tercile.limits) & !is.array(extreme.limits)){ - stop("tercile.limits or extreme.limits are not arrays") - } + # if(!is.array(tercile.limits) & !is.array(extreme.limits)){ + # stop("tercile.limits or extreme.limits are not arrays") + # } + # tercile.limits <- as.data.table(tercile_limits) + # extreme.limits <- as.data.table(extreme_limits) + + for( j in 1:length(tercile.limits)){ if (length(tercile.limits[[j]]) != 2) { -- GitLab From 3398bb075c01a39c16ce0aaadaa3b0e2ae30ffcd Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Wed, 4 Dec 2019 11:50:23 +0100 Subject: [PATCH 04/21] PlotForecastPDF now works in both different terciles and extremes and unique terciles and extremes --- R/PlotForecastPDF.R | 49 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 41 insertions(+), 8 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index f615a2b2..4bed74ae 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -42,6 +42,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N hues <- seq(15, 375, length = n + 1) hcl(h = hues, l = 65, c = 100)[1:n] } + #------------------------ # Define color sets #------------------------ @@ -84,8 +85,15 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # tercile.limits <- as.data.table(tercile_limits) # extreme.limits <- as.data.table(extreme_limits) +if(!is.list(tercile.limits)){ + tercile.limits <- list(tercile.limits) +} + +if(!is.list(extreme.limits)){ + extreme.limits <- list(extreme.limits) +} - + for( j in 1:length(tercile.limits)){ if (length(tercile.limits[[j]]) != 2) { @@ -169,11 +177,23 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } else { stop("Cannot find PANELS in ggp object") } -tercile__limits <- abind(tercile.limits,along=0) #aqui . per _ + +if(length(tercile__limits==1)){ + + 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")) + + +} else { + + tercile__limits <- abind(tercile.limits,along=0) #aqui . per _ tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile__limits[tmp.df$PANEL,1], "Below normal", ifelse(tmp.df$x < tercile__limits[tmp.df$PANEL,2],"Normal", "Above normal")), levels = c("Below normal","Normal", "Above normal")) +} + # if (!is.null(ggp$layout$layout)) { # tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] # tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] @@ -195,12 +215,25 @@ tercile__limits <- abind(tercile.limits,along=0) #aqui . per _ # Compute hatch coordinates for extremes #------------------------ - if (!is.null(extreme.limits)) { -extreme__limits <- abind(extreme.limits,along=0) #aqui . per _ - tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,1], - "Below P10", - ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,2],"Normal", "Above P90")), - levels = c("Below P10","Normal", "Above P90")) +if (!is.null(extreme.limits)) { + +if(length(extreme.limits)==1){ + + 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")) + +} else { + + extreme__limits <- abind(extreme.limits,along=0) #aqui . per _ + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,1], + "Below P10", + ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,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, -- GitLab From 7b04bc0ba84c99e25dc95f5039c929d9f94c4c8b Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Wed, 4 Dec 2019 12:14:39 +0100 Subject: [PATCH 05/21] bug fixed --- R/PlotForecastPDF.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 4bed74ae..8ba81fb9 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -178,7 +178,7 @@ if(!is.list(extreme.limits)){ stop("Cannot find PANELS in ggp object") } -if(length(tercile__limits==1)){ +if(length(tercile__limits)==1){ 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", @@ -580,3 +580,4 @@ if(length(extreme.limits)==1){ } return(do.call("rbind", hatch.ls)) } + -- GitLab From 1a062eda9097f7e295c8d7978b67353a5e9cb319 Mon Sep 17 00:00:00 2001 From: froura Date: Mon, 9 Dec 2019 16:23:10 +0100 Subject: [PATCH 06/21] Input Extreme and tercile limits as array. Check dimensions and consistency. --- R/PlotForecastPDF.R | 196 ++++++++++++++------------------------------ 1 file changed, 63 insertions(+), 133 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 8ba81fb9..22bfd367 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -71,64 +71,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check input arguments #------------------------ - - - # tercile.limits <- as.array(InsertDim(tercile.limits,2,1)) - # extreme.limits <- as.array(InsertDim(extreme.limits,2,1)) - # tercile.limits <- as.data.table(tercile.limits) - # extreme.limits <- as.data.table(extreme.limits) - - #add.ensmemb <- match.arg(add.ensmemb) - # if(!is.array(tercile.limits) & !is.array(extreme.limits)){ - # stop("tercile.limits or extreme.limits are not arrays") - # } - - # tercile.limits <- as.data.table(tercile_limits) - # extreme.limits <- as.data.table(extreme_limits) -if(!is.list(tercile.limits)){ - tercile.limits <- list(tercile.limits) -} - -if(!is.list(extreme.limits)){ - extreme.limits <- list(extreme.limits) -} - - - for( j in 1:length(tercile.limits)){ - - if (length(tercile.limits[[j]]) != 2) { - stop("Parameter 'tercile.limits' should be an array with two limits for delimiting tercile categories") - } - if (tercile.limits[[j]][1] >= tercile.limits[[j]][2]) { - stop("The provided tercile limits are in the wrong order") - } - if (!is.null(extreme.limits[[j]])) { - if (length(extreme.limits[[j]]) != 2) { - stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") - } - if (extreme.limits[[j]][1] >= tercile.limits[[j]][1] | extreme.limits[[j]][2] <= - tercile.limits[[j]][2]) { - stop("The provided extreme limits are not consistent with tercile limits") - } - } - - } - -# 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 (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 #------------------------ @@ -147,6 +89,62 @@ if(!is.list(extreme.limits)){ } else { stop("Parameter 'fcst' should be an array or a data.frame") } + npanels <- ncol(fcst.df) + #------------------------ + # Check tercile limits + #------------------------ + if (!is.array(tercile.limits)) { + stop("Tercile limits should be an array") + } + if (length(dim(tercile.limits)) == 1) { + if (dim(tercile.limits) != 2) { + stop("Provide 2 tercile limits") + } + tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) + } else if (length(dim(tercile.limits)) == 2) { + if (dim(tercile.limits)[2] != 2) { + stop("Provide two tercile limits for each panel") + } + if (dim(tercile.limits)[1] != npanels) { + stop("The number of tercile limits does not match the number of forecasts provided") + } + } else { + stop("tercile limits should have one or two dimensions") + } + # check consistency of tercile limits + if (any(tercile.limits[,1]>=tercile.limits[,2])) { + stop("Inconsistent tercile limits") + } + #------------------------ + # Check extreme limits + #------------------------ + if (!is.null(extreme.limits)) { + if (!is.array(extreme.limits)) { + stop("Extreme limits should be an array") + } + if (length(dim(extreme.limits)) == 1) { + if (dim(extreme.limits) != 2) { + stop("Provide 2 extreme limits") + } + extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) + } else if (length(dim(extreme.limits)) == 2) { + if (dim(extreme.limits)[2] != 2) { + stop("Provide two extreme limits for each panel") + } + if (dim(extreme.limits)[1] != npanels) { + stop("The number of extreme limits does not match the number of forecasts provided") + } + } else { + stop("extreme limits should have one or two dimensions") + } + # Check that extreme limits are consistent with tercile limits + if (any(extreme.limits[,1]>=tercile.limits[,1])) { + stop("Inconsistent lower extreme limits") + } + if (any(extreme.limits[,2]<=tercile.limits[,2])) { + stop("Inconsistent higher extreme limits") + } + } #------------------------ # Set proper fcst names #------------------------ @@ -178,33 +176,11 @@ if(!is.list(extreme.limits)){ stop("Cannot find PANELS in ggp object") } -if(length(tercile__limits)==1){ - - 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")) - - -} else { - - tercile__limits <- abind(tercile.limits,along=0) #aqui . per _ - tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile__limits[tmp.df$PANEL,1], + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL,1], "Below normal", - ifelse(tmp.df$x < tercile__limits[tmp.df$PANEL,2],"Normal", "Above normal")), + ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL,2],"Normal", "Above normal")), levels = c("Below normal","Normal", "Above normal")) -} -# if (!is.null(ggp$layout$layout)) { -# tmp.df <- ggp$data[[1]][, c("x", "ymin", "ymax", "PANEL")] -# 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 #------------------------ @@ -214,26 +190,12 @@ if(length(tercile__limits)==1){ #------------------------ # Compute hatch coordinates for extremes #------------------------ - -if (!is.null(extreme.limits)) { - -if(length(extreme.limits)==1){ - - 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")) - -} else { - - extreme__limits <- abind(extreme.limits,along=0) #aqui . per _ - tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,1], - "Below P10", - ifelse(tmp.df$x < extreme__limits[tmp.df$PANEL,2],"Normal", "Above P90")), + if (!is.null(extreme.limits)) { + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL,1], + "Below P10", + ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL,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, @@ -263,39 +225,7 @@ if(length(extreme.limits)==1){ } max.df <- do.call("rbind", max.ls) } -# 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 #------------------------ -- GitLab From 54e4964ffad769211a8e5ac9eb0ef9e96b25fd52 Mon Sep 17 00:00:00 2001 From: froura Date: Mon, 9 Dec 2019 16:30:07 +0100 Subject: [PATCH 07/21] Corrected bug in labels of tercile/extreme probabilities. Now uses first panel positions for all labels. --- R/PlotForecastPDF.R | 26 ++------------------------ 1 file changed, 2 insertions(+), 24 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 22bfd367..2e308748 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -296,7 +296,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # 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)] @@ -304,45 +303,24 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N 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]][1], mean(tercile.limits[[1]]), tercile.limits[[1]][2]) + lab.pos <- c(tercile.limits[1,1], mean(tercile.limits[1,]), tercile.limits[1,2]) - # 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]][1], NA, extreme.limits[[1]][2]) + lab.pos2 <- c(extreme.limits[1,1], NA, extreme.limits[1,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"))), ] } -# 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 #------------------------ -- GitLab From d2141ef511fce5cacd9085e58dcd610e638af942 Mon Sep 17 00:00:00 2001 From: froura Date: Mon, 9 Dec 2019 16:42:56 +0100 Subject: [PATCH 08/21] Input extreme.limits and tercile.limits as vector or array --- R/PlotForecastPDF.R | 52 +++++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 2e308748..5f5ba451 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -93,23 +93,24 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check tercile limits #------------------------ - if (!is.array(tercile.limits)) { - stop("Tercile limits should be an array") - } - if (length(dim(tercile.limits)) == 1) { - if (dim(tercile.limits) != 2) { + if (is.vector(tercile.limits)) { + if (length(tercile.limits) != 2) { stop("Provide 2 tercile limits") } tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) - } else if (length(dim(tercile.limits)) == 2) { - if (dim(tercile.limits)[2] != 2) { - stop("Provide two tercile limits for each panel") - } - if (dim(tercile.limits)[1] != npanels) { - stop("The number of tercile limits does not match the number of forecasts provided") - } + } else if (is.array(tercile.limits)) { + if (length(dim(tercile.limits)) == 2) { + if (dim(tercile.limits)[2] != 2) { + stop("Provide two tercile limits for each panel") + } + if (dim(tercile.limits)[1] != npanels) { + stop("The number of tercile limits does not match the number of forecasts provided") + } + } else { + stop("tercile limits should have two dimensions") + } } else { - stop("tercile limits should have one or two dimensions") + stop("Tercile limits should be a vector of length two or an array of dimension (nfcsts,2)") } # check consistency of tercile limits if (any(tercile.limits[,1]>=tercile.limits[,2])) { @@ -119,23 +120,24 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Check extreme limits #------------------------ if (!is.null(extreme.limits)) { - if (!is.array(extreme.limits)) { - stop("Extreme limits should be an array") - } - if (length(dim(extreme.limits)) == 1) { - if (dim(extreme.limits) != 2) { + if (is.vector(extreme.limits)) { + if (length(extreme.limits) != 2) { stop("Provide 2 extreme limits") } extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) - } else if (length(dim(extreme.limits)) == 2) { - if (dim(extreme.limits)[2] != 2) { - stop("Provide two extreme limits for each panel") - } - if (dim(extreme.limits)[1] != npanels) { - stop("The number of extreme limits does not match the number of forecasts provided") + } else if (is.array(extreme.limits)) { + if (length(dim(extreme.limits)) == 2) { + if (dim(extreme.limits)[2] != 2) { + stop("Provide two extreme limits for each panel") + } + if (dim(extreme.limits)[1] != npanels) { + stop("The number of extreme limits does not match the number of forecasts provided") + } + } else { + stop("extreme limits should have two dimensions") } } else { - stop("extreme limits should have one or two dimensions") + stop("Extreme limits should be a vector of length two or an array of dimensions (nfcsts,2)") } # Check that extreme limits are consistent with tercile limits if (any(extreme.limits[,1]>=tercile.limits[,1])) { -- GitLab From 4ce76375aff33d3406092b683073155f5e5aab8a Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Thu, 12 Dec 2019 13:10:17 +0100 Subject: [PATCH 09/21] Modification in obs parameter description. Obs format is checked --- R/PlotForecastPDF.R | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 5f5ba451..1f1d6f84 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -6,7 +6,7 @@ #'@param fcst a dataframe or array containing all the ensember members for each forecast. 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 for each panel. Should have 2*npanels dimensions. #'@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 obs (optional) A vector of length qual to the number of forecasts. A number with the observed value is accepted if observations are equal in each forecast. (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 var.name a string with the variable name and units. @@ -91,6 +91,14 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } npanels <- ncol(fcst.df) #------------------------ + # Check observations + #------------------------ + if (is.vector(obs) & length(obs) != 1 ) { + if (length(obs) != length(fcst)) { + stop("The number of observations differ from the number of forecasts.") + } + } + #------------------------ # Check tercile limits #------------------------ if (is.vector(tercile.limits)) { -- GitLab From a5d588e60a8447f1e9750d86c454e5dbcc78805b Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Thu, 12 Dec 2019 17:07:57 +0100 Subject: [PATCH 10/21] Modification in tercile.limits description. tercile.limitbeing 1-D array is now accepted --- R/PlotForecastPDF.R | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 1f1d6f84..9ad21ae4 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -4,7 +4,7 @@ #'@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 forecast. 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 for each panel. Should have 2*npanels dimensions. +#'@param tercile.limits an array with P33 and P66 values that define the tercile categories for each panel. Should have (nforecasts*2) dimensions. An unidimensional array with two elements or a 2-element vector are also accepted. #'@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 vector of length qual to the number of forecasts. A number with the observed value is accepted if observations are equal in each forecast. (Default: observation is not shown). #'@param plotfile (optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved). @@ -101,12 +101,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check tercile limits #------------------------ - if (is.vector(tercile.limits)) { - if (length(tercile.limits) != 2) { - stop("Provide 2 tercile limits") - } - tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) - } else if (is.array(tercile.limits)) { + if (is.array(tercile.limits)) { if (length(dim(tercile.limits)) == 2) { if (dim(tercile.limits)[2] != 2) { stop("Provide two tercile limits for each panel") @@ -114,11 +109,19 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (dim(tercile.limits)[1] != npanels) { stop("The number of tercile limits does not match the number of forecasts provided") } - } else { + } else if (length(dim(tercile.limits == 1))){ + tercile.limits <- as.vector(tercile.limits) + tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) + } else { stop("tercile limits should have two dimensions") + } + } else if (is.vector(tercile.limits)) { + if (length(tercile.limits) != 2) { + stop("Provide 2 tercile limits") } + tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) } else { - stop("Tercile limits should be a vector of length two or an array of dimension (nfcsts,2)") + stop("Tercile limits should be a vector of length two or an array of dimension (1,2) or (nfcsts,2)") } # check consistency of tercile limits if (any(tercile.limits[,1]>=tercile.limits[,2])) { -- GitLab From 17ad54046b2bbdbd0a10300e5bb532e9a01cabc3 Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Mon, 16 Dec 2019 11:51:48 +0100 Subject: [PATCH 11/21] Update of extreme.limits description. 1-dimensional array with 2 elements is now accepted --- R/PlotForecastPDF.R | 48 ++++++++++++++++++++++----------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 9ad21ae4..77ed1437 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -5,7 +5,7 @@ #' #'@param fcst a dataframe or array containing all the ensember members for each forecast. 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 for each panel. Should have (nforecasts*2) dimensions. An unidimensional array with two elements or a 2-element vector are also accepted. -#'@param extreme.limits (optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown). +#'@param extreme.limits (optional) an array with P10 and P90 values that define the extreme categories. Should have (nforecasts*2) dimensions. An unidimensional array with two elements or a 2-element vector are also accepted. (Default: extreme categories are not shown). #'@param obs (optional) A vector of length qual to the number of forecasts. A number with the observed value is accepted if observations are equal in each forecast. (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. @@ -131,29 +131,29 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Check extreme limits #------------------------ if (!is.null(extreme.limits)) { - if (is.vector(extreme.limits)) { - if (length(extreme.limits) != 2) { - stop("Provide 2 extreme limits") - } - extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) - } else if (is.array(extreme.limits)) { - if (length(dim(extreme.limits)) == 2) { - if (dim(extreme.limits)[2] != 2) { - stop("Provide two extreme limits for each panel") - } - if (dim(extreme.limits)[1] != npanels) { - stop("The number of extreme limits does not match the number of forecasts provided") - } - } else { - stop("extreme limits should have two dimensions") - } - } else { - stop("Extreme limits should be a vector of length two or an array of dimensions (nfcsts,2)") - } - # Check that extreme limits are consistent with tercile limits - if (any(extreme.limits[,1]>=tercile.limits[,1])) { - stop("Inconsistent lower extreme limits") - } + if (is.array(extreme.limits)) { + if (length(dim(extreme.limits)) == 2) { + if (dim(extreme.limits)[2] != 2) { + stop("Provide two extreme limits for each panel") + } + if (dim(extreme.limits)[1] != npanels) { + stop("The number of extreme limits does not match the number of forecasts provided") + } + } else if (length(dim(extreme.limits == 1))){ + extreme.limits <- as.vector(extreme.limits) + extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) + } else { + stop("extreme limits should have two dimensions") + } + } else if (is.vector(extreme.limits)) { + if (length(extreme.limits) != 2) { + stop("Provide 2 extreme limits") + } + extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) + } else { + stop("Extreme limits should be a vector of length two or an array of dimension (1,2) or (nfcsts,2)") + } +# check consistency of extreme limits if (any(extreme.limits[,2]<=tercile.limits[,2])) { stop("Inconsistent higher extreme limits") } -- GitLab From fc52bb69eab5fc53de3a7fae646959342222770d Mon Sep 17 00:00:00 2001 From: froura Date: Tue, 17 Dec 2019 10:58:28 +0100 Subject: [PATCH 12/21] Update PlotForecastPDF.R to enhance documentation --- R/PlotForecastPDF.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 77ed1437..b206b9ba 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -4,9 +4,9 @@ #'@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 forecast. 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 for each panel. Should have (nforecasts*2) dimensions. An unidimensional array with two elements or a 2-element vector are also accepted. -#'@param extreme.limits (optional) an array with P10 and P90 values that define the extreme categories. Should have (nforecasts*2) dimensions. An unidimensional array with two elements or a 2-element vector are also accepted. (Default: extreme categories are not shown). -#'@param obs (optional) A vector of length qual to the number of forecasts. A number with the observed value is accepted if observations are equal in each forecast. (Default: observation is not shown). +#'@param tercile.limits an array or vector with P33 and P66 values that define the tercile categories for each panel. Use an array of dimensions (nforecasts,2) to define different terciles for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. +#'@param extreme.limits (optional) an array or vector with P10 and P90 values that define the extreme categories for each panel. Use an array of (nforecasts,2) to define different extreme limits for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. (Default: extreme categories are not shown). +#'@param obs (optional) A vector providing the observed values for each forecast panel or a single value that will be reused for all forecast panels. (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 var.name a string with the variable name and units. -- GitLab From 501f79fdb58981cee407bc2b297514b287c0d77d Mon Sep 17 00:00:00 2001 From: froura Date: Tue, 17 Dec 2019 11:13:13 +0100 Subject: [PATCH 13/21] Update PlotForecastPDF.R comprovation that observations are a vector or a number. --- R/PlotForecastPDF.R | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index b206b9ba..7209be9a 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -93,11 +93,15 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check observations #------------------------ - if (is.vector(obs) & length(obs) != 1 ) { - if (length(obs) != length(fcst)) { - stop("The number of observations differ from the number of forecasts.") - } - } + if (!is.vector(obs) | !is.numeric(obs) ){ + stop("Observations is not a vector or a numeric value") + } else { + if (is.vector(obs) & length(obs) != 1 ) { + if (length(obs) != length(fcst)) { + stop("The number of observations differ from the number of forecasts.") + } + } + } #------------------------ # Check tercile limits #------------------------ -- GitLab From ad6e4134affd6031a6c95edd618a45fa0c77dc27 Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Tue, 17 Dec 2019 11:35:48 +0100 Subject: [PATCH 14/21] If there is a unique forecast, extreme limits and tercile limits must be vectors. Now 1-D array is cot accepted. Observations are proven to be a vector or a numeric value. --- R/PlotForecastPDF.R | 69 ++++++++++++++++++++++----------------------- 1 file changed, 33 insertions(+), 36 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 7209be9a..ac6fb366 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -94,7 +94,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Check observations #------------------------ if (!is.vector(obs) | !is.numeric(obs) ){ - stop("Observations is not a vector or a numeric value") + stop("Observations is not a vector nor a numeric value") } else { if (is.vector(obs) & length(obs) != 1 ) { if (length(obs) != length(fcst)) { @@ -105,7 +105,12 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check tercile limits #------------------------ - if (is.array(tercile.limits)) { + if (is.vector(tercile.limits)) { + if (length(tercile.limits) != 2) { + stop("Provide 2 tercile limits") + } + tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) + } else if (is.array(tercile.limits)) { if (length(dim(tercile.limits)) == 2) { if (dim(tercile.limits)[2] != 2) { stop("Provide two tercile limits for each panel") @@ -113,19 +118,11 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (dim(tercile.limits)[1] != npanels) { stop("The number of tercile limits does not match the number of forecasts provided") } - } else if (length(dim(tercile.limits == 1))){ - tercile.limits <- as.vector(tercile.limits) - tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) - } else { - stop("tercile limits should have two dimensions") - } - } else if (is.vector(tercile.limits)) { - if (length(tercile.limits) != 2) { - stop("Provide 2 tercile limits") + } else { + stop("Tercile limits should have two dimensions") } - tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) } else { - stop("Tercile limits should be a vector of length two or an array of dimension (1,2) or (nfcsts,2)") + stop("Tercile limits should be a vector of length two or an array of dimension (nfcsts,2)") } # check consistency of tercile limits if (any(tercile.limits[,1]>=tercile.limits[,2])) { @@ -135,29 +132,29 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Check extreme limits #------------------------ if (!is.null(extreme.limits)) { - if (is.array(extreme.limits)) { - if (length(dim(extreme.limits)) == 2) { - if (dim(extreme.limits)[2] != 2) { - stop("Provide two extreme limits for each panel") - } - if (dim(extreme.limits)[1] != npanels) { - stop("The number of extreme limits does not match the number of forecasts provided") - } - } else if (length(dim(extreme.limits == 1))){ - extreme.limits <- as.vector(extreme.limits) - extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) - } else { - stop("extreme limits should have two dimensions") - } - } else if (is.vector(extreme.limits)) { - if (length(extreme.limits) != 2) { - stop("Provide 2 extreme limits") - } - extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) - } else { - stop("Extreme limits should be a vector of length two or an array of dimension (1,2) or (nfcsts,2)") - } -# check consistency of extreme limits + if (is.vector(extreme.limits)) { + if (length(extreme.limits) != 2) { + stop("Provide 2 extreme limits") + } + extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) + } else if (is.array(extreme.limits)) { + if (length(dim(extreme.limits)) == 2) { + if (dim(extreme.limits)[2] != 2) { + stop("Provide two extreme limits for each panel") + } + if (dim(extreme.limits)[1] != npanels) { + stop("The number of extreme limits does not match the number of forecasts provided") + } + } else { + stop("extreme limits should have two dimensions") + } + } else { + stop("Extreme limits should be a vector of length two or an array of dimensions (nfcsts,2)") + } + # Check that extreme limits are consistent with tercile limits + if (any(extreme.limits[,1]>=tercile.limits[,1])) { + stop("Inconsistent lower extreme limits") + } if (any(extreme.limits[,2]<=tercile.limits[,2])) { stop("Inconsistent higher extreme limits") } -- GitLab From 4bd71da50b3968335c807c7070c21c597d55d51d Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Mon, 20 Jan 2020 12:13:30 +0100 Subject: [PATCH 15/21] error message has been changed --- R/PlotForecastPDF.R | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index ac6fb366..5514dce3 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -93,8 +93,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check observations #------------------------ - if (!is.vector(obs) | !is.numeric(obs) ){ - stop("Observations is not a vector nor a numeric value") + if (is.array(obs)){ + stop("Observations can not be an array.") } else { if (is.vector(obs) & length(obs) != 1 ) { if (length(obs) != length(fcst)) { @@ -107,7 +107,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ if (is.vector(tercile.limits)) { if (length(tercile.limits) != 2) { - stop("Provide 2 tercile limits") + stop("Parameter 'tercile.limits' should be an array with two limits for delimiting extreme categories") + # stop("Provide 2 tercile limits") } tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) } else if (is.array(tercile.limits)) { @@ -134,7 +135,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (!is.null(extreme.limits)) { if (is.vector(extreme.limits)) { if (length(extreme.limits) != 2) { - stop("Provide 2 extreme limits") + stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") + #stop("Provide 2 extreme limits") } extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) } else if (is.array(extreme.limits)) { -- GitLab From 117eaa0f062e97ad5c92dd89bde7fb7da02b62da Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Tue, 28 Jan 2020 16:15:19 +0100 Subject: [PATCH 16/21] changed testthat error messages for PlotForecastPDF --- R/PlotForecastPDF.R | 8 +++----- tests/testthat/test-PlotForecastPDF.R | 6 +++--- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 5514dce3..07781b2f 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -107,8 +107,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ if (is.vector(tercile.limits)) { if (length(tercile.limits) != 2) { - stop("Parameter 'tercile.limits' should be an array with two limits for delimiting extreme categories") - # stop("Provide 2 tercile limits") + stop("Provide two tercile limits") } tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) } else if (is.array(tercile.limits)) { @@ -134,9 +133,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ if (!is.null(extreme.limits)) { if (is.vector(extreme.limits)) { - if (length(extreme.limits) != 2) { - stop("Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") - #stop("Provide 2 extreme limits") + if (length(extreme.limits) != 2) { + stop("Provide two extreme limits") } extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) } else if (is.array(extreme.limits)) { diff --git a/tests/testthat/test-PlotForecastPDF.R b/tests/testthat/test-PlotForecastPDF.R index dfc35c4d..f9926324 100644 --- a/tests/testthat/test-PlotForecastPDF.R +++ b/tests/testthat/test-PlotForecastPDF.R @@ -2,13 +2,13 @@ context("Generic tests") test_that("Sanity checks", { expect_error( PlotForecastPDF(fcst, tercile.limits), - "object 'tercile.limits' not found") + "object 'fcst' 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") + "object 'fcst' not found") expect_error( PlotForecastPDF(fcst, tercile.limits = c(10, 20)), "object 'fcst' not found") @@ -18,5 +18,5 @@ test_that("Sanity checks", { "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") + "Provide two extreme limits") }) -- GitLab From 57d0a7a9105bf05dd98389d11d4fcdf1f4f94659 Mon Sep 17 00:00:00 2001 From: froura Date: Tue, 28 Jan 2020 16:22:33 +0100 Subject: [PATCH 17/21] Prova RMD vs MD --- vignettes/{MultivarRMSE_vignette.Rmd => MultivarRMSE_vignette.md} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename vignettes/{MultivarRMSE_vignette.Rmd => MultivarRMSE_vignette.md} (100%) diff --git a/vignettes/MultivarRMSE_vignette.Rmd b/vignettes/MultivarRMSE_vignette.md similarity index 100% rename from vignettes/MultivarRMSE_vignette.Rmd rename to vignettes/MultivarRMSE_vignette.md -- GitLab From 512c198d89134eaa726c381fd50812086ea4d1cb Mon Sep 17 00:00:00 2001 From: froura Date: Tue, 28 Jan 2020 16:24:10 +0100 Subject: [PATCH 18/21] Desfer prova chungui --- vignettes/{MultivarRMSE_vignette.md => MultivarRMSE_vignette.Rmd} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename vignettes/{MultivarRMSE_vignette.md => MultivarRMSE_vignette.Rmd} (100%) diff --git a/vignettes/MultivarRMSE_vignette.md b/vignettes/MultivarRMSE_vignette.Rmd similarity index 100% rename from vignettes/MultivarRMSE_vignette.md rename to vignettes/MultivarRMSE_vignette.Rmd -- GitLab From 9ec5d7ceadac41fa414298afd2087f5555f9c4c7 Mon Sep 17 00:00:00 2001 From: FRANCESC ROURA ADSERIAS Date: Tue, 28 Jan 2020 17:32:01 +0100 Subject: [PATCH 19/21] Added a vignette with three examples --- vignettes/PlotForecastPDF.Rmd | 47 +++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 vignettes/PlotForecastPDF.Rmd diff --git a/vignettes/PlotForecastPDF.Rmd b/vignettes/PlotForecastPDF.Rmd new file mode 100644 index 00000000..c6d89618 --- /dev/null +++ b/vignettes/PlotForecastPDF.Rmd @@ -0,0 +1,47 @@ +--- +author: "Francesc Roura and Llorenç Lledó" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Plot Forecast PDFs} + %\usepackage[utf8]{inputenc} +--- + +Plot Forecast PDFs (Probability Distibution Functions) +------------------------------------------ + +Library *CSTools*, should be installed from CRAN and loaded: + + +```{r,warning=FALSE,message=FALSE,error=FALSE} +library(CSTools) +``` + +### 1.- A simple example + +The first step is to put your forecasts in an appropriate format. For this vignette we generate some random values from two normal distributions. The PlotForecastPDF by default will plot the ensemble members, the estimated density distributions and the tercile probabilities. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=23,sd=4.5,n=30)) +PlotForecastPDF(fcst,tercile.limits=c(20,26)) +``` + +![Example 1](./Figures/PlotForecastPDF_ex1.png) + +### 2.- Some useful parameters +Changing the title, the forecast labels or the units will be needed in most cases. + +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=23,sd=4.5,n=30)) +PlotForecastPDF(fcst,tercile.limits=c(20,26),var.name="Temperature (ºC)",title="Forecasts valid on 2019-01-01 at Sunny Hills",fcst.names = c("model a","model b")) +``` +![Example 2](./Figures/PlotForecastPDF_ex2.png) + +### 3.- Adding extremes and observed values +We can add the probability of extreme values and the observed values. The tercile and extreme limits can be specified for each panel separately, as well as the observed values. +```{r,fig.show = 'hide',warning=F} +fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=28,sd=4.5,n=30),fcst3=rnorm(mean=17,sd=3,n=30)) +PlotForecastPDF(fcst,tercile.limits=rbind(c(20,26),c(22,28),c(15,22)),var.name="Temperature (ºC)",title="Forecasts at Sunny Hills",fcst.names = c("January","February","March"),obs=c(21,24,17),extreme.limits = rbind(c(18,28),c(20,30),c(12,24))) +``` +![Example 3](./Figures/PlotForecastPDF_ex3.png) -- GitLab From cf29303d046399a760e16fe7b104041f52050180 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lloren=C3=A7=20Lled=C3=B3?= Date: Fri, 31 Jan 2020 15:34:56 +0100 Subject: [PATCH 20/21] Fromatting and cleanining issues --- R/PlotForecastPDF.R | 69 +++++++++++++++++++-------------------------- 1 file changed, 29 insertions(+), 40 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 07781b2f..78ea53c2 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -1,10 +1,10 @@ #'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. 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. +#'@description This function plots the probability distribution function of several ensemble forecasts. Separate panels are used to plot forecasts valid or initialized at different times or by different models or even at different locations. 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 forecast. 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 or vector with P33 and P66 values that define the tercile categories for each panel. Use an array of dimensions (nforecasts,2) to define different terciles for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. +#'@param tercile.limits an array or vector with P33 and P66 values that define the tercile categories for each panel. Use an array of dimensions (nforecasts,2) to define different terciles for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. #'@param extreme.limits (optional) an array or vector with P10 and P90 values that define the extreme categories for each panel. Use an array of (nforecasts,2) to define different extreme limits for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. (Default: extreme categories are not shown). #'@param obs (optional) A vector providing the observed values for each forecast panel or a single value that will be reused for all forecast panels. (Default: observation is not shown). #'@param plotfile (optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved). @@ -23,6 +23,7 @@ #'@importFrom reshape2 melt #'@importFrom plyr . #'@importFrom plyr dlply +#'@importFrom s2dverification InsertDim #' #'@examples #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), @@ -42,11 +43,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N hues <- seq(15, 375, length = n + 1) hcl(h = hues, l = 65, c = 100)[1:n] } - #------------------------ # Define color sets #------------------------ - # color.set <- match.arg(color.set) + color.set <- match.arg(color.set) if (color.set == "s2s4e") { colorFill <- rev(c("#FF764D", "#b5b5b5", "#33BFD1")) colorHatch <- c("deepskyblue3", "indianred3") @@ -71,6 +71,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check input arguments #------------------------ + add.ensmemb <- match.arg(add.ensmemb) #------------------------ # Check fcst type and convert to data.frame if needed #------------------------ @@ -93,13 +94,11 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # Check observations #------------------------ - if (is.array(obs)){ - stop("Observations can not be an array.") - } else { - if (is.vector(obs) & length(obs) != 1 ) { - if (length(obs) != length(fcst)) { - stop("The number of observations differ from the number of forecasts.") - } + if (!is.null(obs)) { + if (!is.vector(obs)){ + stop("Observations should be a vector") + } else if (!length(obs) %in% c(1, npanels)) { + stop("The number of observations should equal one or the number of forecasts") } } #------------------------ @@ -107,9 +106,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ if (is.vector(tercile.limits)) { if (length(tercile.limits) != 2) { - stop("Provide two tercile limits") + stop("Provide two tercile limits") } - tercile.limits <- s2dverification::InsertDim(tercile.limits,1,npanels) + tercile.limits <- InsertDim(tercile.limits, 1, npanels) } else if (is.array(tercile.limits)) { if (length(dim(tercile.limits)) == 2) { if (dim(tercile.limits)[2] != 2) { @@ -125,7 +124,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N stop("Tercile limits should be a vector of length two or an array of dimension (nfcsts,2)") } # check consistency of tercile limits - if (any(tercile.limits[,1]>=tercile.limits[,2])) { + if (any(tercile.limits[, 1] >= tercile.limits[, 2])) { stop("Inconsistent tercile limits") } #------------------------ @@ -133,10 +132,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ if (!is.null(extreme.limits)) { if (is.vector(extreme.limits)) { - if (length(extreme.limits) != 2) { - stop("Provide two extreme limits") + if (length(extreme.limits) != 2) { + stop("Provide two extreme limits") } - extreme.limits <- s2dverification::InsertDim(extreme.limits,1,npanels) + extreme.limits <- InsertDim(extreme.limits, 1, npanels) } else if (is.array(extreme.limits)) { if (length(dim(extreme.limits)) == 2) { if (dim(extreme.limits)[2] != 2) { @@ -152,10 +151,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N stop("Extreme limits should be a vector of length two or an array of dimensions (nfcsts,2)") } # Check that extreme limits are consistent with tercile limits - if (any(extreme.limits[,1]>=tercile.limits[,1])) { + if (any(extreme.limits[, 1] >= tercile.limits[, 1])) { stop("Inconsistent lower extreme limits") } - if (any(extreme.limits[,2]<=tercile.limits[,2])) { + if (any(extreme.limits[, 2] <= tercile.limits[, 2])) { stop("Inconsistent higher extreme limits") } } @@ -180,7 +179,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # 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)] @@ -189,12 +187,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } else { stop("Cannot find PANELS in ggp object") } - - tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL,1], - "Below normal", - ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL,2],"Normal", "Above normal")), - levels = c("Below normal","Normal", "Above normal")) - + tmp.df$tercile <- factor(ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 1], "Below normal", + ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 2], "Normal", "Above normal")), levels = c("Below normal", + "Normal", "Above normal")) #------------------------ # Get the height and width of a panel #------------------------ @@ -205,11 +200,9 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # Compute hatch coordinates for extremes #------------------------ if (!is.null(extreme.limits)) { - tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL,1], - "Below P10", - ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL,2],"Normal", "Above P90")), - levels = c("Below P10","Normal", "Above P90")) - + tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 1], "Below P10", + ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 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, @@ -239,7 +232,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } max.df <- do.call("rbind", max.ls) } - #------------------------ # Compute jitter space for ensemble members #------------------------ @@ -317,8 +309,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N 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,1], mean(tercile.limits[1,]), tercile.limits[1,2]) - + pct$lab.pos <- as.vector(apply(tercile.limits, 1, function(x) {c(min(x), mean(x), max(x))})) #------------------------ # Compute probability for extremes #------------------------ @@ -328,13 +319,12 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N 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,1], NA, extreme.limits[1,2]) + pct2$lab.pos <- as.vector(apply(extreme.limits, 1, function(x) {c(x[1], NA, x[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 #------------------------ @@ -345,16 +335,16 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N labpos <- 0 vjust <- -0.5 } - plot <- plot + geom_text(data = pct, aes(x = lab.pos[as.integer(tercile)], y = labpos, + plot <- plot + geom_text(data = pct, aes(x = lab.pos, 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)], + angle = -90, size = 3.2) + geom_text(data = pct[MLT == T, ], aes(x = lab.pos, 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)], + plot <- plot + geom_text(data = pct2[extremes != "Normal", ], aes(x = lab.pos, 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])) @@ -502,4 +492,3 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N } return(do.call("rbind", hatch.ls)) } - -- GitLab From 24a31d2d35497989c41e49dd9482ac7b5cf530ce Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 7 Feb 2020 18:05:21 +0100 Subject: [PATCH 21/21] defining global variables andfixing dependency declaration with s2dverification --- R/PlotForecastPDF.R | 5 +++-- man/PlotForecastPDF.Rd | 10 +++++----- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 78ea53c2..c40cc32c 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -23,7 +23,7 @@ #'@importFrom reshape2 melt #'@importFrom plyr . #'@importFrom plyr dlply -#'@importFrom s2dverification InsertDim +#'@import s2dverification #' #'@examples #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), @@ -38,7 +38,8 @@ 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 + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- NULL + y <- xend <- yend <- yjitter <- MLT <- lab.pos <- NULL ggColorHue <- function(n) { hues <- seq(15, 375, length = n + 1) hcl(h = hues, l = 65, c = 100)[1:n] diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd index 3d1d65c3..d7b95b08 100644 --- a/man/PlotForecastPDF.Rd +++ b/man/PlotForecastPDF.Rd @@ -10,13 +10,13 @@ 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 forecast. 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{tercile.limits}{an array or vector with P33 and P66 values that define the tercile categories for each panel. Use an array of dimensions (nforecasts,2) to define different terciles for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels.} -\item{extreme.limits}{(optional) an array with P10 and P90 values that define the extreme categories. (Default: extreme categories are not shown).} +\item{extreme.limits}{(optional) an array or vector with P10 and P90 values that define the extreme categories for each panel. Use an array of (nforecasts,2) to define different extreme limits for each forecast panel, or a vector with two elements to reuse the same tercile limits for all forecast panels. (Default: extreme categories are not shown).} -\item{obs}{(optional) a number with the observed value. (Default: observation is not shown).} +\item{obs}{(optional) A vector providing the observed values for each forecast panel or a single value that will be reused for all forecast panels. (Default: observation is not shown).} \item{plotfile}{(optional) a filename (pdf, png...) where the plot will be saved. (Default: the plot is not saved).} @@ -34,7 +34,7 @@ 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. 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. +This function plots the probability distribution function of several ensemble forecasts. Separate panels are used to plot forecasts valid or initialized at different times or by different models or even at different locations. 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(10), fcst2 = rnorm(10, 0.5, 1.2), -- GitLab