Commit d3c16597 authored by Reinhard Furrer's avatar Reinhard Furrer

Version 1.0 as on CRAN

parent c1a13db4
^.Rhistory
^\.git/*$
^.gitignore$
^.*\.tar\.gz$
^PRSim-manual\.pdf$
^README\.md$
^.*\.Rproj$
^\.Rproj\.user$
Package: PRSim
Type: Package
Title: Stochastic Simulation of Streamflow Time Series using Phase Randomization
Version: 1.0
Date: 2019-03-10
Authors@R: c(person("Manuela", "Brunner", role = c("aut", "cre"),
email = "manuela.brunner@wsl.ch",
comment = c(ORCID = "0000-0001-8824-877X")),
person("Reinhard", "Furrer", role = c("aut"),
email = "reinhard.furrer@math.uzh.ch",
comment = c(ORCID = "0000-0002-6319-2332")))
Author: Manuela Brunner [aut, cre] (<https://orcid.org/0000-0001-8824-877X>),
Reinhard Furrer [aut] (<https://orcid.org/0000-0002-6319-2332>)
Maintainer: Manuela Brunner <manuela.brunner@wsl.ch>
Description: Provides a simulation framework to simulate streamflow time series with similar
main characteristics as observed data. These characteristics include the distribution of daily
streamflow values and their temporal correlation as expressed by short- and long-range
dependence. The approach is based on the randomization of the phases of the Fourier
transform. We further use the flexible four-parameter Kappa distribution, which allows
for the extrapolation to yet unobserved low and high flows.
URL: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
BugReports: https://git.math.uzh.ch/reinhard.furrer/PRSim-devel
License: GPL-3
Encoding: UTF-8
LazyData: true
Depends: homtest
Suggests: lattice
Imports: stats
importFrom("homtest", "Lmoments", "rand.kappa", "par.kappa")
importFrom("stats", "rnorm", "runif", "fft", "ks.test")
export("prsim")
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: knitr
LaTeX: pdfLaTeX
BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
PackageCheckArgs: --as-cran
This diff is collapsed.
# PRSim-devel
The R package `PRSim` provides a simulation framework to simulate streamflow time series with similar main characteristics as observed data. These characteristics include the distribution of daily streamflow values and their temporal correlation as expressed by short- and long-range dependence. The approach is based on the randomization of the phases of the Fourier transform. We further use the flexible four-parameter Kappa distribution, which allows for the extrapolation to yet unobserved low and high flows.
This is the development repository. A recent stable version is posted on https://CRAN.R-project.org/package=PRSim
File added
PRSim Simulating runoff
PRSim-validate Validation of simulated runoff
# short demo
data( simulations)
sim <- simulations$simulation
# periodogram of deseasonalized
kern <- kernel("modified.daniell",c(10,10))
sp1 <- spec.pgram(sim$Qobs, k=kern, taper=0, log="no", plot=FALSE)
sp2 <- spec.pgram(sim$des, k=kern, taper=0, log="no", plot=FALSE)
plot(sp1, xlim=c(0,.05))
plot( sp2, add=T, col=2)
# Peaks correspond to the following cycles:
1/sp1$freq[head(order(sp1$spec, decreasing=TRUE))]
# compare periodogram of simulated series
plot(sp1, xlim=c(0,.05)) # would be nice to identify the peaks...
for (i in grep("r",names(sim))) {
spi <- spec.pgram(sim[,i], k=kern, taper=0, log="no", plot=FALSE)
plot( spi, add=T, col="gray")
}
sp3 <- spec.pgram(sim$Qobs, taper=0, log="no", plot=FALSE)
1/sp3$freq[head(order(sp3$spec, decreasing=TRUE))]
# Annual, 6 months and 4 months
### plot mean regime for each simulation run and compare to observed regime
### 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)
year <- unique(sim$YYYY)
### compute mean runoff hydrograph
sim$day_id <- rep(seq(1:365),times=length(year))
mean_hydrograph_obs <- aggregate(sim$Qobs, by=list(sim$day_id), FUN=mean,simplify=FALSE)
plot(unlist(mean_hydrograph_obs[,2]), lty=1, lwd=1, col="black", ylab=expression(paste("Discharge [m"^3,"/s]")),
xlab="Time [d]", main="Mean hydrographs", ylim=c(0,max(unlist(mean_hydrograph_obs[,2]))*1.5),type="l")
### add mean runoff hydrographs
for(r in 7:(length(names(sim))-1)){
mean_hydrograph <- aggregate(sim[,r], by=list(sim$day_id), FUN=mean,simplify=FALSE)
lines(mean_hydrograph, lty=1, lwd=1, col=col_sim)
}
### redo observed mean
lines(mean_hydrograph_obs, lty=1, lwd=1, col="black")
### autocorrelation
acf_mare <- list()
acf_obs <- acf(sim$Qobs, plot=FALSE)
plot(acf_obs$acf, type="l", xlab="Lag", main="Autocorrelation", ylab="ACF")
for(r in 7:(length(names(sim))-2)){
acf_sim <- acf(sim[,r], plot=FALSE)
lines(acf_sim$acf, col=col_sim, type="l")
### compute mean relative error in the acf
acf_mare[[r]]<- mean(abs((acf_obs$acf-acf_sim$acf)/acf_obs$acf))
}
lines(acf_obs$acf)
### partial autocorrelation function
pacf_obs <- pacf(sim$Qobs, plot=FALSE)
pacf_mare <- list()
plot(pacf_obs$acf, type="l", xlab="Lag", main="Partial autocorrelation", ylab="PACF")
for(r in 7:(length(names(sim))-2)){
pacf_sim <- pacf(sim[,r], plot=FALSE)
lines(pacf_sim$acf, col=col_sim, type="l")
### compute mean relative error in the acf
pacf_mare[[r]] <- mean(abs((pacf_obs$acf-pacf_sim$acf)/pacf_obs$acf))
}
lines(pacf_obs$acf)
### compute seasonal statistics
### Q50,Q05,Q95, boxplots
### define seasons: Winter:12,1,2; spring:3,4,5; summer: 6,7,8; fall: 9,10,11
sim$season <- "winter"
sim$season[which(sim$MM%in%c(3,4,5))] <- "spring"
sim$season[which(sim$MM%in%c(6,7,8))] <- "summer"
sim$season[which(sim$MM%in%c(9,10,11))] <- "fall"
### all simulated series show the same seasonal statistics. plot only one
boxplot(sim$Qobs[which(sim$season=="winter")], sim$r1[which(sim$season=="winter")],
sim$Qobs[which(sim$season=="spring")], sim$r1[which(sim$season=="spring")],
sim$Qobs[which(sim$season=="summer")], sim$r1[which(sim$season=="summer")],
sim$Qobs[which(sim$season=="fall")], sim$r1[which(sim$season=="fall")],
border=c("black", col_sim, "black", col_sim, "black", col_sim, "black", col_sim),
xaxt="n", main="Seasonal statistics", outline=FALSE)
mtext(side=1, text=c("Winter", "Spring", "Summer", "Fall"), at=c(1.5,3.5,5.5,7.5))
# short demo
data(runoff)
out <- prsim(runoff, 1, 10)
sim <- out$simulation
par(mai=c(.9,.9,.1,.1))
plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))
matlines(sim$timestamp[1:1000], sim[1:1000, grep("r", names(sim))],
lty=1, col="gray")
\name{PRSim-package}
\alias{PRSim-package}
\alias{PRSim}
\docType{package}
\title{
\packageTitle{PRSim}
}
\description{
\packageDescription{PRSim}
}
\details{
The DESCRIPTION file:
\packageDESCRIPTION{PRSim}
\packageIndices{PRSim}
Simulation in the frequency domain is based on the randomization of the phases of the Fourier transform. We here combine phase randomization simulation with the flexible, four-parameter Kappa distribution, which allows for the extrapolation to yet unobserved low and high flows. The
simulation approach consists of eight steps: 1) fitting of theoretical Kappa distribution, 2) normalization and deseasonalization, 3) Fourier transformation, 4) Fourier phases computation, 5) random phase generation, 6) inverse Fourier transformation, 7) back transformation, and 8) simulation.
}
\author{
\packageAuthor{PRSim}
Maintainer: \packageMaintainer{PRSim}
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
}
\keyword{ package }
\examples{
\dontrun{
demo("PRSim")
demo("PRSim-validate")
}
}
\name{pRsim}
\alias{PRsim}
\alias{prsim}
\title{Simulate for one station}
\description{
Applies the algorithm to a single station
}
\usage{
prsim(data, station_id="Qobs", number_sim=1, win_h_length=15,
verbose=TRUE, kappapar=TRUE, suppWarn=FALSE, KStest=FALSE)
}
\arguments{
\item{data}{data frame containing the time indications and runoff of at least one station. See \sQuote{Details}.}
\item{station_id}{identifies the station in case several runoffs are present in \code{data}. See \sQuote{Details}.}
\item{win_h_length}{(half-)length of moving window size.}
\item{number_sim}{number of simulations to be carried out.}
\item{verbose}{logical. Should progress be reported?}
\item{kappapar}{logical. Should the estimated parameters of the kappa distribution be returned?}
\item{suppWarn}{logical. See \sQuote{Details}.}
\item{KStest}{logical. Should a \code{ks.test} for daily data be performed and the p-values be passed back?}
}
\details{
Time can be given with three columns named \code{"YYYY"}, \code{"MM"}, \code{"DD"}, or as in POSIXct format \code{YYYY-MM-DD}.
All leap days (Feb 29th) will be omitted from the analysis, but no missing observations are allowed.
Stations are identified by column name (default \code{"Qobs"}), or by column index.
The function \code{homtest::par.kappa} might issue quite a few warnings of type \code{In fn(par, ...) : value out of range in 'gammafn'}. The argument \code{suppWarn} allows to silence warnings for the specific function call via \code{suppressWarnings()}. Of course, a subsequent check via \code{warnings()} is recommended.
}
\value{A list with elements
\item{simulation }{A data frame with time information, observations, deseaonalized observations and
\code{number_sim} columns containing the simulated runoff.}
\item{kappa_pars}{A matrix containing the estimated parameters of the kappa distribution (if \code{kappapar}).}
\item{p_val}{A vector containing the p-values of \code{ks.test} applied to the daily detrended data (if \code{KStest}).}
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
}
\author{
Manuela Brunner
}
%\note{
%% ~~further notes~~
%}
%% ~Make other sections like Warning with \section{Warning }{....} ~
%\seealso{
%% ~~objects to See Also as \code{\link{help}}, ~~~
%}
\examples{
data( runoff)
out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
# warnings() # as a follow-up to `suppWarn=TRUE`
}
\keyword{ts}
\name{runoff}
\alias{runoff}
\docType{data}
\title{
Sample runoff of a catchment
}
\description{
Artifical runoff data based on actual and simulated observations.
}
\usage{data("runoff")}
\format{
A data frame with 15695 observations of the following 4 variables.
\describe{
\item{\code{YYYY}}{a numeric vector, year}
\item{\code{MM}}{a numeric vector, month}
\item{\code{DD}}{a numeric vector, day}
\item{\code{Qobs}}{a numeric vector, synthetic observed runoff}
}
}
\details{
The data mimiks the runoff of the river Plessur at the gauging station Chur, Switzerland. The the flow regime of the river is melt dominated. More information is given in the reference below.
}
\source{
The provided data is a weighted average of the acutually observed values and a particular simulated runoff.
The actual discharge data can be ordered from
\url{http://www.bafu.admin.ch/wasser/13462/13494/15076/index}.
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
}
\examples{
data(runoff)
str(runoff)
runoff$timestamp <- paste(runoff$YYYY, runoff$MM, runoff$DD, sep=" ")
runoff$timestamp <- as.POSIXct(strptime(runoff$timestamp,
format="\%Y \%m \%d", tz="GMT"))
plot(runoff$timestamp[1:1000], runoff$Qobs[1:1000], type="l",
xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))
}
\keyword{datasets}
\name{simulations}
\alias{simulations}
\docType{data}
\title{
Simulated runoff
}
\description{
The dataset is generated with the package own routines and represent 50 series of 18 years of runoff
}
\usage{data("simulations")}
\format{
A list of three elements, containing (i) a data frame with 6570 observations of the following 56 variables
\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{deseasonalized}}{deseasonalized time series}
\item{\code{Qobs}}{observed runoff}
\item{\code{r1},\dots,\code{r50}}{50 simulated runoff series}
}
(ii) a data frame with the daily fitted kappa parameters and (iii) p-values of the daily \code{ks.test}.
}
\details{
The data is included to illustrate the validation and visualization routines in \code{demo("PRSim-validate")}.
}
\source{
The data has been generated with
\code{set.seed(14); prsim( runoff[ runoff$YYYY>1999,], number_sim=50,
KStest=TRUE)}
(default values for all other arguments).
}
\references{
Brunner, Bardossy, Furrer (2019) Technical note: Stochastic simulation of streamflow time series using phase randomization. Submitted.
}
\examples{
data(simulations)
names(simulations)
sim <- simulations$simulation
dim(sim)
sim$day_id <- rep(seq(1:365), times=length(unique(sim$YYYY)))
mean_obs <- aggregate(sim$Qobs, by=list(sim$day_id), FUN=mean, simplify=FALSE)
plot(unlist(mean_obs[,2]),lty=1,lwd=1,col="black", ylab="Discharge [m3/s]",
xlab="Time [d]", main="Mean hydrographs", ylim=c(0,22), type="l")
for(r in 7:(length(names(sim))-1)){
mean_hydrograph <- aggregate(sim[,r], by=list(sim$day_id), FUN=mean, simplify=FALSE)
lines(mean_hydrograph, lty=1, lwd=1, col="gray")
}
lines( mean_obs, lty=1, lwd=1, col="black")
}
\keyword{datasets}
R Under development (unstable) (2019-03-24 r76263) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> pkgname <- "PRSim"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('PRSim')
Loading required package: homtest
>
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("PRSim-package")
> ### * PRSim-package
>
> flush(stderr()); flush(stdout())
>
> ### Name: PRSim-package
> ### Title: Stochastic Simulation of Streamflow Time Series using Phase
> ### Randomization
> ### Aliases: PRSim-package PRSim
> ### Keywords: package
>
> ### ** Examples
>
> ## Not run:
> ##D demo("PRSim")
> ##D demo("PRSim-validate")
> ## End(Not run)
>
>
>
> cleanEx()
> nameEx("fun_stoch_sim")
> ### * fun_stoch_sim
>
> flush(stderr()); flush(stdout())
>
> ### Name: pRsim
> ### Title: Simulate for one station
> ### Aliases: PRsim prsim
> ### Keywords: ts
>
> ### ** Examples
>
> data( runoff)
> out <- prsim( runoff[ runoff$YYYY<1980, ], "Qobs", 1, suppWarn=TRUE)
Detrending with (half-)length 15...
Starting 1 simulations:
.
Finished.
> # warnings() # as a follow-up to `suppWarn=TRUE`
>
>
>
>
> cleanEx()
> nameEx("runoff")
> ### * runoff
>
> flush(stderr()); flush(stdout())
>
> ### Name: runoff
> ### Title: Sample runoff of a catchment
> ### Aliases: runoff
> ### Keywords: datasets
>
> ### ** Examples
>
> data(runoff)
> str(runoff)
'data.frame': 15706 obs. of 4 variables:
$ YYYY: int 1975 1975 1975 1975 1975 1975 1975 1975 1975 1975 ...
$ MM : int 1 1 1 1 1 1 1 1 1 1 ...
$ DD : int 1 2 3 4 5 6 7 8 9 10 ...
$ Qobs: num 2.05 1.75 1.62 1.58 1.47 ...
> runoff$timestamp <- paste(runoff$YYYY, runoff$MM, runoff$DD, sep=" ")
> runoff$timestamp <- as.POSIXct(strptime(runoff$timestamp,
+ format="%Y %m %d", tz="GMT"))
> plot(runoff$timestamp[1:1000], runoff$Qobs[1:1000], type="l",
+ xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))
>
>
>
> cleanEx()
> nameEx("simulations")
> ### * simulations
>
> flush(stderr()); flush(stdout())
>
> ### Name: simulations
> ### Title: Simulated runoff
> ### Aliases: simulations
> ### Keywords: datasets
>
> ### ** Examples
>
> data(simulations)
> names(simulations)
[1] "simulation" "kappa_pars" "p_val"
> sim <- simulations$simulation
> dim(sim)
[1] 6570 56
> sim$day_id <- rep(seq(1:365), times=length(unique(sim$YYYY)))
> mean_obs <- aggregate(sim$Qobs, by=list(sim$day_id), FUN=mean, simplify=FALSE)
> plot(unlist(mean_obs[,2]),lty=1,lwd=1,col="black", ylab="Discharge [m3/s]",
+ xlab="Time [d]", main="Mean hydrographs", ylim=c(0,22), type="l")
>
> for(r in 7:(length(names(sim))-1)){
+ mean_hydrograph <- aggregate(sim[,r], by=list(sim$day_id), FUN=mean, simplify=FALSE)
+ lines(mean_hydrograph, lty=1, lwd=1, col="gray")
+ }
> lines( mean_obs, lty=1, lwd=1, col="black")
>
>
>
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 1.34 0.032 1.377 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')
# some simple testing commands...
# testthat would be an overshoot...
require(PRSim)
# raw demos
demo( "PRSim", ask=FALSE)
demo( "PRSim-validate", ask=FALSE)
# testing input
data(runoff)
unique(runoff$YYYY)
try( prsim( runoff[1:130,] )) # At least one year of data required.
try( prsim( runoff[1:730,] )) # At least one year of data required.
try( prsim( runoff[1:1445,] )) # No missing values allowed. Some days are missing.
try( prsim( runoff[runoff$YYYY<1976,] )) # At least one year of data required.
suppressWarnings( out <- prsim( runoff[runoff$YYYY<1977,] ) )
runof <- runoff[runoff$YYYY<1980,]
set.seed(1)
out1 <- prsim( runof, kappapar=FALSE, suppWarn=TRUE)
runo <- runof
names( runo) <- tolower( names(runof))
try( prsim( runo, kappapar=FALSE, suppWarn=TRUE)) # Wrong column for observations selected.
runo <- runof[,4:1]
set.seed(1)
out3 <- prsim( runo, kappapar=FALSE, suppWarn=TRUE) # ok
identical(out1,out3)
runo <- runof[,4:1]
set.seed(1)
out4 <- prsim( runo, station_id=1, kappapar=FALSE, suppWarn=TRUE) # ok
identical(out1,out4)
tmp <- paste(runof$YYYY, runof$MM, runof$DD,sep=" ")
runo <- data.frame(time=as.POSIXct(strptime(tmp, format="%Y %m %d", tz="GMT")), Qobs=runof$Qobs)
set.seed(1)
out5 <- prsim( runo, kappapar=FALSE, suppWarn=TRUE) # ok
identical(out1,out5)
#
\ No newline at end of file
R Under development (unstable) (2019-03-24 r76263) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> # some simple testing commands...
> # testthat would be an overshoot...
>
> require(PRSim)
Loading required package: PRSim
Loading required package: homtest
>
> # raw demos
> demo( "PRSim", ask=FALSE)
demo(PRSim)
---- ~~~~~
> # short demo
> data(runoff)
> out <- prsim(runoff, 1, 10)
Detrending with (half-)length 15...
Starting 10 simulations:
..........
Finished.
> sim <- out$simulation
> par(mai=c(.9,.9,.1,.1))
> plot(sim$timestamp[1:1000], sim$Qobs[1:1000], type="l",
+ xlab="Time [d]", ylab=expression(paste("Discharge [m"^3,"/s]")))
> matlines(sim$timestamp[1:1000], sim[1:1000, grep("r", names(sim))],
+ lty=1, col="gray")
> demo( "PRSim-validate", ask=FALSE)
demo(PRSim-validate)
---- ~~~~~~~~~~~~~~
> # short demo
> data( simulations)
> sim <- simulations$simulation
> # periodogram of deseasonalized
> kern <- kernel("modified.daniell",c(10,10))
> sp1 <- spec.pgram(sim$Qobs, k=kern, taper=0, log="no", plot=FALSE)
> sp2 <- spec.pgram(sim$des, k=kern, taper=0, log="no", plot=FALSE)
> plot(sp1, xlim=c(0,.05))
> plot( sp2, add=T, col=2)
> # Peaks correspond to the following cycles:
> 1/sp1$freq[head(order(sp1$spec, decreasing=TRUE))]
[1] 355.2632 375.0000 337.5000 321.4286 397.0588 306.8182
> # compare periodogram of simulated series
> plot(sp1, xlim=c(0,.05)) # would be nice to identify the peaks...
> for (i in grep("r",names(sim))) {
+ spi <- spec.pgram(sim[,i], k=kern, taper=0, log="no", plot=FALSE)
+ plot( spi, add=T, col="gray")
+ }
> sp3 <- spec.pgram(sim$Qobs, taper=0, log="no", plot=FALSE)
> 1/sp3$freq[head(order(sp3$spec, decreasing=TRUE))]
[1] 375.0000 355.2632 182.4324 122.7273 337.5000 421.8750
> # Annual, 6 months and 4 months
>
>
> ### plot mean regime for each simulation run and compare to observed regime
> ### 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)
> year <- unique(sim$YYYY)
> ### compute mean runoff hydrograph
> sim$day_id <- rep(seq(1:365),times=length(year))
> mean_hydrograph_obs <- aggregate(sim$Qobs, by=list(sim$day_id), FUN=mean,simplify=FALSE)
> plot(unlist(mean_hydrograph_obs[,2]), lty=1, lwd=1, col="black", ylab=expression(paste("Discharge [m"^3,"/s]")),
+ xlab="Time [d]", main="Mean hydrographs", ylim=c(0,max(unlist(mean_hydrograph_obs[,2]))*1.5),type="l")
> ### add mean runoff hydrographs
> for(r in 7:(length(names(sim))-1)){
+ mean_hydrograph <- aggregate(sim[,r], by=list(sim$day_id), FUN=mean,simplify=FALSE)
+ lines(mean_hydrograph, lty=1, lwd=1, col=col_sim)
+ }
> ### redo observed mean
> lines(mean_hydrograph_obs, lty=1, lwd=1, col="black")
> ### autocorrelation
> acf_mare <- list()
> acf_obs <- acf(sim$Qobs, plot=FALSE)
> plot(acf_obs$acf, type="l", xlab="Lag", main="Autocorrelation", ylab="ACF")
> for(r in 7:(length(names(sim))-2)){
+ acf_sim <- acf(sim[,r], plot=FALSE)
+ lines(acf_sim$acf, col=col_sim, type="l")
+ ### compute mean relative error in the acf
+ acf_mare[[r]]<- mean(abs((acf_obs$acf-acf_sim$acf)/acf_obs$acf))
+ }
> lines(acf_obs$acf)
> ### partial autocorrelation function
> pacf_obs <- pacf(sim$Qobs, plot=FALSE)
> pacf_mare <- list()
> plot(pacf_obs$acf, type="l", xlab="Lag", main="Partial autocorrelation", ylab="PACF")
> for(r in 7:(length(names(sim))-2)){
+ pacf_sim <- pacf(sim[,r], plot=FALSE)
+ lines(pacf_sim$acf, col=col_sim, type="l")
+ ### compute mean relative error in the acf
+ pacf_mare[[r]] <- mean(abs((pacf_obs$acf-pacf_sim$acf)/pacf_obs$acf))
+ }
> lines(pacf_obs$acf)
> ### compute seasonal statistics
> ### Q50,Q05,Q95, boxplots
> ### define seasons: Winter:12,1,2; spring:3,4,5; summer: 6,7,8; fall: 9,10,11
> sim$season <- "winter"
> sim$season[which(sim$MM%in%c(3,4,5))] <- "spring"
> sim$season[which(sim$MM%in%c(6,7,8))] <- "summer"
> sim$season[which(sim$MM%in%c(9,10,11))] <- "fall"
> ### all simulated series show the same seasonal statistics. plot only one
> boxplot(sim$Qobs[which(sim$season=="winter")], sim$r1[which(sim$season=="winter")],
+ sim$Qobs[which(sim$season=="spring")], sim$r1[which(sim$season=="spring")],
+ sim$Qobs[which(sim$season=="summer")], sim$r1[which(sim$season=="summer")],
+ sim$Qobs[which(sim$season=="fall")], sim$r1[which(sim$season=="fall")],
+ border=c("black", col_sim, "black", col_sim, "black", col_sim, "black", col_sim),
+ xaxt="n", main="Seasonal statistics", outline=FALSE)
> mtext(side=1, text=c("Winter", "Spring", "Summer", "Fall"), at=c(1.5,3.5,5.5,7.5))
>
>
>
> # testing input
> data(runoff)
> unique(runoff$YYYY)