diff --git a/NEWS.md b/NEWS.md index f0886f8ef4f3e2e78c3f2d3142555d60ad669add..1ecc1bfe4619c480b346763a5328e0a73abc0d06 100644 --- a/NEWS.md +++ b/NEWS.md @@ -21,8 +21,8 @@ + CST_SplitDims returns ordered output following ascending order provided in indices when it is numeric + qmap library moved from Imports to Depends + CST_QuantileMapping correctly handles exp_cor - + Figures resize option from vignettes has been removed. - + + Figures resize option from vignettes has been removed + + Vignette PlotForecastPDF updated plots ### CSTools 3.1.0 **Submission date to CRAN: 02-07-2020** diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R index 4cf6de6c82c370de6a370326f78873b3416edf2c..6311628916972bdd150e2f3682520b1e0fe2433a 100644 --- a/R/PlotForecastPDF.R +++ b/R/PlotForecastPDF.R @@ -48,23 +48,23 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ color.set <- match.arg(color.set) if (color.set == "s2s4e") { - colorFill <- c("#FF764D", "#b5b5b5", "#33BFD1") - colorHatch <- c("deepskyblue3", "indianred3") + colorFill <- c("#FF764D", "#b5b5b5", "#33BFD1") # AN, N, BN fill colors + colorHatch <- c("indianred3", "deepskyblue3") # AP90, BP10 line colors colorMember <- c("#ffff7f") colorObs <- "purple" - colorLab <- c("blue", "red") + colorLab <- c("red", "blue") # AP90, BP10 text colors } else if (color.set == "hydro") { colorFill <- c("#41CBC9", "#b5b5b5", "#FFAB38") - colorHatch <- c("darkorange1", "deepskyblue3") + colorHatch <- c("deepskyblue3", "darkorange1") colorMember <- c("#ffff7f") colorObs <- "purple" - colorLab <- c("darkorange3", "blue") + colorLab <- c("blue", "darkorange3") } else if (color.set == "ggplot") { colorFill <- ggColorHue(3) - colorHatch <- c("deepskyblue3", "indianred1") + colorHatch <- c("indianred3", "deepskyblue3") colorMember <- c("#ffff7f") colorObs <- "purple" - colorLab <- c("blue", "red") + colorLab <- c("red", "blue") } else { stop("Parameter 'color.set' should be one of ggplot/s2s4e/hydro") } @@ -108,7 +108,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (length(tercile.limits) != 2) { stop("Provide two tercile limits") } - tercile.limits <- InsertDim(tercile.limits, 1, npanels) + tercile.limits <- InsertDim(tercile.limits, 1, npanels, name = "new") } else if (is.array(tercile.limits)) { if (length(dim(tercile.limits)) == 2) { if (dim(tercile.limits)[2] != 2) { @@ -135,7 +135,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N if (length(extreme.limits) != 2) { stop("Provide two extreme limits") } - extreme.limits <- InsertDim(extreme.limits, 1, npanels) + extreme.limits <- InsertDim(extreme.limits, 1, npanels, name = "new") } else if (is.array(extreme.limits)) { if (length(dim(extreme.limits)) == 2) { if (dim(extreme.limits)[2] != 2) { @@ -170,10 +170,11 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # 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 <- reshape2::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 @@ -187,9 +188,10 @@ 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("Above normal", "Normal", "Below normal")) #------------------------ # Get the height and width of a panel #------------------------ @@ -200,9 +202,10 @@ 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("Above P90", "Normal", "Below P10")) 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, @@ -236,10 +239,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # 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) + jitter.df <- reshape2::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) { + jitter.df$x <- reshape2::melt(data.frame(dlply(melt.df, .(init), function(x) { sort(x$value, na.last = T) })), value.name = "x", id.vars = NULL)$x } @@ -257,8 +260,10 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N #------------------------ # 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) + plot <- plot + + geom_ribbon(data = tmp.df, + aes(x = x, ymin = ymin, ymax = ymax, fill = tercile), + alpha = 0.7) #------------------------ # Add hatches for extremes #------------------------ @@ -267,37 +272,55 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N 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)) + plot <- plot + + geom_segment(data = hatch.df[hatch.df$extremes != "Normal", ], + aes(x = x, y = y, + xend = xend, yend = yend, color = extremes)) } } #------------------------ # Add obs line #------------------------ if (!is.null(obs)) { - plot <- plot + geom_vline(data = obs.dt, aes(xintercept = value), linetype = "dashed", - color = colorObs) + 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")) + 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, + aes(x = x, + y = -pan.height / 10 - magic.ratio * yjitter, + shape = "Ensemble members"), + color = "black", fill = colorMember, alpha = 1) + } 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")) + plot <- plot + + geom_point(data = jitter.df, + aes(x = x, + y = 0.7 * magic.ratio * yjitter, + shape = "Ensemble members"), + color = "black", fill = colorMember, alpha = 1) + } #------------------------ # 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) + 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 @@ -308,14 +331,15 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # include potentially missing groups pct <- merge(pct, CJ(init = factor(levels(pct$init), levels = levels(pct$init)), tercile = factor(c("Below normal", "Normal", "Above normal"), - levels = c("Below normal", "Normal", "Above normal"))), + levels = c("Above normal", "Normal", "Below normal"))), by = c("init", "tercile"), all.y = T) pct[is.na(pct),"pct"] <- 0 tot <- pct[, .(tot = sum(pct)), by = init] pct <- merge(pct, tot, by = "init") - pct$pct <- round(100 * pct$pct/pct$tot, 0) + pct$pct <- round(100 * pct$pct / pct$tot, 0) pct$MLT <- pct[, .(MLT = pct == max(pct)), by = init]$MLT - pct$lab.pos <- as.vector(apply(tercile.limits, 1, function(x) {c(min(x), mean(x), max(x))})) + pct <- pct[order(init, tercile)] + pct$lab.pos <- as.vector(apply(tercile.limits, 1, function(x) {c(max(x), mean(x), min(x))})) #------------------------ # Compute probability for extremes #------------------------ @@ -325,13 +349,14 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N # include potentially missing groups pct2 <- merge(pct2, CJ(init = factor(levels(pct2$init), levels = levels(pct2$init)), extremes = factor(c("Below P10", "Normal", "Above P90"), - levels = c("Below P10", "Normal", "Above P90"))), + levels = c("Above P90", "Normal", "Below P10"))), by = c("init", "extremes"), all.y=T) pct2[is.na(pct),"pct"] <- 0 tot2 <- pct2[, .(tot = sum(pct)), by = init] pct2 <- merge(pct2, tot2, by = "init") - pct2$pct <- round(100 * pct2$pct/pct2$tot, 0) - pct2$lab.pos <- as.vector(apply(extreme.limits, 1, function(x) {c(x[1], NA, x[2])})) + pct2$pct <- round(100 * pct2$pct / pct2$tot, 0) + pct2 <- pct2[order(init, extremes)] + pct2$lab.pos <- as.vector(apply(extreme.limits, 1, function(x) {c(x[2], NA, x[1])})) pct2 <- merge(pct2, max.df, by = c("init", "extremes"), all.x = T) } #------------------------ @@ -344,35 +369,53 @@ 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, 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, - y = labpos, label = "*", hjust = as.integer(tercile) * 3.5 - 5), vjust = 0.1, - angle = -90, size = 7, color = "black") + plot <- plot + + geom_text(data = pct, + aes(x = lab.pos, y = labpos, label = paste0(pct, "%"), + hjust = as.integer(tercile) * -1.5 + 3.5), + vjust = vjust, angle = -90, size = 3.2) + + geom_text(data = pct[MLT == T, ], + aes(x = lab.pos, y = labpos, label = "*", + hjust = as.integer(tercile) * -3.5 + 9), + 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.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])) + 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 + 3.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)) + plot <- plot + + theme_minimal() + + scale_fill_manual(name = "Probability of\nterciles", + 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), + shape = guide_legend(order = 3, label = F), + size = guide_legend(order = 4, label = F)) #------------------------ # Save to plotfile if needed, and return plot #------------------------ diff --git a/vignettes/Figures/PlotForecastPDF_ex2.png b/vignettes/Figures/PlotForecastPDF_ex2.png index 67b4af1e1b523c564047fb545baecc11a4dbac2c..e9623e5f5ae2131b88ffffa8712bfd456888cbaf 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex2.png and b/vignettes/Figures/PlotForecastPDF_ex2.png differ diff --git a/vignettes/Figures/PlotForecastPDF_ex3.png b/vignettes/Figures/PlotForecastPDF_ex3.png index 10b1e63b4766e8aff8d10419b849d8677f197c58..50b0f847ea03fc3962965ee6a73abcae76a3f1f6 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex3.png and b/vignettes/Figures/PlotForecastPDF_ex3.png differ diff --git a/vignettes/Figures/PlotForecastPDF_ex4.png b/vignettes/Figures/PlotForecastPDF_ex4.png index 34df30d73e73f1d1d8e5ff4c34e57c2d487f4c2f..f7723a8da66f2751dc6c2e53793f95824254546b 100644 Binary files a/vignettes/Figures/PlotForecastPDF_ex4.png and b/vignettes/Figures/PlotForecastPDF_ex4.png differ diff --git a/vignettes/PlotForecastPDF.Rmd b/vignettes/PlotForecastPDF.Rmd index 98b2ae144003e7d876689a827325cdf0eadf9fdc..757f4047440b7ad59f25aaaa310557767214f1fc 100644 --- a/vignettes/PlotForecastPDF.Rmd +++ b/vignettes/PlotForecastPDF.Rmd @@ -30,20 +30,28 @@ 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. +Input data can also be provided as an two-dimensional array, as far as one of the dimensions is named 'members': + +```{r,fig.show = 'hide',warning=F} +fcst <- array(rnorm(mean=25, sd=2, n=90), dim=c(members=30, 3)) +PlotForecastPDF(fcst, tercile.limits = c(23, 27)) +``` + +### 2.- Customizing the appearance of your plots +Some parameters allow to customize your plot by changing the title, the forecast labels, the variable name and units, or the colors. ```{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")) + title = "Forecasts valid for 2019-01-01 at Sunny Hills", + fcst.names = c("model a", "model b"), + color.set = "s2s4e") ``` ![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. +Optionally, we can include the probability of extreme values or the actually observed values. The tercile limits, extreme limits and observation values can be specified for each panel separately. ```{r,fig.show = 'hide',warning=F} fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30), @@ -51,30 +59,31 @@ fcst <- data.frame(fcst1 = rnorm(mean = 25, 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))) + extreme.limits = rbind(c(18, 28), c(20, 30), c(12, 24)), + color.set="s2s4e") ``` -The same example using a forecast in array format is provided. -```{r,fig.show = 'hide',warning=F} -fcst <- array(cbind(cbind(rnorm(mean = 25, sd = 3, n = 30), - rnorm(mean = 23, sd = 4.5, n = 30)), rnorm(mean = 17, sd = 3, n = 30)), - dim = c(members = 30, 3)) -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) -### 4.- Example using lonlat_data -An example using the lonlat_data from CSTools is provided. +### 4.- Saving your plot to a file +PlotForecastPDF uses ggplot2, so you can save the output of the function to a variable and operate with it as a ggplot2 object. For instance, you can save it to a file: + +``` +library(ggplot2) +fcst <- array(rnorm(mean=25, sd=2, n=90), dim=c(members=30, 3)) +plot <-PlotForecastPDF(fcst, tercile.limits = c(23, 27)) +ggsave("outfile.pdf", plot, width=7, height=5) +``` + +### 5.- A reproducible example using lonlat_data +This final example uses the sample lonlat data from CSTools. It is suitable for checking reproducibility of results. ```{r,fig.show = 'hide',warning=F} fcst <- data.frame(fcst1 = lonlat_data$exp$data[1,,1,1,1,1] - 273.15, fcst2 = lonlat_data$exp$data[1,,1,2,1,1] - 273.15) PlotForecastPDF(fcst, tercile.limits = c(5, 7), extreme.limits = c(4, 8), var.name = "Temperature (ºC)", - title = "Forecasts valid on 2000-11 at sample mediterranean region", + title = "Forecasts initialized on Nov 2000 at sample Mediterranean region", fcst.names = c("November", "December")) ``` ![Example 4](./Figures/PlotForecastPDF_ex4.png)