diff --git a/R/VizForecastPDF.R b/R/VizForecastPDF.R index c477bd6bb56447586e562a0d56ca94963f9ce6ab..39bc2068de4e06fac07de451fb9e15b561777efe 100644 --- a/R/VizForecastPDF.R +++ b/R/VizForecastPDF.R @@ -54,10 +54,27 @@ #'@importFrom s2dv InsertDim #'@export VizForecastPDF <- 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"), + title = "Set a title", + var.name = "Varname (units)", + title.legend = 'Probability of terciles', + ensemble.legend = 'Ensemble members', + obs.title = 'Observations', + title.extremes = 'Probability of extremes', + strings.extremes = c('Below p10', 'Above p90'), + strings.legend = c('Below normal', 'Near normal', 'Above normal'), + strings.obs = NULL, + xlab.title = 'Probability density', + title.cex = 1, + labs.cex = 1, + fcst.names = NULL, + fcst.names.cex = 1, + obs.lines = TRUE, + obs.size = 3, + add.ensmemb = c("above", "below", "no"), color.set = c("ggplot", "s2s4e", "hydro", "vitigeoss"), - memb_dim = 'member') { + memb_dim = 'member', + height = 5, width = 6, res = 300, plotfile = NULL) { + value <- init <- extremes <- x <- ymin <- ymax <- tercile <- NULL y <- xend <- yend <- yjitter <- MLT <- lab.pos <- NULL ggColorHue <- function(n) { @@ -122,11 +139,25 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU # Check observations #------------------------ 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") + if (is.array(obs) & length(dim(obs)) == 1) { + obs <- as.vector(obs) } + if (is.vector(obs)) { + if (!length(obs) %in% c(1, npanels)) { + stop("The number of observations should equal one or the number of forecasts") + } + if (!is.null(strings.obs) & !length(obs) == length(strings.obs)) { + stop('The lenghts of string.obs and obs must coincide') + } + } else if (is.array(obs)) { + if (!length(obs[1, ]) %in% c(1, npanels)) { + stop("The number of observations (second dimension) should equal one or the number of forecasts") + } + if (!is.null(strings.obs) & !length(obs[, 1]) == length(strings.obs)) { + stop('The lenghts of string.obs and obs must coincide') + } + colorObs <- rep(c(colorObs, brewer.pal(length(obs[, 1]) - 1, "Accent")), as.numeric(dim(obs)[2])) + } else {stop('Observations should be a vector or an array')} } #------------------------ # Check tercile limits @@ -199,9 +230,9 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU #------------------------ 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))) + 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 @@ -216,35 +247,36 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU 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("Above normal", "Normal", "Below normal")) + strings.legend[1], + ifelse(tmp.df$x < tercile.limits[tmp.df$PANEL, 2], + strings.legend[2], strings.legend[3])), + levels = rev(strings.legend)) #------------------------ # Get the height and width of a panel #------------------------ pan.width <- diff(range(tmp.df$x)) pan.height <- max(tmp.df$ymax) - magic.ratio <- 9 * pan.height/pan.width + magic.ratio <- 9 * pan.height / pan.width #------------------------ # Compute hatch coordinates for extremes #------------------------ if (!is.null(extreme.limits)) { tmp.df$extremes <- factor(ifelse(tmp.df$x < extreme.limits[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")) + strings.extremes[1], ifelse(tmp.df$x < extreme.limits[tmp.df$PANEL, 2], + strings.legend[2], strings.extremes[2])), + levels = c(strings.extremes[2], strings.legend[2], strings.extremes[1])) 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)) + 0)) # compute the hatches for this polygon hatches <- .polygon.fullhatch(tmp.df2$x, tmp.df2$y, angle = 60, density = 10, - width.units = pan.width, height.units = pan.height) + width.units = pan.width, height.units = pan.height) # add bottom segment end1 <- data.frame(x = x$x[1], y = x$ymax[1], xend = x$x[1], yend = 0) # add top segment end2 <- data.frame(x = x$x[length(x$x)], y = x$ymax[length(x$x)], xend = x$x[length(x$x)], - yend = 0) + yend = 0) return(rbind(hatches, end1, end2)) }) attr <- attr(hatch.ls, "split_labels") @@ -279,8 +311,8 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU #------------------------ if (!is.null(obs)) { tmp.dt <- data.table(tmp.df, key = c("init", "x")) - obs.dt <- data.table(init = factor(colnames(fcst.df), levels = colnames(fcst.df)), - value = rep(obs, dim(fcst.df)[2])) + obs.dt <- data.table(init = factor(colnames(fcst.df), levels = colnames(fcst.df)), + value = rep(obs, dim(fcst.df)[2])) setkey(obs.dt, init, value) obs.xy <- tmp.dt[obs.dt, roll = "nearest"] } @@ -295,24 +327,24 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU # Add hatches for extremes #------------------------ if (!is.null(extreme.limits)) { - if (nrow(hatch.df[hatch.df$extremes != "Normal", ]) == 0) { + if (nrow(hatch.df[hatch.df$extremes != strings.legend[2], ]) == 0) { warning("The provided extreme categories are outside the plot bounds. The extremes will not be drawn.") extreme.limits <- NULL } else { plot <- plot + - geom_segment(data = hatch.df[hatch.df$extremes != "Normal", ], - aes(x = x, y = y, - xend = xend, yend = yend, color = extremes)) + geom_segment(data = hatch.df[hatch.df$extremes != strings.legend[2], ], + aes(x = x, y = y, + xend = xend, yend = yend, color = extremes)) } } #------------------------ # Add obs line #------------------------ - if (!is.null(obs)) { + if (!is.null(obs) & isTRUE(obs.lines)) { plot <- plot + - geom_vline(data = obs.dt, - aes(xintercept = value), - linetype = "dashed", color = colorObs) + geom_vline(data = unique(obs.dt), + aes(xintercept = value), + linetype = "dashed", color = colorObs) } #------------------------ # Add ensemble members @@ -321,46 +353,45 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU 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) + + 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) - + aes(x = x, + y = -pan.height / 10 - magic.ratio * yjitter, + shape = ensemble.legend), + color = "black", fill = colorMember, alpha = 1) + } else if (add.ensmemb == "above") { 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) - + aes(x = x, + y = 0.7 * magic.ratio * yjitter, + shape = ensemble.legend), + 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) - } - #------------------------ + # if (!is.null(obs)) { + # plot <- plot + + # # this adds the obs diamond + # geom_point(data = obs.xy, + # aes(x = x, y = ymax, size = obs.title), + # shape = 23, color = "black", fill = colorObs, show.legend = F) + # } + # #------------------------ # 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)] + by = .(init, tercile)] # 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("Above normal", "Normal", "Below normal"))), - by = c("init", "tercile"), all.y = T) - pct[is.na(pct),"pct"] <- 0 + tercile = factor(strings.legend, levels = rev(strings.legend))), + 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) @@ -372,13 +403,13 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU #------------------------ if (!is.null(extreme.limits)) { pct2 <- tmp.dt[, .(pct = integrate(approxfun(x, ymax), lower = min(x), upper = max(x))$value), - by = .(init, extremes)] + by = .(init, extremes)] # 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("Above P90", "Normal", "Below P10"))), - by = c("init", "extremes"), all.y=T) - pct2[is.na(pct),"pct"] <- 0 + extremes = factor(c(strings.extremes[1], strings.legend[2], strings.extremes[2]), + levels = c(strings.extremes[2], strings.legend[2], strings.extremes[1]))), + 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) @@ -398,58 +429,79 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU } 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) + + 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") + 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 + 3.5), - vjust = -0.5, angle = -90, size = 3.2, - color = rep(colorLab, dim(fcst.df)[2])) + geom_text(data = pct2[extremes != strings.legend[2], ], + 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", - 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)) + + scale_fill_manual(name = title.legend, + values = colorFill, drop = F) + + scale_color_manual(name = title.extremes, + values = colorHatch) + + scale_shape_manual(name = ensemble.legend, + values = c(21)) + + # scale_size_manual(name = obs.title, + # values = c(3)) + labs(x = var.name, - y = "Probability density\n(total area=1)", - title = title) + + y = xlab.title, + 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"), - panel.background = element_rect(fill = "white"), - plot.background = element_rect(fill = "white", color = NA)) + + plot.title = element_text(size = 10*title.cex), + axis.title.x = element_text(size = 8*labs.cex), + axis.title.y = element_text(size = 8*labs.cex), + legend.title = element_text(size = 8*labs.cex), + legend.text = element_text(size = 8*labs.cex), + panel.grid.minor.x = element_blank(), + legend.key.size = unit(0.2, "in"), + panel.border = element_rect(fill = NA, color = "black"), + strip.background = element_rect(colour = "black", fill = "gray80"), + strip.text = element_text(size = 8*fcst.names.cex), + panel.spacing = unit(0.2, "in"), + panel.grid.major.x = element_line(color = "grey93"), + panel.background = element_rect(fill = "white"), + plot.background = element_rect(fill = "white", color = NA)) + 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)) + color = guide_legend(order = 2), + shape = guide_legend(order = 3, label = F), + size = guide_legend(order = 4, label = F)) + #------------------------ + # Add obs diamond + #------------------------ + if (!is.null(obs)) { + obs.xy <- obs.xy[match(obs, obs.xy$x), ] + obs.xy$obs.factor <- factor(obs.xy$x, levels = obs) + obs.xy$obs.color <- colorObs + plot <- plot + new_scale(new_aes = "fill") + + geom_point(data = obs.xy, aes(x = x, y = ymax, fill = obs.color), + shape = 23, size = obs.size, color = "black") + + scale_fill_manual(name = obs.title, + values = unique(colorObs), + breaks = unique(colorObs), + labels = strings.obs) + } #------------------------ # Save to plotfile if needed, and return plot #------------------------ if (!is.null(plotfile)) { - ggsave(plotfile, plot) + ggsave(filename = plotfile, plot = plot, height = height, width = width, dpi = res) } return(plot) } @@ -465,7 +517,7 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU lev[is.na(lev)] <- -1 level <- 1 while (any(lev == 0)) { - last <- -1/0 + last <- -1 / 0 for (i in 1:length(x)) { if (lev[i] != 0) { next @@ -476,7 +528,7 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU } } level <- level + 1 - last <- 1/0 + last <- 1 / 0 for (i in seq(length(x), 1, -1)) { if (lev[i] != 0) { next @@ -489,7 +541,7 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU level <- level + 1 } lev[lev == -1] <- NA - return(lev * thr * sqrt(3)/2) + return(lev * thr * sqrt(3) / 2) } .polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven = F) { halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) @@ -502,13 +554,13 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU y1 <- y[-length(y)][does.cross] x2 <- x[-1L][does.cross] y2 <- y[-1L][does.cross] - t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * - (x2 - x1))) + t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1)) / (xd * (y2 - y1) - yd * + (x2 - x1))) o <- order(t) tsort <- t[o] crossings <- cumsum(cross[does.cross][o]) if (fillOddEven) { - crossings <- crossings%%2 + crossings <- crossings %% 2 } drawline <- crossings != 0 lx <- x0 + xd * tsort @@ -520,11 +572,11 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU return(data.frame(x = lx1, y = ly1, xend = lx2, yend = ly2)) } .polygon.fullhatch <- function(x, y, density, angle, width.units, height.units, inches = c(5, - 1)) { + 1)) { x <- c(x, x[1L]) y <- c(y, y[1L]) - angle <- angle%%180 - upi <- c(width.units, height.units)/inches + angle <- angle %% 180 + upi <- c(width.units, height.units) / inches if (upi[1L] < 0) { angle <- 180 - angle } @@ -532,8 +584,8 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU angle <- 180 - angle } upi <- abs(upi) - xd <- cos(angle/180 * pi) * upi[1L] - yd <- sin(angle/180 * pi) * upi[2L] + xd <- cos(angle / 180 * pi) * upi[1L] + yd <- sin(angle / 180 * pi) * upi[2L] hatch.ls <- list() i <- 1 if (angle < 45 || angle > 135) { @@ -544,7 +596,7 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU first.x <- min(x) last.x <- max(x) } - y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) + y.shift <- upi[2L] / density / abs(cos(angle/180 * pi)) x0 <- 0 y0 <- floor((min(y) - first.x * yd/xd)/y.shift) * y.shift y.end <- max(y) - last.x * yd/xd @@ -573,4 +625,3 @@ VizForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NU } return(do.call("rbind", hatch.ls)) } -