Commit d1b8215b authored by Manuela's avatar Manuela

checking and running prsim.wave

documentation and validation examples
new multivariate dataset for four example catchments
parent bef00ff8
......@@ -82,7 +82,11 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### generate a day index
# data[[l]]$index <- rep(c(1:365), times=length(unique(data[[l]]$YYYY)))
data[[l]]$index <- as.numeric(format(data[[l]]$timestamp,format='%j'))
}
### replace empty index positions
if(length(which(is.na(data[[l]]$index))>0)){
data[[l]]$index[which(is.na(data[[l]]$index))] <- rep(c(1:365), times=length(unique(data[[l]]$YYYY)))[which(is.na(data[[l]]$index))]
}
}
if (verbose) cat(paste0("Detrending with (half-)length ",win_h_length,"...\n"))
......@@ -105,7 +109,9 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### TO RESOLVE
### data[[l]]$index is somehow overwritten
if(marginal=='empirical'){
marginal_list[[l]]<-'empirical'
}
if(marginal=="kappa"){
marginal_list[[l]] <- 'kappa'
p_vals <- numeric(365)
......@@ -187,6 +193,7 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
}
}
}
par_day_list[[l]] <- par_day
}
### use either a predefined distribution in R or define own function
......@@ -210,8 +217,8 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
data_random <- rCDF(n=length(data_window), theta)
# density_gengam[[d]] <- density(data_gengam)
hist(data_window)
hist(data_random,add=T,col="red")
# hist(data_window)
# hist(data_random,add=T,col="red")
if (tolower(GoFtest)=="ks"){
p_vals[d] <- ks.test(data_window,data_random)$p.value
}
......@@ -221,9 +228,11 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
### store parameters
par_day[d,] <- theta
}
par_day_list[[l]] <- par_day
}
par_day_list[[l]] <- par_day
}
### replace NA values by mean if necessary: otherwise, problems with transformation
for(l in 1:length(data)){
# if(length(which(is.na(data[[l]]$Qobs)))>0){
......@@ -235,15 +244,15 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
}
### repeat stochastic simulation several times
### list for storing results
data_sim <- list()
if(verbose) cat(paste0("Starting ",number_sim," simulations:\n"))
### run through all stations
for(l in 1:length(data)){
### list for storing results
data_sim <- list()
### simulate n series
for (r in 1:number_sim){
for (r in c(1:number_sim)){
###===============================###===============================###
### use the R-package wmtsa, which relates to the book by Percival and Walden
### on wavelet methods for time series analysis
......@@ -273,9 +282,6 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
modulus <- Mod(morlet_mat)
### extract phases (argument)
phases <- Arg(morlet_mat)
### generate r surrogates
# n_series <- 10
data_sim <- list()
### use the noise matrix corresponding to this run
noise_mat <- noise_mat_r[[r]]
......@@ -401,7 +407,7 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
data_new$simulated_seasonal[which(data_new$index%in%c(d))] <- data_ordered$cdf[data_new$rank[which(data[[l]]$index%in%c(d))]]
}
} # end for loop
data_sim[[r]] <- data_new$simulated_seasonal
data_sim[[r]] <- data_new$simulated_seasonal
if(verbose) cat(".")
### next simulation run
......@@ -430,7 +436,8 @@ prsim.wave <- function(data, station_id="Qobs", number_sim=1, win_h_length=15,
}
}else{
# return(list(simulation=data_stoch))
save(file=paste(l,'_stoch_sim.rda',sep=''),simulation=data_stoch)
out_list <- list(simulation=data_stoch, pars=NULL, p_val=NULL)
save(file=paste(l,'_stoch_sim.rda',sep=''),simulation=out_list)
}
### to to next station
}
......
# short demo
### load data for four stations
data(simulations_multi_sites)
sim <- simulations_multi_sites
### loads simulations_multi_sites
# sim <- list()
# load("~/PRSim-devel/data/simulations_multi_sites/1_stoch_sim.rda")
# sim[[1]] <- out_list$simulation
# load("~/PRSim-devel/data/simulations_multi_sites/2_stoch_sim.rda")
# sim[[2]] <- out_list$simulation
# load("~/PRSim-devel/data/simulations_multi_sites/3_stoch_sim.rda")
# sim[[3]] <- out_list$simulation
# load("~/PRSim-devel/data/simulations_multi_sites/3_stoch_sim.rda")
# sim[[4]] <- out_list$simulation
# setwd("~/PRSim-devel/data")
# save(simulations_multi_sites,file='simulations_multi_sites.rda')
### define plotting colors
col_sim <- adjustcolor("#fd8d3c",alpha=0.8)
col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)
col_obs <- adjustcolor( "black", alpha.f = 0.2)
### greys
col_vect_obs <- c('#cccccc','#969696','#636363','#252525')
### oranges
col_vect_sim <- c('#fdbe85','#fd8d3c','#e6550d','#a63603')
### plot time series for multiple sites
par(mfrow=c(2,1),mar=c(3,3,2,1))
### determine ylim
ylim_max <- max(sim[[1]]$Qobs)*1.5
### observed
plot(sim[[1]]$Qobs[1:1000],ylab=expression(bold(paste("Specific discharge [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_obs[1],ylim=c(0,ylim_max),main='Observations')
for(l in 2:4){
lines(sim[[l]]$Qobs[1:1000],col=col_vect_obs[l])
}
legend('topleft',legend=c('Station 1','Station 2','Station 3','Station 4'),lty=1,col=col_vect_obs[1:4])
### simulated (one run)
plot(sim[[1]]$r1[1:1000],ylab=expression(bold(paste("Specific discharge [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_sim[1],ylim=c(0,ylim_max),main='Stochastic simulations')
for(l in 2:4){
lines(sim[[l]]$r1[1:1000],col=col_vect_sim[l])
}
### plot cross-correlation function
par(mfrow=c(4,4),mar=c(0,0,0,0))
### run through each station comtination
for(j in 1:4){
for(i in 1:4){
### ccf of observations
data_mat <- matrix(unlist(lapply(sim, "[", , "Qobs")),ncol=length(sim))
ccf_obs <- ccf(data_mat[,i],data_mat[,j],plot=FALSE)
### plot ccfs of observations
plot(ccf_obs$lag,ccf_obs$acf,col=col_obs,type="l",ylim=c(0,1),main='',xaxt='n',yaxt='n')
### simulated ccf
### run through each simulation run
for(r in 1:2){
data_mat_sim <- matrix(unlist(lapply(sim, "[", , paste("r",r,sep=""))),ncol=length(sim))
ccf_sim <- ccf(na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,1],na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,2],plot=FALSE)
### add one ccf plot per simulation run
lines(ccf_obs$lag,ccf_sim$acf,col=col_sim)
}
### overplot observations again
lines(ccf_obs$lag,ccf_obs$acf,col="black",lwd=2)
}
}
# short demo
data(runoff_multi_sites)
### for empirical distribution
prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="empirical",out_dir="~/PRSim-devel/data/simulations_multi_sites")
### for kappa distribution
prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="kappa", GoFtest = "KS",out_dir="~/PRSim-devel/data/simulations_multi_sites")
### for GEV distribution
require("evd")
require("ismev")
rGEV <- function(n, theta) rgev(n, theta[1], theta[2], theta[3])
pGEV <- function(x, theta) pgev(x, theta[1], theta[2], theta[3])
GEV_fit <- function( xdat, ...) gev.fit( xdat, show=FALSE, ...)$mle
### GEV
prsim.wave(data=runoff_multi_sites, number_sim=1, marginal="GEV", GoFtest = "KS", n_par=3, out_dir="~/PRSim-devel/data/simulations_multi_sites")
### load stochastically simulated time series
### station 2
load("~/PRSim-devel/data/simulations_multi_sites/3_stoch_sim.rda")
sim <- out_list$simulation
par(mai=c(.9,.9,.1,.1))
### observed time series
plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))
### add simulations
matlines(sim$timestamp[1:1000], sim[1:1000, grep("r", names(sim))],
lty=1, col="gray")
\name{simulations_multi_sites}
\alias{simulations.multi.sites}
\docType{data}
\title{
Simulated runoff for four catchments
}
\description{
The dataset is generated with the package own routines and represent 10 series of 38 years of runoff for four catchments
}
\usage{data("simulations_multi_sites")}
\format{
A list of four elements (one per catchment), containing a data frame each holding information about the observed time series and the stochastic simulations
\describe{
\item{\code{YYYY}}{a numeric vector, year}
\item{\code{MM}}{a numeric vector, month}
\item{\code{DD}}{a numeric vector, day}
\item{\code{timestamp}}{\code{POSIXct} vector of the daily runoff}
\item{\code{Qobs}}{observed runoff}
\item{\code{r1},\dots,\code{r10}}{10 simulated runoff series}
}
}
\details{
The data is included to illustrate the validation and visualization routines in \code{demo("PRSim_wave-validate")}.
}
\source{
The data has been generated with
\code{prsim.wave(data=runoff_multi_sites, number_sim=10, marginal="kappa", GoFtest = "KS",out_dir="~/PRSim-devel/data/simulations_multi_sites")}
(default values for all other arguments).
}
\references{
Brunner and Gilleland to be submitted
}
\examples{
data(simulations_multi_sites)
sim <- simulations_multi_sites
dim(sim[[1]])
### plot time series for multiple sites
par(mfrow=c(2,1),mar=c(3,3,2,1))
### determine ylim
ylim_max <- max(sim[[1]]$Qobs)*1.5
### observed
plot(sim[[1]]$Qobs[1:1000],ylab=expression(bold(paste("Specific discharge [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_obs[1],ylim=c(0,ylim_max),main='Observations')
for(l in 2:4){
lines(sim[[l]]$Qobs[1:1000],col=col_vect_obs[l])
}
legend('topleft',legend=c('Station 1','Station 2','Station 3','Station 4'),lty=1,col=col_vect_obs[1:4])
### simulated (one run)
plot(sim[[1]]$r1[1:1000],ylab=expression(bold(paste("Specific discharge [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_sim[1],ylim=c(0,ylim_max),main='Stochastic simulations')
for(l in 2:4){
lines(sim[[l]]$r1[1:1000],col=col_vect_sim[l])
}
\keyword{datasets}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment