diff --git a/NAMESPACE b/NAMESPACE index a27aaaa260e57e255e7714262bc2a1be9aaa3e32..696a3e89078ea52b1a70f32e6cb00738166e76d1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -30,6 +30,7 @@ export(MultiEOF) export(PlotCombinedMap) export(PlotForecastPDF) export(PlotMostLikelyQuantileMap) +export(PlotPDFsOLE) export(PlotTriangles4Categories) export(RFSlope) export(RainFARM) diff --git a/NEWS.md b/NEWS.md index f3ebd2b3f04785e3fb3b26dacec8004912aaf1fb..208fec3928a3834ae63edfc05f78da7354981a6f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,8 @@ + CST_Anomaly allows to smooth the climatology with filter.span parameter + PlotTriangles4Categories new plotting function to convert any 3-d numerical array to a grid of coloured triangles. + CST_WeatherRegimes/WeatherRegimes and CST_RegimeAssign/RegimeAssign + + PlotPDFsOLE plots two probability density gaussian functions and the optimal linear estimation + - Fixes + CST_Anomaly handles exp, obs or both + PlotForecastPDF vignette displays figures correctly diff --git a/R/PlotPDFsOLE.R b/R/PlotPDFsOLE.R new file mode 100644 index 0000000000000000000000000000000000000000..c7239c73a3089ab8389c91f677b18767cd97371f --- /dev/null +++ b/R/PlotPDFsOLE.R @@ -0,0 +1,246 @@ +#' Plotting two probability density gaussian functions and the optimal linear +#' estimation (OLE) as result of combining them. +#' +#' @author Eroteida Sanchez-Garcia - AEMET, //email{esanchezg@aemet.es} +#' +#' @description This function plots two probability density gaussian functions +#' and the optimal linear estimation (OLE) as result of combining them. +#' +#' @param pdf_1 A numeric array with a dimension named 'statistic', containg +#' two parameters: mean' and 'standard deviation' of the first gaussian pdf +#' to combining. +#' @param pdf_2 A numeric array with a dimension named 'statistic', containg +#' two parameters: mean' and 'standard deviation' of the second gaussian pdf +#' to combining. +#' @param nsigma (optional) A numeric value for setting the limits of X axis. +#' (Default nsigma = 3). +#' @param plotfile (optional) A filename where the plot will be saved. +#' (Default: the plot is not saved). +#' @param width (optional) A numeric value indicating the plot width in +#' units ("in", "cm", or "mm"). (Default width = 30). +#' @param height (optional) A numeric value indicating the plot height. +#' (Default height = 15). +#' @param units (optional) A character value indicating the plot size +#' unit. (Default units = 'cm'). +#' @param dpi (optional) A numeric value indicating the plot resolution. +#' (Default dpi = 300). +#' +#' @return PlotPDFsOLE() returns a ggplot object containing the plot. +#' +#' @import ggplot2 +#' +#' @examples +#' # Example 1 +#' pdf_1 <- c(1.1,0.6) +#' attr(pdf_1, "name") <- "NAO1" +#' dim(pdf_1) <- c(statistic = 2) +#' pdf_2 <- c(1,0.5) +#' attr(pdf_2, "name") <- "NAO2" +#' dim(pdf_2) <- c(statistic = 2) +#' +#' PlotPDFsOLE(pdf_1, pdf_2) +#' +#'@export +PlotPDFsOLE <- function(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) { + + + if(!is.null(plotfile)){ + if (!is.numeric(dpi)) { + stop("Parameter 'dpi' must be numeric.") + } + if (length(dpi) > 1) { + warning("Parameter 'dpi' has length greater than 1 and ", + "only the first element will be used.") + dpi <- dpi[1] + } + if (!is.character(units)) { + stop("Parameter 'units' must be character") + } + if (length(units) > 1) { + warning("Parameter 'units' has length greater than 1 and ", + "only the first element will be used.") + units <- units[1] + } + if(!(units %in% c("in", "cm", "mm"))) { + stop("Parameter 'units' must be equal to 'in', 'cm' or 'mm'.") + } + if (!is.numeric(height)) { + stop("Parameter 'height' must be numeric.") + } + if (length(height) > 1) { + warning("Parameter 'height' has length greater than 1 and ", + "only the first element will be used.") + height <- height[1] + } + if (!is.numeric(width)) { + stop("Parameter 'width' must be numeric.") + } + if (length(width) > 1) { + warning("Parameter 'width' has length greater than 1 and ", + "only the first element will be used.") + width <- width[1] + } + if (!is.character(plotfile)) { + stop("Parameter 'plotfile' must be a character string ", + "indicating the path and name of output png file.") + } + } + if (!is.numeric(nsigma)) { + stop("Parameter 'nsigma' must be numeric.") + } + if (length(nsigma) > 1) { + warning("Parameter 'nsigma' has length greater than 1 and ", + "only the first element will be used.") + nsigma <- nsigma[1] + } + if (!is.array(pdf_1)) { + stop("Parameter 'pdf_1' must be an array.") + } + if (!is.array(pdf_2)) { + stop("Parameter 'pdf_2' must be an array.") + } + if (!is.numeric(pdf_1)) { + stop("Parameter 'pdf_1' must be a numeric array.") + } + if (!is.numeric(pdf_2)) { + stop("Parameter 'pdf_2' must be a numeric array.") + } + if (is.null(names(dim(pdf_1))) || + is.null(names(dim(pdf_2)))) { + stop("Parameters 'pdf_1' and 'pdf_2' ", + "should have dimmension names.") + } + if(!('statistic' %in% names(dim(pdf_1)))) { + stop("Parameter 'pdf_1' must have dimension 'statistic'.") + } + if(!('statistic' %in% names(dim(pdf_2)))) { + stop("Parameter 'pdf_2' must have dimension 'statistic'.") + } + if (length(dim(pdf_1)) != 1) { + stop("Parameter 'pdf_1' must have only dimension 'statistic'.") + } + if (length(dim(pdf_2)) != 1) { + stop("Parameter 'pdf_2' must have only dimension 'statistic'.") + } + if ((dim(pdf_1)['statistic'] != 2) || (dim(pdf_2)['statistic'] != 2)) { + stop("Length of dimension 'statistic'", + "of parameter 'pdf_1' and 'pdf_2' must be equal to 2.") + } + if(!is.null(attr(pdf_1, "name"))){ + if(!is.character(attr(pdf_1, "name"))){ + stop("The 'name' attribute of parameter 'pdf_1' must be a character ", + "indicating the name of the variable of parameter 'pdf_1'.") + } + } + if(!is.null(attr(pdf_2, "name"))){ + if(!is.character(attr(pdf_2, "name"))){ + stop("The 'name' attribute of parameter 'pdf_2' must be a character ", + "indicating the name of the variable of parameter 'pdf_2'.") + } + } + if(is.null(attr(pdf_1, "name"))){ + name1 <- "variable 1" + } else { + name1 <- attr(pdf_1, "name") + } + if(is.null(attr(pdf_2, "name"))){ + name2 <- "Variable 2" + } else { + name2 <- attr(pdf_2, "name") + } + + #----------------------------------------------------------------------------- + # Set parameters of gaussian distributions (mean and sd) + #----------------------------------------------------------------------------- + mean1 <- pdf_1[1] + sigma1 <- pdf_1[2] + mean2 <- pdf_2[1] + sigma2 <- pdf_2[2] + pdfBest <- CombinedPDFs(pdf_1, pdf_2) + meanBest <- pdfBest[1] + sigmaBest <- pdfBest[2] + + + #----------------------------------------------------------------------------- + # Plot the gaussian distributions + #----------------------------------------------------------------------------- + nameBest <- paste0(name1, " + ", name2) + graphicTitle <- "OPTIMAL LINEAR ESTIMATION" + xlimSup <- max(nsigma*sigmaBest+meanBest, nsigma*sigma1+mean1, + nsigma*sigma2+mean2) + xlimInf <- min(-nsigma*sigmaBest+meanBest, -nsigma*sigma1+mean1, + -nsigma*sigma2+mean2) + # deltax <- 0.02 + deltax <- (xlimSup - xlimInf)/10000 + + x <- seq(xlimInf,xlimSup,deltax) + df1 <- data.frame(x=x,y=dnorm(x,mean=mean1, sd=sigma1), + type=name1) + df2 <- data.frame(x=x,y=dnorm(x,mean=mean2, sd=sigma2), + type=name2) + df3 <- data.frame(x=x,y=dnorm(x,mean=meanBest, sd=sigmaBest), + type=nameBest) + df123 <- rbind(df1,df2,df3) + label1 <- paste0(name1, ": N(mean=",round(mean1,2),", sd=", round(sigma1,2), + ")") + label2 <- paste0(name2, ": N(mean=",round(mean2,2),", sd=", round(sigma2,2), + ")") + labelBest <- paste0(nameBest, ": N(mean=",round(meanBest,2),", sd=", + round(sigmaBest,2), ")") + cols <- c("#DC3912", "#13721A", "#1F5094") + names(cols) <- c(name1, name2, nameBest) + g <- ggplot(df123)+geom_line(aes(x,y,colour=type), size=rel(1.2)) + + g <- g + scale_colour_manual(values = cols, + limits = c(name1, name2, nameBest), + labels = c(label1, label2, labelBest)) + g <- g + theme(plot.title=element_text(size=rel(1.1), colour="black", + face= "bold"), + axis.text.x = element_text(size=rel(1.2)), + axis.text.y = element_text(size=rel(1.2)), + axis.title.x = element_blank(), + legend.title = element_blank(), + legend.position = c(1,1), legend.justification = c(1,1), + legend.text = element_text(face = "bold")) + g <- g + ggtitle(graphicTitle) + g <- g + labs(y="probability", size=rel(1.9)) + g <- g + stat_function(fun = dnorm_limit, args = list(mean=mean1, sd=sigma1), + fill = cols[name1], alpha=0.2, geom="area") + g <- g + stat_function(fun = dnorm_limit, args = list(mean=mean2, sd=sigma2), + fill = cols[name2], alpha=0.2, geom="area") + g <- g + stat_function(fun = dnorm_limit, args = list(mean=meanBest, + sd=sigmaBest), + fill = cols[nameBest], alpha=0.2, geom="area") + + + #----------------------------------------------------------------------------- + # Save to plotfile if needed, and return plot + #----------------------------------------------------------------------------- + if (!is.null(plotfile)) { + ggsave(plotfile, g, width = width, height = height, units = units, dpi = dpi) + } + return(g) +} + +# Auxiliar function to plot +CombinedPDFs <- function(pdf_1, pdf_2) { + mean_1 <- pdf_1[1] + sigma_1 <- pdf_1[2] + mean_2 <- pdf_2[1] + sigma_2 <- pdf_2[2] + a_1 <- (sigma_2^2)/((sigma_1^2)+(sigma_2^2)) + a_2 <- (sigma_1^2)/((sigma_1^2)+(sigma_2^2)) + pdf_mean <- a_1*mean_1 + a_2*mean_2 + pdf_sigma <- sqrt((sigma_1^2)*(sigma_2^2)/((sigma_1^2)+(sigma_2^2))) + data <- c(pdf_mean, pdf_sigma) + dim(data) <- c(statistic = 2) + return(data) +} + +dnorm_limit <- function(x,mean,sd){ + y <- dnorm(x,mean,sd) + y[x mean+sd] <- NA + return(y) +} diff --git a/man/PlotPDFsOLE.Rd b/man/PlotPDFsOLE.Rd new file mode 100644 index 0000000000000000000000000000000000000000..f2e2be8c8c75cd70dabba63accaa5aafbcc37cc3 --- /dev/null +++ b/man/PlotPDFsOLE.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/PlotPDFsOLE.R +\name{PlotPDFsOLE} +\alias{PlotPDFsOLE} +\title{Plotting two probability density gaussian functions and the optimal linear +estimation (OLE) as result of combining them.} +\usage{ +PlotPDFsOLE( + pdf_1, + pdf_2, + nsigma = 3, + plotfile = NULL, + width = 30, + height = 15, + units = "cm", + dpi = 300 +) +} +\arguments{ +\item{pdf_1}{A numeric array with a dimension named 'statistic', containg +two parameters: mean' and 'standard deviation' of the first gaussian pdf +to combining.} + +\item{pdf_2}{A numeric array with a dimension named 'statistic', containg +two parameters: mean' and 'standard deviation' of the second gaussian pdf + to combining.} + +\item{nsigma}{(optional) A numeric value for setting the limits of X axis. +(Default nsigma = 3).} + +\item{plotfile}{(optional) A filename where the plot will be saved. +(Default: the plot is not saved).} + +\item{width}{(optional) A numeric value indicating the plot width in +units ("in", "cm", or "mm"). (Default width = 30).} + +\item{height}{(optional) A numeric value indicating the plot height. +(Default height = 15).} + +\item{units}{(optional) A character value indicating the plot size +unit. (Default units = 'cm').} + +\item{dpi}{(optional) A numeric value indicating the plot resolution. +(Default dpi = 300).} +} +\value{ +PlotPDFsOLE() returns a ggplot object containing the plot. +} +\description{ +This function plots two probability density gaussian functions +and the optimal linear estimation (OLE) as result of combining them. +} +\examples{ +# Example 1 +pdf_1 <- c(1.1,0.6) +attr(pdf_1, "name") <- "NAO1" +dim(pdf_1) <- c(statistic = 2) +pdf_2 <- c(1,0.5) +attr(pdf_2, "name") <- "NAO2" +dim(pdf_2) <- c(statistic = 2) + +PlotPDFsOLE(pdf_1, pdf_2) + +} +\author{ +Eroteida Sanchez-Garcia - AEMET, //email{esanchezg@aemet.es} +} diff --git a/tests/testthat/test-PlotPDFsOLE.R b/tests/testthat/test-PlotPDFsOLE.R new file mode 100644 index 0000000000000000000000000000000000000000..bd1279d9df403c795910d5d9b921f00b6df38592 --- /dev/null +++ b/tests/testthat/test-PlotPDFsOLE.R @@ -0,0 +1,121 @@ +context("Generic tests") + +test_that("Sanity checks", { + pdf_1 <- c(1.1,0.6) + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(statistic = 2) + pdf_2 <- c(1,0.5) + attr(pdf_2, "name") <- "NAO2" + dim(pdf_2) <- c(statistic = 2) + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = "plot.png", + width = 30, height = 15, + units = "cm", dpi = '300') , + "Parameter 'dpi' must be numeric.") + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = "plot.png", + width = 30, height = 15, + units = 20, dpi = 300) , + "Parameter 'units' must be character") + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = "plot.png", + width = 30, height = 15, + units = "dm", dpi = 300) , + "Parameter 'units' must be equal to 'in', 'cm' or 'mm'.") + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = "plot.png", + width = 30, height = '15', + units = "cm", dpi = 300) , + "Parameter 'height' must be numeric.") + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = "plot.png", + width = list(30), height = 15, + units = "cm", dpi = 300) , + "Parameter 'width' must be numeric.") + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = 0, + width = 30, height = 15, + units = "cm", dpi = 300) , + paste0("Parameter 'plotfile' must be a character string ", + "indicating the path and name of output png file.")) + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = '3', plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + "Parameter 'nsigma' must be numeric.") + + pdf_1 <- list(1.1,0.6) + attr(pdf_1, "name") <- "NAO1" + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + "Parameter 'pdf_1' must be an array.") + + pdf_1 <- c('1.1','0.6') + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(statistic = 2) + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + "Parameter 'pdf_1' must be a numeric array.") + + pdf_1 <- c(1.1,0.6) + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(2) + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + paste0("Parameters 'pdf_1' and 'pdf_2' ", + "should have dimmension names.")) + + pdf_1 <- c(1.1,0.6) + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(statisti = 2) + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + "Parameter 'pdf_1' must have dimension 'statistic'.") + + pdf_1 <- c(1.1,0.6) + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(statistic = 2, model = 1) + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + "Parameter 'pdf_1' must have only dimension 'statistic'.") + + pdf_1 <- c(1.1, 0.6, 0.2) + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(statistic = 3) + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + paste0("Length of dimension 'statistic'", + "of parameter 'pdf_1' and 'pdf_2' must be equal to 2.")) + + pdf_1 <- c(1.1, 0.6) + attr(pdf_1, "name") <- 12 + dim(pdf_1) <- c(statistic = 2) + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + paste0("The 'name' attribute of parameter 'pdf_1' must be a character ", + "indicating the name of the variable of parameter 'pdf_1'.")) + + pdf_1 <- c(1.1,0.6) + attr(pdf_1, "name") <- "NAO1" + dim(pdf_1) <- c(statistic = 2) + pdf_2 <- c(1,0.5) + attr(pdf_2, "name") <- 12 + dim(pdf_2) <- c(statistic = 2) + + expect_error(PlotPDFsOLE(pdf_1, pdf_2, nsigma = 3, plotfile = NULL, + width = 30, height = 15, + units = "cm", dpi = 300) , + paste0("The 'name' attribute of parameter 'pdf_2' must be a character ", + "indicating the name of the variable of parameter 'pdf_2'.")) + + + + +}) diff --git a/vignettes/BestEstimateIndex_vignette.Rmd b/vignettes/BestEstimateIndex_vignette.Rmd index e577eeb2cffbfe417333d123d76b48c3f52d8179..e549533f389fec6b3eb074a22cdf9d61109cfdf3 100644 --- a/vignettes/BestEstimateIndex_vignette.Rmd +++ b/vignettes/BestEstimateIndex_vignette.Rmd @@ -1,134 +1,174 @@ ---- -author: "Eroteida Sánchez-García" -date: "`r Sys.Date()`" -output: rmarkdown::html_vignette -vignette: > - %\VignetteEngine{knitr::knitr} - %\VignetteIndexEntry{Achiving Best Estimate Index} - %\usepackage[utf8]{inputenc} ---- - -Achiving the Precipitation Best prediction giving the NAO index -------------------------------------------------------------------- - -The boreal winter precipitation forecast, acumulated from November to March, can be improved by considering the NAO index. The first step is to find the best estimation of winter NAO, given by two Seasonal Forecast System (SFS). The second step is to employ the enhanced NAO index pdf to produce weights for a SFS (it could be the same or a different SFS from the previous one). The third step is to apply this weights to a precipitation field. The methodology has been proved to improve the skill on precipitation forecast in the iberian peninsula given the relation between the winter precipitation and the NAO index in seasonal time scale (Sánchez-García, E., Voces-Aboy, J., Navascués, N., & Rodríguez-Camino, E. (2019). Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Adv. Sci. Res., 16, 165174, ). - -This document is aim to ilustrate a practical use of the functions included in CSTools applying this methodology. - - -## Loading packages and data - -Open an R sesion and load library CSTools: - -``` -library(CSTools) -``` -The required data to applied this methodology are: -- the winter index NAO for two (or three) different SFS in a reference period (hindcast) and in a future simulation (forecast), -- the observed (reconstructed) NAO index in the reference period and -- the acumulated precipitation field from a SFS that aims to be improved (hindcast or hindcast and forecast) - -Given the memory limitations, the following example uses synthetic data. - -The first SFS is a dynamical model containing 25 ensemble members and its output will be save in the object `NAO_hind1` for the 20 years reference period and `NAO_fcst1` for a next season forecast. -The second SFS is a empirical model, so, the `NAO_hind2` and `NAO_fcst2` are characterized by a mean and standard deviation saved in the 'statistics' dimension. - -The synthetic data is created by running the following lines: - - -``` -NAO_obs <- rnorm(3, sd=3) -dim(NAO_obs) <- c(time = 3) -NAO_hind1 <- rnorm(3 * 6, sd=3) -dim(NAO_hind1) <- c(time = 3, member = 6) -NAO_hind2_mean <- rnorm(3, sd=3) -NAO_hind2_sd <- rnorm(3, mean=5, sd=1) -NAO_hind2 <- cbind(NAO_hind2_mean,NAO_hind2_sd) -dim(NAO_hind2) <- c(time=3, statistic=2) - -``` - -The acumulated precipiation field (only march for memory limiations) is loaded by running: - -``` -lonlat_prec$data <- apply(lonlat_prec$data, c(1, 2, 3, 5, 6), sum) -``` - - -## 1- Best Estimate Index NAO - -The function `PDFBest` performes the next steps: - - improves the PDF NAO index for each SFS (it could use bias correction method) and - - does a statistical combination of these improved NAO indexes. -Its output is an array containing the parameters (mean and starndar deviation) of the PDF for the reference period. - - -``` -pdf_hind_best <- BEI_PDFBest(NAO_obs, NAO_hind1, NAO_hind2, index_fcst1 = NULL, - index_fcst2 = NULL, method_BC = 'none', - time_dim_name = 'time', na.rm = FALSE) -``` - - -## 2- Compute weights using the Best Estimation of Index NAO - - -An array of weights is calculated for a SFS. This SFS could be a different SFS than the ones in section 1. -The function WeightIndex computes these weights for each ensemble's member based on the best NAO PDF estimate. - - -``` -weights_hind <- BEI_Weights(NAO_hind1, pdf_hind_best) -names(dim(weights_hind))[1] <- 'sdate' -``` - -The expected dimensions of these weights are 'member' and temporal dimension. - - -## 3- Apply weights to a precipitation field - - -The function `CST_BEI_Weighting` computes the ensemble mean or the terciles probabilities for a climate variable. -The corrected precipitation field is obtained by running: - -``` -em_prec <- CST_BEI_Weighting(lonlat_prec, weights_hind, type = 'ensembleMean', time_dim_name = 'sdate') -prob_prec <- CST_BEI_Weighting(lonlat_prec, weights_hind, type = 'probs', ) -``` - -## Comparison and visualization - -The original model output can be compared against the BEI corrected field. -To do this a equiprobability weights are created and they are applied to the original precipitation field: - - -``` -aweights_raw <- rep(1/6, 3 * 6) -dim(aweights_raw) <- dim(weights_hind) - -em_prec_raw <- CST_BEI_Weighting(lonlat_prec, aweights_raw, type = 'ensembleMean', time_dim_name = 'sdate') -prob_prec_raw <- CST_BEI_Weighting(lonlat_prec, aweights_raw, type = 'probs', time_dim_name = 'sdate') - -``` - -A map with the probability that the total precipitation will be in the lower/normal/upper tercile based on the Best Estimate Index NAO could be obtained using the 'PlotEquiMap' function or 'PlotMostLikelyQuantileMap' function from 'CSToools' package. - -The following figures show the probabilities for lower tercile for precipitation from November to March 2012/13 for ECMWF S5 system applying the methodology exposed or not, obtained using real data: -- NAO indices from the ECMWF-S5 dynamical model and the S-ClimWaRe empirical model from AEMET from 1997 to 2016 for computing the Best Estimation of Index NAO fro this hindcast period. -- The winter precipitation (from November to March) from 1997 to 2016 over Iberia Peninsula from he ECMWF-S5 dynamical model with resolution 0.5º x 0.5º, to weighting with the previous Best Estimation of Index NAO. - - -![](./Figures/BestEstimateIndex_fig1.png){width=70%} - -![](./Figures/BestEstimateIndex_fig3.png){width=70%} - - - -In a similar way we can plot the map with the probability that the total precipitation from November 2012 to -March 2013, for example, will be in the lower tercile from ECMWF Seasonal Forecast System 5 (raw) to compare results running the code: - - - -![](./Figures/BestEstimateIndex_fig2.png){width=70%} - +--- +author: "Eroteida Sánchez-García" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{Achiving Best Estimate Index} + %\usepackage[utf8]{inputenc} +--- + +Achiving the Precipitation Best prediction giving the NAO index +------------------------------------------------------------------- + +The boreal winter precipitation forecast, acumulated from November to March, can be improved by considering the NAO index. The first step is to find the best estimation of winter NAO, given by two Seasonal Forecast System (SFS). The second step is to employ the enhanced NAO index pdf to produce weights for a SFS (it could be the same or a different SFS from the previous one). The third step is to apply this weights to a precipitation field. The methodology has been proved to improve the skill on precipitation forecast in the iberian peninsula given the relation between the winter precipitation and the NAO index in seasonal time scale (Sánchez-García, E., Voces-Aboy, J., Navascués, N., & Rodríguez-Camino, E. (2019). Regionally improved seasonal forecast of precipitation through Best estimation of winter NAO, Adv. Sci. Res., 16, 165174, ). + +This document is aim to ilustrate a practical use of the functions included in CSTools applying this methodology. + + +## Loading packages and data + +Open an R sesion and load library CSTools: + +``` +library(CSTools) +``` +The required data to applied this methodology are: +- the observed (reconstructed) NAO index in the reference period +- the winter index NAO for two different SFSs (SFS1 and SFS2, to combine both of them) in a reference period (hindcast) and in a future simulation (forecast) +- the winter index NAO and the acumulated precipitation field from a SFS that aims to be improved (hindcast and forecast) + +Given the memory limitations, the following example uses synthetic data. + +The SFS1 system is a dynamical model containing 25 ensemble members and its output will be save in the object `NAO_hind1` for the 20 years reference period and `NAO_fcst1` for a next season forecast. +The second SFS, SFS2, is a empirical model, so, the `NAO_hind2` and `NAO_fcst2` are characterized by a mean and standard deviation saved in the 'statistics' dimension. +The model for improving is a dynamical model containing 25 ensemble members. +The synthetic data is created by running the following lines: + + +``` +# observations +NAO_obs <- rnorm(20, sd=3) +dim(NAO_obs) <- c(time = 20) + +# hindcast and forecast of a dynamical SFS 1 +NAO_hind1 <- rnorm(20 * 2 * 25, sd=2.5) +dim(NAO_hind1) <- c(time = 20, member = 50) +NAO_fcst1 <- rnorm(2*51, sd=2.5) +dim(NAO_fcst1) <- c(time = 1, member = 102) + +# hindcast and forecast of an empirical SFS 2 +NAO_hind2_mean <- rnorm(20, sd=3) +NAO_hind2_sd <- rnorm(20, mean=5, sd=1) +NAO_hind2 <- cbind(NAO_hind2_mean, NAO_hind2_sd) +dim(NAO_hind2) <- c(time=20, statistic=2) +NAO_fcst2_mean <- rnorm(1, sd=3) +NAO_fcst2_sd <- rnorm(1, mean=5, sd=1) +NAO_fcst2 <- cbind(NAO_fcst2_mean, NAO_fcst2_sd) +dim(NAO_fcst2) <- c(time=1, statistic=2) + +``` + +The winter index NAO and the acumulated precipiation field from the dynamical SFS that aims to be improved could be created by running: + +``` +# NAO index of a SFS to compute weights for each ensemble's member +NAO_hind <- rnorm(20 * 25, sd=2.5) +dim(NAO_hind) <- c(time = 20, member = 25) +NAO_fcst <- rnorm(51, sd=2.5) +dim(NAO_fcst) <- c(time = 1, member = 51) + +# The acumulated precipiation field +prec_hind <- rnorm(20 * 25 * 21 * 31, mean=30, sd=10) +dim(prec_hind) <- c(time = 20, member = 25, lat = 21, lon = 31) +prec_hind <- list(data = prec_hind) +class(prec_hind) <- 's2dv_cube' +prec_fcst <- rnorm(51 * 21 * 31, mean=25,sd=8) +dim(prec_fcst) <- c(time = 1, member = 51, lat = 21, lon = 31) +prec_fcst <- list(data = prec_fcst) +class(prec_fcst) <- 's2dv_cube' + +``` + + +## 1- Best Estimate Index NAO + +The function `PDFBest` performes the next steps: + - improves the PDF NAO index for each SFS (it could use bias correction method) and + - does a statistical combination of these improved NAO indexes. +Its output is an array containing the parameters (mean and starndar deviation) of the PDF for the reference period (hindcast) or forecast period. + + +``` +# for hindcast +pdf_hind_best <- BEI_PDFBest(NAO_obs, NAO_hind1, NAO_hind2, index_fcst1 = NULL, + index_fcst2 = NULL, method_BC = 'none', + time_dim_name = 'time', na.rm = FALSE) +# for forecast +pdf_fcst_best <- BEI_PDFBest (NAO_obs, NAO_hind1, NAO_hind2, index_fcst1 = NAO_fcst1, + index_fcst2 = NAO_fcst2, method_BC = 'none', + time_dim_name = 'time', na.rm = FALSE) +``` + + +## 2- Compute weights using the Best Estimation of Index NAO + + +An array of weights is calculated for the SFS. This SFS could be the same or different SFS than the ones in section 1. +The function WeightIndex computes these weights for each ensemble's member based on the best NAO PDF estimate. + + +``` +# for hindcast +weights_hind <- BEI_Weights(NAO_hind, pdf_hind_best) +# for forecast +weights_fcst <- BEI_Weights(NAO_fcst, pdf_fcst_best) + +``` + +The expected dimensions of these weights are 'member' and temporal dimension. + + +## 3- Apply weights to a precipitation field + + +The function `CST_BEI_Weighting` computes the ensemble mean or the terciles probabilities for a climate variable. +The ensemble mean and the probabilities of terciles from the weighted precipitation field is obtained by running: + +``` +# for hindcast +em_prec_hind <- CST_BEI_Weighting(prec_hind, weights_hind, type = 'ensembleMean') +prob_prec_hind <- CST_BEI_Weighting(prec_hind, weights_hind, type = 'probs') + +# for forecast +em_prec_fcst <- CST_BEI_Weighting(prec_fcst, weights_fcst, type = 'ensembleMean') +prob_prec_fcst <- CST_BEI_Weighting(prec_fcst, weights_fcst, type = 'probs') + +``` + +## Comparison and visualization + +The original model output can be compared against the BEI corrected field. +To do this a equiprobability weights are created for hindcast period and they are applied to the original precipitation field: + + +``` +aweights_raw <- rep(1/25, 20 * 25) +dim(aweights_raw) <- dim(weights_hind) + +em_prec_raw <- CST_BEI_Weighting(prec_hind, aweights_raw, type = 'ensembleMean') +prob_prec_raw <- CST_BEI_Weighting(prec_hind, aweights_raw, type = 'probs') + +``` + +A map with the probability that the total precipitation will be in the lower/normal/upper tercile based on the Best Estimate Index NAO could be obtained using the 'PlotEquiMap' function or 'PlotMostLikelyQuantileMap' function from 'CSToools' package. + +The following figures show the probabilities for lower tercile for precipitation from November to March 2012/13 for ECMWF S5 system applying the methodology exposed or not, obtained using real data: +- NAO indices from the ECMWF-S5 dynamical model and the S-ClimWaRe empirical model from AEMET from 1997 to 2016 for computing the Best Estimation of Index NAO fro this hindcast period. +- The winter precipitation (from November to March) from 1997 to 2016 over Iberia Peninsula from he ECMWF-S5 dynamical model with resolution 0.5º x 0.5º, to weighting with the previous Best Estimation of Index NAO. + + +![](./Figures/BestEstimateIndex_fig1.png){width=70%} + +![](./Figures/BestEstimateIndex_fig3.png){width=70%} + + + +In a similar way we can plot the map with the probability that the total precipitation from November 2012 to +March 2013, for example, will be in the lower tercile from ECMWF Seasonal Forecast System 5 (raw) to compare results running the code: + + + +![](./Figures/BestEstimateIndex_fig2.png){width=70%} + ![](./Figures/BestEstimateIndex_fig4.png){width=70%} \ No newline at end of file