######################################################### # Barcelona Supercomputing Center # # llledo@bsc.es - September 2018 # # Functions to plot several probabilistic forecasts # ######################################################### library(data.table) library(ggplot2) library(reshape2) library(plyr) #======================== # Function to plot several forecasts for an event, # either initialized at different moments or by different models. # Probabilities for extreme categories can be added as hatches. # Ensemble members can be added as jittered points. # The observed value can be added as a diamond over the pdf. #======================== plot_forecasts <- function(fcst_df,tercile_limits,extreme_limits=NULL,obs=NULL,plotfile=NULL,title="Set a title",varname="Varname (units)",fcst_names=NULL,add_ensmemb=c("above","below","no"),colors=c("ggplot","dst","hydro")) { #------------------------ # Color definitions #------------------------ colors <- match.arg(colors) if(colors=="dst") { colorFill <- rev(c("#FF764D","#b5b5b5","#33BFD1")) colorHatch <- c("deepskyblue3","indianred3") colorMember <- c("#ffff7f") colorObs <- "purple" colorLab <- c("blue","red") } else if(colors=="hydro") { colorFill <- rev(c("#41CBC9","#b5b5b5","#FFAB38")) colorHatch <- c("darkorange1","deepskyblue3") colorMember <- c("#ffff7f") colorObs <- "purple" colorLab <- c("darkorange3","blue") } else if(colors=="ggplot") { gg_color_hue <- function(n) { hues = seq(15, 375, length = n + 1) hcl(h = hues, l = 65, c = 100)[1:n] } colorFill <- rev(gg_color_hue(3)) colorHatch <- c("deepskyblue3","indianred1") colorMember <- c("#ffff7f") colorObs <- "purple" colorLab <- c("blue","red") } else { stop("Unknown color set") } #------------------------ # Check args #------------------------ add_ensmemb <- match.arg(add_ensmemb) if(length(tercile_limits)!=2) stop("Provide two limits for delimiting tercile categories") if(tercile_limits[1]>=tercile_limits[2]) stop("The provided tercile limits are in the wrong order") if(!is.null(extreme_limits)) { if(length(extreme_limits)!=2) stop("Provide two limits for delimiting extreme categories") if(extreme_limits[1]>=tercile_limits[1] | extreme_limits[2]<=tercile_limits[2]) stop("The provided extreme limits are not consistent with tercile limits") } #------------------------ # Set proper fcst names #------------------------ if(!is.null(fcst_names)) { colnames(fcst_df) <- factor(fcst_names,levels=fcst_names) } #------------------------ # Produce a first plot with the pdf for each init in a panel #------------------------ melt_df <- melt(fcst_df,variable.name="init") plot <- ggplot(melt_df,aes(x=value)) + geom_density(alpha=1,na.rm=T) + coord_flip() + facet_wrap(~init,strip.position="top",nrow=1) + xlim(range(c(obs,density(melt_df$value,na.rm=T)$x))) ggp <- ggplot_build(plot) #------------------------ # Gather the coordinates of the plots # together with init and corresponding terciles #------------------------ tmp.df <- ggp$data[[1]][,c("x","ymin","ymax","PANEL")] tmp.df$init <- ggp$layout$layout$init[as.numeric(tmp.df$PANEL)] tmp.df$tercile <- factor(ifelse(tmp.df$xthr) { lev[i] <- level last <- x[i] } } level <- level+1 last <- 1/0 for(i in seq(length(x),1,-1)) { if (lev[i]!=0) { next } if (last-x[i]>thr) { lev[i] <- level last <- x[i] } } level <- level+1 } return(lev*thr*sqrt(3)/2) } #================== # Hatch one polygon # Based on base polygon function #================== polygon.onehatch <- function(x, y, x0, y0, xd, yd, fillOddEven=F) { halfplane <- as.integer(xd * (y - y0) - yd * (x - x0) <= 0) cross <- halfplane[-1L] - halfplane[-length(halfplane)] does.cross <- cross != 0 if (!any(does.cross)) return() x1 <- x[-length(x)][does.cross] y1 <- y[-length(y)][does.cross] x2 <- x[-1L][does.cross] y2 <- y[-1L][does.cross] t <- (((x1 - x0) * (y2 - y1) - (y1 - y0) * (x2 - x1))/(xd * (y2 - y1) - yd * (x2 - x1))) o <- order(t) tsort <- t[o] crossings <- cumsum(cross[does.cross][o]) if (fillOddEven) crossings <- crossings%%2 drawline <- crossings != 0 lx <- x0 + xd * tsort ly <- y0 + yd * tsort lx1 <- lx[-length(lx)][drawline] ly1 <- ly[-length(ly)][drawline] lx2 <- lx[-1L][drawline] ly2 <- ly[-1L][drawline] return(data.frame(x=lx1,y=ly1,xend=lx2,yend=ly2)) } #================== # Hatch one polygon # Based on base polygon function #================== polygon.fullhatch <- function(x, y, density, angle, width_units, height_units, inches=c(5,1)) { x <- c(x, x[1L]) y <- c(y, y[1L]) angle <- angle%%180 #usr <- par("usr") #inches <- par("pin") upi <- c(width_units, height_units)/inches if (upi[1L] < 0) angle <- 180 - angle if (upi[2L] < 0) angle <- 180 - angle upi <- abs(upi) xd <- cos(angle/180 * pi) * upi[1L] yd <- sin(angle/180 * pi) * upi[2L] hatch_ls <- list() i <- 1 if (angle < 45 || angle > 135) { if (angle < 45) { first.x <- max(x) last.x <- min(x) } else { first.x <- min(x) last.x <- max(x) } y.shift <- upi[2L]/density/abs(cos(angle/180 * pi)) x0 <- 0 y0 <- floor((min(y) - first.x * yd/xd)/y.shift) * y.shift y.end <- max(y) - last.x * yd/xd while (y0 < y.end) { hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) i <- i+1 y0 <- y0 + y.shift } } else { if (angle < 90) { first.y <- max(y) last.y <- min(y) } else { first.y <- min(y) last.y <- max(y) } x.shift <- upi[1L]/density/abs(sin(angle/180 * pi)) x0 <- floor((min(x) - first.y * xd/yd)/x.shift) * x.shift y0 <- 0 x.end <- max(x) - last.y * xd/yd while (x0 < x.end) { hatch_ls[[i]] <- polygon.onehatch(x, y, x0, y0, xd, yd) i <- i+1 x0 <- x0 + x.shift } } return(do.call("rbind",hatch_ls)) }