NULL
#'
#' The comprehensive Precipitation Generator
#'
#' @param station character vector of the IDs of the considered meteorological stations
#' @param prec_all data frame containing daily precipitation of all meteorological stations. See \code{\link{PRECIPITATION}} defined in the \code{\link{trentino}} dataset for formatting.
#' @param mean_climate_prec a matrix containing monthly mean daily precipitation for the considered station. If it is \code{NULL}, it is calculated. See input of \code{\link{is.monthly.climate}}
#' @param year_max start year of the recorded (calibration) period
#' @param year_min end year of the recorded (calibration) period
#' @param leap logical variables. If it is \code{TRUE} (default)(recommended), leap years are considered, otherwise all years have 365 days
#' @param nmonth number of months in one year (default is 12)
#' @param verbose logical variable
#' @param cpf see \code{\link{normalizeGaussian_severalstations}}
#' @param sample,extremes,qnull,valmin see \code{\link{normalizeGaussian_severalstations}}
#' @param step see \code{\link{normalizeGaussian_severalstations}}. Default is 0.
#' @param p,type,lag.max,ic,activateVARselect see respective input parameter on \code{\link{getVARmodel}}
#' @param year_max_sim last year of the simulation period. Default is equal to \code{year_max}
#' @param year_min_sim first year of the simulation period. Default is equal to \code{year_min}
#' @param mean_climate_prec_sim a matrix containing monthly mean daily precipitation for the simulation period. If is \code{NULL} (Default), it is set equal to \code{mean_climate_prec}.
#' @param n_GPCA_iteration number of iterations of Gaussianization process for data. Default is 0 (no Gaussianization)
#' @param n_GPCA_iteration_residuals number of iterations of Gaussianization process for VAR residuals. Default is 0 (no Gaussianization)
#' @param exogen data frame or matrix containing the (normalized or not) exogenous variables (predictors) for the recorded (calibration) period.
#' @param exogen_sim data frame or matrix containing the (normalized or not) exogenous variables (predictors) for the simulation period. Default is \code{NULL}. If it is \code{NULL}, it is replaced with \code{exogen} within the function.
#' @param is_exogen_gaussian logical value. If \code{TRUE}, \code{exogen_sim} and \code{exogen} are given as already normalized variables, otherwhise they are not normalized. Default is \code{FALSE}
#' @param onlygeneration logical value. If \code{TRUE} the VAR model \code{varmodel} is given as input and only random generation is done, otherwise (default) is calculated from measured data
#' @param varmodel the comprehensinve VAR model as a \code{\link{varest2}} S4 object or a \code{NULL} object. If \code{NULL} (default), the comprehensinve VAR is estimated from measured data within the function, otherwise it is given as input and only random generation is done.
#' @param type_quantile see \code{type} on \code{\link{quantile}}
#### @param option_temp integer value. If 0 (default), exogenous variables (\code{exogen} and \code{exogen_sim} ) are not considered, if 1 exogenous variables are considered.
#' @param exogen_all data frame containing exogenous variable formatted like \code{prec_all}. Default is \code{NULL}.
#' It is alternative to \code{exogen} and if it not \code{NULL},\code{is_exogen_gaussian} is automatically set \code{FALSE}
#' @param exogen_all_col vector of considered columns of \code{exogen_all}. Default is \code{station}.
#' @param no_spline logical value. See \code{\link{splineInterpolateMonthlytoDailyforSeveralYears}}. Default is \code{TRUE}.
#' @param nscenario number of generated scenarios for daily maximum and minimum temperature
#' @param seed seed for stochastic random generation see \code{\link{set.seed}}.
#' @param noise stochastic noise to add for variabile generation. Default is \code{NULL}. See \code{\link{newVARmultieventRealization}}. Not used in case that \code{nscenario>1}.
#'
#'
#'
#' @export
#'
#'
#' @author Emanuele Cordano, Emanuele Eccel
#'
#' @seealso \code{\link{splineInterpolateMonthlytoDailyforSeveralYears}}
#'
#'
#'
#' @note It pre-processes and generates a multi-site precipitation fields. It uses \code{\link{getVARmodel}}. Detailed examples can be viewed of this function in \href{https://docs.google.com/file/d/0B8xDtMCnW3dJU2JIemVqMnpKTHc/edit}{this presentation}.
#' Unfortunately, using this approach, the spatial correlations are underestimated. This is due to the persinstence of zeros in the precipitation records.
#' This problem is known in literature and can be solved in the future versions of RMAWGEN.
#' See the R code for further details
#'
#'
#'
#'
#' @return A list of the following variables:
#'
#'
#' \code{prec_mes} matrix containing measured daily precipitation (the data is copied by the measured data given as input for the period and the station considered for \code{varmodel} estimation)
#'
#' \code{prec_spline} matrix containing climatic "spline-interpolated" daily preciptation from \code{mean_climate_prec}
#'
#' \code{data_prec} matrix containing normalized measured precipitation variable
#'
#' \code{prec_gen} matrix containing generated daily precipitation [mm]
#'
#' \code{prec_spline_sim} matrix containing climatic "spline-interpolated" daily preciptation from \code{mean_climate_prec_sim}
#'
#' \code{data_prec_gen} matrix containing normalized generated precipitation variable
#'
#' \code{mean_climate_prec} matrix containing monthly means of daily precipitation (historical scenario)
#'
#' \code{mean_climate_prec_sim} matrix containing monthly means of daily precipitation (predicted/simulated scenario)
#'
#' \code{var} a varest object containing the used VAR model
#'
#'
#'
#'
#' @examples
#'
#' data(trentino)
#' set.seed(1222) # set the seed for random generations!
#'year_max <- 1990
#'year_min <- 1961
#'year_max_sim <- 1982
#'year_min_sim <- 1981
#'
#'n_GPCA_iter <- 2
#' p <- 1
#'nscenario=1
#'station <- c("T0090","T0083")
#'\donttest{
#' ## Not Run: the call to ComprehensivePrecipitationGenerator may elapse too
#' ## long time (more than 5 eseconds) and is not executed by default CRAN check.
#' ## Please uncomment the following line to run the example on your own PC.
#' generation00 <- ComprehensivePrecipitationGenerator(station=station,
#' prec_all=PRECIPITATION,year_min=year_min,year_max=year_max,
#' year_min_sim=year_min_sim,year_max_sim=year_max_sim,p=p,
#' n_GPCA_iteration=n_GPCA_iter,n_GPCA_iteration_residuals=0,
#' sample="monthly",nscenario=nscenario,no_spline=TRUE)
#'
#' }
#'
ComprehensivePrecipitationGenerator <-
function(
station=c("T0001","T0010","T0099"),
prec_all,
mean_climate_prec=NULL,
year_max=1990,
year_min=1961,
leap=TRUE,
nmonth=12,
cpf=NULL,
verbose=TRUE,
p=1,
type="none",
lag.max=NULL,
ic="AIC",
activateVARselect=FALSE,
exogen=NULL,
exogen_sim=NULL,
is_exogen_gaussian=FALSE,
year_max_sim=year_max,
year_min_sim=year_min,
mean_climate_prec_sim=NULL,
onlygeneration=FALSE,
varmodel=NULL,
type_quantile=3,
qnull=NULL,
valmin=0.5,
step=0,
n_GPCA_iteration=0,
n_GPCA_iteration_residuals=n_GPCA_iteration,
sample=NULL,
extremes=TRUE,
exogen_all=NULL,
exogen_all_col=station,
### option_temp=0,
no_spline=FALSE,
nscenario=1,
seed=NULL,
noise=NULL
) {
useVAR=TRUE
origin <- paste(year_min,"1","1",sep="/") # Must start from Jan 1
origin_sim <- paste(year_min_sim,"1","1",sep="/") # Must start from Jan 1
prec_mes <- as.data.frame(extractyears(prec_all,year_min=year_min,year_max=year_max,station=station))
nyear <- year_max-year_min+1
if (!is.monthly.climate(mean_climate_prec,nstation=length(station),nmonth=nmonth,verbose=verbose)) mean_climate_prec <- getMonthlyMean(prec_all,year_min=year_min,year_max=year_max,station=station)
MEAN_CLIMATE_prec_SAVED <- mean_climate_prec
# TO DO SOME MODIFICATIONS
prec_spline <- as.data.frame(splineInterpolateMonthlytoDailyforSeveralYears(val=mean_climate_prec,start_year=year_min,nyear=nyear,leap=leap,no_spline=no_spline))
names(prec_spline) <- names(mean_climate_prec)
if (min(prec_spline)<=0) {
print("Error in Precipitation generator: spline interpolation of daily precipitation less than 0!")
return(-1)
}
data_prec <- normalizeGaussian_severalstations(x=prec_mes,data=prec_mes,sample=sample,cpf=cpf,step=step,origin_x=origin,origin_data=origin,extremes=extremes)
# data_prec <- normalizeGaussian_severalstations_prec(x=prec_mes,data=prec_mes,sample=sample,cpf=cpf,qnull=qnull,valmin=valmin,origin_x=origin,origin_data=origin,extremes=extremes)
if (!onlygeneration){
if (!is.null(exogen_all)) {
exogen <- as.data.frame(extractyears(exogen_all,year_min=year_min,year_max=year_max,station=exogen_all_col))
is_exogen_gaussian=FALSE
}
if (is.null(exogen_sim)) exogen_sim <- exogen
if (!is_exogen_gaussian){
exogen0 <- exogen
if (!is.null(exogen)) exogen <- normalizeGaussian_severalstations(x=exogen0,data=exogen0,sample=sample,cpf=cpf,origin_x=origin,origin_data=origin)
}
var <- getVARmodel(data=data_prec,suffix=NULL,sep="",p=p,type=type,exogen=exogen,lag.max=lag.max,ic=ic,activateVARselect=activateVARselect,n_GPCA_iteration_residuals=n_GPCA_iteration_residuals,n_GPCA_iteration=n_GPCA_iteration,extremes=extremes)
if (activateVARselect) return(var)
} else {
var <- varmodel
}
# THIS LINE IS TEMPORARY AND ONLY FOR TESTING
if (!is.null(noise)) {if (noise=="residuals") noise <- residuals(var) }
if (is.null(mean_climate_prec_sim)) mean_climate_prec_sim <- mean_climate_prec
nyear_sim <- year_max_sim-year_min_sim+1
nyear_max <- max(nyear_sim,nyear)
## print(nyear_sim)
## print(nyear_max)
prec_spline_sim <- as.data.frame(splineInterpolateMonthlytoDailyforSeveralYears(val=mean_climate_prec_sim,start_year=year_min_sim,nyear=nyear_max,leap=leap,no_spline=no_spline))
prec_spline_sim2 <- as.data.frame(splineInterpolateMonthlytoDailyforSeveralYears(val=mean_climate_prec_sim,start_year=year_min_sim,nyear=nyear_sim,leap=leap,no_spline=no_spline))# This is a spline interpotation only for the years on which generatiion is requested
names(prec_spline_sim) <- colnames(mean_climate_prec_sim)
names(prec_spline_sim2) <- colnames(mean_climate_prec_sim)
if (is.null(exogen_sim)) exogen_sim <- exogen
if (!is.null(exogen_sim) & (!is_exogen_gaussian)) {
exogen_sim0 <- exogen_sim
exogen_sim <- normalizeGaussian_severalstations(x=exogen_sim0,data=exogen_sim0,sample=sample,cpf=cpf,origin_x=origin_sim,origin_data=origin_sim,extremes=extremes)
}
if (!is.null(seed)) set.seed(seed)
data_prec_gen <- newVARmultieventRealization(var,exogen=exogen_sim,nrealization=nrow(prec_spline_sim2),extremes=extremes,type=type_quantile,noise=noise)
# data_prec_gen <- newVARmultieventRealization(var,exogen=exogen_sim,nrealization=nyear_sim,extremes=extremes,type=type_quantile)
precrows <- 1:(min(c(nrow(prec_mes),nrow(prec_spline),nrow(prec_spline_sim))))
prec_mes_rescaled <- prec_mes[precrows,]/prec_spline[precrows,]*prec_spline_sim[precrows,]
prec_gen <- as.data.frame(normalizeGaussian_severalstations(x=data_prec_gen,data=prec_mes_rescaled,inverse=TRUE,type=type_quantile,step=step,sample=sample,origin_x=origin_sim,origin_data=origin,extremes=extremes))
# data_prec <- normalizeGaussian_severalstations_prec(x=prec_mes,data=prec_mes,sample=sample,cpf=cpf,qnull=qnull,valmin=valmin,origin_x=origin,origin_data=origin,extremes=extremes)
names(prec_gen) <- names(prec_spline_sim)
colnames(data_prec_gen) <- names(prec_spline_sim) # added on may 4 2012
# names(prec_gen) <- names(prec_mes)
out <- NULL
if (onlygeneration) {
names_out <- c("prec_gen","prec_spline_sim","data_prec_gen","mean_climate_prec_sim","prec_mes","prec_spline","prec_mes_rescaled")
for (it in names_out) { if (!exists(it)) assign(it,NULL)}
out <- list(prec_gen,prec_spline_sim,data_prec_gen,mean_climate_prec_sim,prec_mes,prec_spline,prec_mes_rescaled)# May 2012
# out <- list(prec_gen,prec_spline_sim,data_prec,prec_gen,prec_spline_sim,data_prec_gen,mean_climate_prec,mean_climate_prec_sim,var)
names(out) <- names_out
} else {
names_out <- c("prec_mes","prec_spline","data_prec","prec_gen","prec_spline_sim","data_prec_gen","mean_climate_prec","mean_climate_prec_sim","var")
for (it in names_out) { if (!exists(it)) assign(it,NULL)}
out <- list(prec_mes,prec_spline,data_prec,prec_gen,prec_spline_sim,data_prec_gen,mean_climate_prec,mean_climate_prec_sim,var)
names(out) <- names_out
}
if (nscenario>1) {
for (kk in 1:nscenario) {
### 20200727
data_prec_gen <- newVARmultieventRealization(var,exogen=exogen_sim,nrealization=nrow(prec_spline_sim2),extremes=extremes,type=type_quantile)
colnames(data_prec_gen) <- names(prec_spline_sim2)### 20200727 end
# ntn_rows <- 1:nrow()
# prec_mes_rescaled <- prec_mes/prec_spline*prec_spline_sim
precrows <- 1:(min(c(nrow(prec_mes),nrow(prec_spline),nrow(prec_spline_sim))))
prec_mes_rescaled <- prec_mes[precrows,]/prec_spline[precrows,]*prec_spline_sim[precrows,]
prec_gen <- as.data.frame(normalizeGaussian_severalstations(x=data_prec_gen,data=prec_mes_rescaled,inverse=TRUE,type=type_quantile,step=step,sample=sample,origin_x=origin_sim,origin_data=origin,extremes=extremes))
# prec_gen <- as.data.frame(normalizeGaussian_severalstations_prec(x=data_prec_gen,data=prec_mes_rescaled,inverse=TRUE,type=type_quantile,qnull=qnull,valmin=valmin,sample=sample,origin_x=origin_sim,origin_data=origin,extremes=extremes))
names(prec_gen) <- names(prec_spline_sim)
prec_index <- sprintf("prec_gen%05d",kk)
out[[prec_index]] <- prec_gen
}
}
return(out)
# return(list(prec_mes=prec_mes,prec_spline=prec_spline,data_prec=data_prec,prec_gen=prec_gen,prec_spline_sim=prec_spline_sim,data_prec_gen=data_prec_gen,mean_climate_prec=mean_climate_prec,mean_climate_prec_sim=mean_climate_prec_sim,var=var))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.