Commit 9838eed6 authored by Christoph Ritschel's avatar Christoph Ritschel
Browse files

Version_1.2

parent 096be90e
Package: IDF
Type: Package
Title: Estimation and Plotting of IDF Curves
Version: 1.1
Date: 2017-09-29
Version: 1.2
Date: 2018-01-22
Author: Christoph Ritschel, Carola Detring, Sarah Joedicke
Maintainer: Christoph Ritschel <christoph.ritschel@met.fu-berlin.de>
Description: Intensity-duration-frequency (IDF) curves are a widely used analysis-tool in hydrology to assess extreme values of precipitation [e.g. Mailhot et al., 2007, <doi:10.1016/j.jhydrol.2007.09.019>]. The package 'IDF' provides a function to read precipitation data from German weather service (DWD) 'webwerdis' <http://www.dwd.de/EN/ourservices/webwerdis/webwerdis.html> files and Berlin station data from 'Stadtmessnetz' <http://www.geo.fu-berlin.de/en/met/service/stadtmessnetz/index.html> files, and additionally IDF parameters can be estimated also from a given data.frame containing a precipitation time series. The data is aggregated to given levels yearly intensity maxima are calculated either for the whole year or given months. From these intensity maxima IDF parameters are estimated on the basis of a duration-dependent generalised extreme value distribution [Koutsoyannis et al., 1998, <doi:10.1016/S0022-1694(98)00097-3>]. IDF curves based on these estimated parameters can be plotted.
......
export(IDF.read,IDF.fit,IDF.plot,TS.acc)
export(IDF.read,IDF,IDF.plot,TS.acc)
import(ismev,
stats4,
......
......@@ -80,9 +80,9 @@ IDF.read <- function(file,type){
colnames(Tab_end) <- c("year","mon","day","hour","min","RR")
attr(Tab_end,"accumulation time (min)") <- as.numeric(difftime(new_timect[2],new_timect[1], units="mins"))
# Liste <- list(t1=Tab_end)
# Liste <- list(t1=Tab_end)
Liste <- Tab_end
if (type == "webwerdis"){
# WEBWERDIS:
attr(Liste,"StationName") <- as.character(Tab$Stationname[1])
......@@ -135,25 +135,25 @@ TS.acc <- function(x,acc.val=2,moving.sum="FALSE") {
if(acc.val<1) cat(paste("Warning: accumulation value acc.val too small for accumulation of the time series \n"))
if(moving.sum){
x.acc <- as.numeric(filter(x,filter=rep(1,acc.val),method="convolution",sides=1))
}else{
l.new <- length(x)%/%acc.val ## calculate new length of accumulated time series
l.rest <- length(x)%%acc.val ## calculate values left over
if(l.rest==0) {
x.acc <- apply(matrix(x,nrow=l.new,byrow=T),1,sum)
}else{
x.acc <- apply(matrix(x[1:(length(x)-l.rest)],nrow=l.new,byrow=T),1,sum)
#cat(paste("Warning: ",l.rest,"time steps left and not used for accumulation \n"))
}
l.new <- length(x)%/%acc.val ## calculate new length of accumulated time series
l.rest <- length(x)%%acc.val ## calculate values left over
if(l.rest==0) {
x.acc <- apply(matrix(x,nrow=l.new,byrow=T),1,sum)
}else{
x.acc <- apply(matrix(x[1:(length(x)-l.rest)],nrow=l.new,byrow=T),1,sum)
#cat(paste("Warning: ",l.rest,"time steps left and not used for accumulation \n"))
}
}
## return accumulated time series
return(x.acc)
} # End of function TS.acc
#####################################################################################
......@@ -230,7 +230,7 @@ rgev.d <- function(n,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1) {
x[which(is.nan(x))] <- NA
return(x)
}
}
#######################################################################
#' @title Negativ log-likelihood of modified GEV
......@@ -264,21 +264,21 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
if(DEBUG) debug.values <- c(mu,sigma,xi,theta,eta)
if(sum(is.nan(sigma.d))==0) {
## Weibull und Frechet
if(xi!=0){
C <- 1 + xi * (x/sigma.d - mu )
nll <- switch((sum(C<0,na.rm=T)>0)+1,
sum(log(sigma.d),na.rm=T)+(1+1/xi)*sum(log(C),na.rm=T)+sum((C)^(-1/xi),na.rm=T),
NA)
# + penalty*(sum(C[C<0]^2))
## Gumbel
}else if(xi==0){# & sigma<1 & eta<1)
Y <- x/sigma.d-mu
nll <- -(-sum(log(sigma.d),na.rm=T)-sum((Y),na.rm=T)-sum(exp(-Y),na.rm=T))
}
}else{ nll <- NA}
## Weibull und Frechet
if(xi!=0){
C <- 1 + xi * (x/sigma.d - mu )
nll <- switch((sum(C<0,na.rm=T)>0)+1,
sum(log(sigma.d),na.rm=T)+(1+1/xi)*sum(log(C),na.rm=T)+sum((C)^(-1/xi),na.rm=T),
NA)
# + penalty*(sum(C[C<0]^2))
## Gumbel
}else if(xi==0){# & sigma<1 & eta<1)
Y <- x/sigma.d-mu
nll <- -(-sum(log(sigma.d),na.rm=T)-sum((Y),na.rm=T)-sum(exp(-Y),na.rm=T))
}
}else{ nll <- NA}
if(DEBUG){
cat(debug.values,nll,"\n")
options(digits.secs=6)
......@@ -289,7 +289,7 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
return(nll/length(x))
} # end of function IDF.nll
} # end of function IDF.nll
######################################################################################################
#' @title Fitting function to optimize IDF model parameters
......@@ -314,7 +314,7 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,method="Nelder-Mead",upper=Inf,lower=-Inf) {
use.log=use.log
if(use.log) {
......@@ -326,13 +326,13 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
eta <- log(eta)
if(method=="L-BFGS-B") {
upper[2] <- log(upper[2])
upper[4] <- log(upper[4])
upper[5] <- log(upper[5])
lower[2] <- log(lower[2])
lower[4] <- log(lower[4])
lower[5] <- log(lower[5])
upper[2] <- log(upper[2])
upper[4] <- log(upper[4])
upper[5] <- log(upper[5])
lower[2] <- log(lower[2])
lower[4] <- log(lower[4])
lower[5] <- log(lower[5])
}
}
......@@ -393,49 +393,140 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
} ## end of function fit.fun
##################################################################################
#' @title Data aggregation for IDF parameter estimation
#' @description The function \code{IDF.agg} aggregates a data.frame of observations \code{data} with temporal inforamtion (at least years) and values of precipitation
#' at a given temporal resoultion at given aggregation levels \code{agg.lev} and yearly maxima of intensity are caluclated for a specific month or the whole year/dataset.
#' @param data a \code{data,frame}, preferably generated by function \code{IDF.read}. It should at least contain a \code{$RR} and \code{$year} element for the
#' function tow work properly. Also an option to use \code{moving.sum} is given. The function returns a vector of intensities and durations as well as the number of years of data.
#' @param agg.lev a vector of aggregation levels used to fit the IDF curves.
#' @param month \code{integer} value specifying the month to be used for estimating the IDF parameters. Type "all" for all months or if
#' the whole time series should be fitted.
#' @param moving.sum \code{logical} specifying if moving sum filtering should be applied for time series aggregation.
#' @return $ints.vec vector of sorted intensities for selected aggregation levels
#' @return $durs vector of sorted aggregation levels
#' @return $n.y number of years of data
#' @examples
#' RR <- rgamma(10*30*24,shape=1)
#' year <- sort(rep(1:(10),30*24))
#' data <- data.frame(RR,year)
#' data.agg <- IDF.agg(data)
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.agg <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FALSE",DEBUG=FALSE) {
RR <- data$RR ## get precipitation time series from data.frame
years <- unique(data$year) # get years from data.frame
n.y <- length(years) # number of years
n.a <- length(agg.lev) # number of aggregation times
## initilise arrays
agg.1 <- array(NA,dim=c(n.y))
ints <- array(NA,dim=c(n.y*n.a))
###loop over years
for(y in 1:n.y) {
RR <- data$RR ## get precipitation time series from data.frame
years <- unique(data$year) # get years from data.frame
n.y <- length(years) # number of years
n.a <- length(agg.lev) # number of aggregation times
## initilise arrays
agg.1 <- array(NA,dim=c(n.y))
ints <- array(NA,dim=c(n.y*n.a))
###loop over years
for(y in 1:n.y) {
if(month[1]=="all") {
index <- which(data$year==years[y])
}else if(is.integer(month) | is.numeric(month)) {
index <- which(data$year==years[y] & data$mon >= min(month) & data$mon <= max(month))
}
if(length(index)>0) {
RR.year <- RR[index]
agg.1[y] <- max(RR.year,na.rm=T)
###loop over agg.lev
for(a in 1:n.a) {
if(month[1]=="all") {
index <- which(data$year==years[y])
}else if(is.integer(month) | is.numeric(month)) {
index <- which(data$year==years[y] & data$mon >= min(month) & data$mon <= max(month))
}
if(length(index)>0) {
RR.year <- RR[index]
agg.1[y] <- max(RR.year,na.rm=T)
###loop over agg.lev
for(a in 1:n.a) {
ints[y+((a-1)*n.y)] <- max(TS.acc(RR.year,agg.lev[a],moving.sum=moving.sum),na.rm=T)/agg.lev[a]
} # end for all aggregation times
} # end if lenght
} # end for all years
## vector of all intensities
int.vec <- c(agg.1,ints)
## vector of all durations (single)
d.all <- c(1,agg.lev)
## long vector of all durations (repeated for each year to have same length as intensity vector)
durs <- rep(d.all,each=n.y)
return(list(int.vec=int.vec,durs=durs,n.y=n.y))
}
ints[y+((a-1)*n.y)] <- max(TS.acc(RR.year,agg.lev[a],moving.sum=moving.sum),na.rm=T)/agg.lev[a]
} # end for all aggregation times
} # end if lenght
} # end for all years
## vector of all intensities
int.vec <- c(agg.1,ints)
## vector of all durations (single)
d.all <- c(1,agg.lev)
## long vector of all durations (repeated for each year to have same length as intensity vector)
durs <- rep(d.all,each=n.y)
return(list(int.vec=int.vec,durs=durs,n.y=n.y))
} #
###############################################################################
#' @title Estimation of initial values for IDF fitting.
#' @description The function \code{IDF.init} estimates inital values for \code{mu,sigma,xi and eta} assuming \code{theta}
#' equals zero. A generalized extreme value distribution is fitted individually for each year and then the inital values
#' for the duration dependent gev fit are estimated from those by applying a linear regression to the scale parameters of each year.
#' @param int.vec a \code{vector} of yearly maxima of intensity sorted by year and aggregatin level
#' @param durs a \code{vector} of durations used to fit the model. Has to have same length and order as \code{int.vec}
#' @param n.y \code{integer} value specifying the number of years of data.
#' @param method \code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e
#' @return $mu initial estimation of location parameter
#' @return $sigma initial estimation of scale parameter
#' @return $xi inital estimation of shape parameter
#' @return $eta intial estimation of slope parameter for sigma-power law.
#' @examples
#' RR <- rgamma(10*30*24,shape=1)
#' year <- sort(rep(1:(10),30*24))
#' data <- data.frame(RR,year)
#' data.agg <- IDF.agg(data,agg.lev=c(2,6,12,24))
#' pars.init <- IDF.init(data.agg$int.vec,data.agg$durs,data.agg$n.y)
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.init <- function(int.vec,durs,n.y,method="Nelder-Mead") {
## Fit a generalized extreme value distribution to the maximum intensities of each year for a single
## aggregation level and write the estimated parameters in an array for further analyisis.
d.all <- unique(durs)
ints.all <- matrix(int.vec,nrow=n.y) ## sort intensities in a matrix, rows are years, columns are aggregation levels
pars <- array(NA,dim=c(3,length(d.all)))
## In case of NA values the optimization fails, therefore years with NA values need to be removed.
ints.all <- matrix(ints.all[rowSums(!is.na(ints.all)) == length(d.all)],ncol=length(d.all))
if(nrow(ints.all)<3) {
cat("Warning: optimization did not converge and no parameters were estimated. Time Series contains less than 3 years of valid data. \n")
mu=NA
sigma=NA
xi=NA
eta=NA
}else{
## loop over all aggregation levels
for(d in 1:length(d.all)) {
#fit <- fit.fun.emp(obs=ints.all[,d],mu=mu,sigma=sigma,xi=xi,use.log=use.log,
# DEBUG=DEBUG,method=method,upper=upper,lower=lower)
fit <- gev.fit(xdat=ints.all[,d],method=method,show=F)
pars[,d] <- fit$mle
} ## end loop over aggregation levels
#############################################################
### Derive starting parameters for duration-dependent GEV ###
#############################################################
## Fit a linear model to the individual sigmas for individual aggregation times in a log-log environment
## The slope coefficient is an estimate for the slope in the duration-dependent GEV, namely parameter eta
## The intersection is an estimation of the starting parameter sigma
## Parameter mu is estimated as mean value of individual mus divided by indiviudal sigmas
## The initial value for xi will be the mean of all individual xi, since it is approximately independent of duration
formel <- lm(log(pars[2,]) ~ log(d.all))
sigma <- as.numeric(exp(formel$coefficients[1]))
mu <- mean(pars[1,]/pars[2,])
eta <- as.numeric(-formel$coefficients[2])
xi <- max(0,mean(pars[3,],na.rm=T))
}
return(list("mu"=mu,"sigma"=sigma,"xi"=xi,"eta"=eta))
} # EOF
#################################################################################
#' @title Fitting IDF model parameters to observations at different durations
#' @description The function \code{IDF.fit} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
......@@ -451,6 +542,10 @@ IDF.agg <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=
#' @param month \code{integer} value specifying the month to be used for estimating the IDF parameters. Type "all" for all months or if
#' the whole time series should be fitted.
#' @param moving.sum \code{logical} specifying if moving sum filtering should be applied for time series aggregation.
#' @param mu.init initial estimation of location parameter, default is NA. Initial value estimated by fitting individual gev parameters
#' @param sigma.init initial estimation of scale parameter,default is NA. Initial value estimated by fitting individual gev parameters
#' @param xi.init inital estimation of shape parameter, default is NA. Initial value estimated by fitting individual gev parameters
#' @param eta.init intial estimation of slope parameter for sigma-power law, default is NA. Initial value estimated by fitting individual gev parameters
#' @param theta.init inital value defining the curvature of the IDF, default is zero, it is not recommended to change it
#' @param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
#' @param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
......@@ -475,151 +570,119 @@ IDF.agg <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=
#' pars <- fit$par
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.fit <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FALSE",theta.init=0,
use.log=FALSE,DEBUG=FALSE,method="Nelder-Mead",upper=Inf,lower=-Inf,plot=FALSE,
probs=c(0.5,0.9,0.99),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),station.name="Berlin",data.name="obs") {
IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FALSE",mu.init=NA,sigma.init=NA,xi.init=NA,theta.init=0,eta.init=NA,
use.log=FALSE,DEBUG=FALSE,method="Nelder-Mead",upper=Inf,lower=-Inf,plot=FALSE,
probs=c(0.5,0.9,0.99),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),station.name="Berlin",data.name="obs") {
#########################################################################
### Calculate extreme values for each year and each aggregation level ###
#########################################################################
dummy.list <- IDF.agg(data,agg.lev,month,moving.sum,DEBUG=FALSE)
int.vec <- dummy.list$int.vec
durs <- dummy.list$durs
n.y <- dummy.list$n.y
###############################################
### Estimate Parameters for single duration ###
###############################################
## Fit a generalized extreme value distribution to the maximum intensities of each year for a single
## aggregation level and write the estimated parameters in an array for further analyisis.
ints.all <- matrix(int.vec,nrow=n.y) ## sort intensities in a matrix, rows are years, columns are aggregation levels
pars <- array(NA,dim=c(3,length(d.all)))
## In case of NA values the optimization fails, therefore years with NA values need to be removed.
ints.all <- matrix(ints.all[rowSums(!is.na(ints.all)) == length(d.all)],ncol=length(d.all))
if(length(ints.all)==0) {
cat("Warning: optimization did not converge and no parameters were estimated. Time Series contains too many NA values. \n")
}else{
## loop over all aggregation levels
for(d in 1:length(d.all)) {
#fit <- fit.fun.emp(obs=ints.all[,d],mu=mu,sigma=sigma,xi=xi,use.log=use.log,
# DEBUG=DEBUG,method=method,upper=upper,lower=lower)
if(method=="L-BFGS-B"){
fit <- gev.fit(xdat=ints.all[,d],method=method,upper=upper,lower=lower,show=F)
}else{
fit <- gev.fit(xdat=ints.all[,d],method=method,show=F)
}
dummy.list <- IDF.agg(data,agg.lev,month,moving.sum,DEBUG=FALSE)
int.vec <- dummy.list$int.vec
durs <- dummy.list$durs
n.y <- dummy.list$n.y
d.all <- unique(durs)
###################################################################################
### Estimate Parameters for single duration if not given initial values by user ###
###################################################################################
if(is.na(mu.init) | is.na(sigma.init) | is.na(xi.init) | is.na(eta.init)) {
pars[,d] <- fit$mle
pars.init <- IDF.init(int.vec,durs,n.y,method)
mu.init <- pars.init$mu
sigma.init <- pars.init$sigma
xi.init <- pars.init$xi
eta.init <- pars.init$eta
} ## end loop over aggregation levels
#############################################################
### Derive starting parameters for duration-dependent GEV ###
#############################################################
## Fit a linear model to the individual sigmas for individual aggregation times in a log-log environment
## The slope coefficient is an estimate for the slope in the duration-dependent GEV, namely parameter eta
## The intersection is an estimation of the starting parameter sigma
## Parameter mu is estimated as mean value of individual mus divided by indiviudal sigmas
## The initial value for xi will be the mean of all individual xi, since it is approximately independent of duration
formel <- lm(log(pars[2,]) ~ log(d.all))
sigma.est <- as.numeric(exp(formel$coefficients[1]))
mu.est <- mean(pars[1,]/pars[2,])
eta.est <- as.numeric(-formel$coefficients[2])
xi.est <- max(0,mean(pars[3,],na.rm=T))
}
######################################################
### Estimate parameters for duration-dependent GEV ###
######################################################
fit <- fit.fun(obs=int.vec,dur=durs,mu=mu.est,sigma=sigma.est,xi=xi.est,theta=theta.init,eta=eta.est,use.log=use.log,
DEBUG=DEBUG,method=method,upper=upper,lower=lower)
######################################################
### success? Than plot! ###
######################################################
if(plot&& !is.na(fit$min)) {
ds <- sort(rep(d.all,length(int.vec)/length(d.all)))
IDF.plot(pars=fit$par,probs,st.name=station.name,dt.name=data.name,ints=int.vec,ds=durs)
}
}
if(!plot && is.na(fit$min)) {
cat("Warning: optimization did not converge and no parameters were estimated. \n")
}
if(plot && is.na(fit$min)) {
cat("Warning: optimization did not converge and no parameters were estimated. Plot not possible. \n")
if(!is.na(mu.init) | !is.na(sigma.init) | !is.na(xi.init) | !is.na(eta.init)) {
fit <- fit.fun(obs=int.vec,dur=durs,mu=mu.init,sigma=sigma.init,xi=xi.init,theta=theta.init,eta=eta.init,use.log=use.log,
DEBUG=DEBUG,method=method,upper=upper,lower=lower)
}else {
cat("Warning: Optimization not carried out due to invalid initial values. \n")
fit.min <- NA
}
return(list("ints"=int.vec,"durs"=durs,"min"=fit$min,"par"=fit$par))
} ## End of function IDF.fit
######################################################################################################################
########################################################################################################
#' @title Plotting IDF curves
#' @description The function \code{IDF.plot} plots a set of IDF curves with given IDF model parameters \code{pars} for
#' several probability levels \code{probs} at given durations \code{dur}. The colors of the curves can be defined with
#' parameter \code{cols} (need to have same length as \code{probs}). The \code{station.name} will be printed in the legend.
#' @param pars a vector of IDF model parameters mu,sigma,xi,eta,theta
#' @param probs a vector of probabilities for which the IDF curves are calculated
#' @param dur a vector of durations at which the IDF curves are calculated
#' @param cols a vector of colors for the seperate IDF curves, needs same length as \code{probs}
#' @param st.name \code{character} overall naming of the IDF plot, e.g. name of location or model name
#' @param dt.name \code{character} naming the data points, e.g. obs or model name
#' @param ints \code{vector} of observational intensities (surted by durations)
#' @param ds \code{vector} of durations (same length as intensities)
#' @examples
#' RR <- rgamma(10*30*24,shape=1)
#' year <- sort(rep(1:(10),30*24))
#' data <- data.frame(RR,year)
#' fit <- IDF.fit(data)
#' param <- fit$par
#' IDF.plot(pars=param,st.name="example",dt.name="rgamma")
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.plot <- function(pars,probs=c(0.5,0.9,0.99),dur=c(0.5,1,2,3,6,12,24,48,72,96),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),st.name="Berlin-Dahlem",dt.name="obs",ints=NA,ds=NA) {
## initialize array for IDF values at different durations and for different probabilities
idf.array <- array(NA,dim=c(length(dur),length(probs)))
## loop over probabilities
for(i in 1:length(probs)) {
######################################################
### success? Than plot! ###
######################################################
## calculate IDF values for given probability at all durations
idf.array[,i] <- qgev.d(probs[i],mu=pars[1],sigma=pars[2],xi=pars[3],theta=pars[4],eta=pars[5],d=dur)
if(plot&& !is.na(fit$min)) {
ds <- sort(rep(d.all,length(int.vec)/length(d.all)))
IDF.plot(pars=fit$par,probs,st.name=station.name,dt.name=data.name,ints=int.vec,ds=durs)
}
} ## end of loop over probs
## initiialize plot window with limits of IDF values
plot(NA,axes=F,xlim=c(min(dur,na.rm=T),max(dur,na.rm=T)),ylim=c(min(idf.array[,1],na.rm=T),max(idf.array[,3],na.rm=T)),xlab="duration [h]",ylab="intensity [mm/h]",log="xy")
axis(1,at=dur,labels=dur)
axis(2)
points(ds,ints,pch=16,col=rgb(0,0,0,0.5))
## loop over probabilities
## plot IDF curve
for(i in 1:length(probs)) {
points(dur,idf.array[,i],type="l",col=cols[i],lwd=1.5)
}
legend.text.2 <- "quantile"
if(!plot && is.na(fit$min)) {
cat("Warning: optimization did not converge and no parameters were estimated. \n")
}
if(plot && is.na(fit$min)) {
cat("Warning: optimization did not converge and no parameters were estimated. Plot not possible. \n")
}
return(list("ints"=int.vec,"durs"=durs,"min"=fit$min,"par"=fit$par))
} ## End of function IDF.fit
######################################################################################################################
########################################################################################################
#' @title Plotting IDF curves
#' @description The function \code{IDF.plot} plots a set of IDF curves with given IDF model parameters \code{pars} for
#' several probability levels \code{probs} at given durations \code{dur}. The colors of the curves can be defined with
#' parameter \code{cols} (need to have same length as \code{probs}). The \code{station.name} will be printed in the legend.
#' @param pars a vector of IDF model parameters mu,sigma,xi,eta,theta
#' @param probs a vector of probabilities for which the IDF curves are calculated
#' @param dur a vector of durations at which the IDF curves are calculated
#' @param cols a vector of colors for the seperate IDF curves, needs same length as \code{probs}
#' @param st.name \code{character} overall naming of the IDF plot, e.g. name of location or model name
#' @param dt.name \code{character} naming the data points, e.g. obs or model name
#' @param ints \code{vector} of observational intensities (surted by durations)
#' @param ds \code{vector} of durations (same length as intensities)
#' @examples
#' RR <- rgamma(10*30*24,shape=1)
#' year <- sort(rep(1:(10),30*24))
#' data <- data.frame(RR,year)
#' fit <- IDF.fit(data)
#' param <- fit$par
#' IDF.plot(pars=param,st.name="example",dt.name="rgamma")
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.plot <- function(pars,probs=c(0.5,0.9,0.99),dur=c(0.5,1,2,3,6,12,24,48,72,96),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),st.name="Berlin-Dahlem",dt.name="obs",ints=NA,ds=NA) {
## initialize array for IDF values at different durations and for different probabilities
idf.array <- array(NA,dim=c(length(dur),length(probs)))
## loop over probabilities
for(i in 1:length(probs)) {
## calculate IDF values for given probability at all durations
idf.array[,i] <- qgev.d(probs[i],mu=pars[1],sigma=pars[2],xi=pars[3],theta=pars[4],eta=pars[5],d=dur)
} ## end of loop over probs
## initiialize plot window with limits of IDF values
plot(NA,axes=F,xlim=c(min(dur,na.rm=T),max(dur,na.rm=T)),ylim=c(min(idf.array[,1],na.rm=T),max(idf.array[,3],na.rm=T)),xlab="duration [h]",ylab="intensity [mm/h]",log="xy")
axis(1,at=dur,labels=dur)
axis(2)
points(ds,ints,pch=16,col=rgb(0,0,0,0.5))
## loop over probabilities
## plot IDF curve
for(i in 1:length(probs)) {
points(dur,idf.array[,i],type="l",col=cols[i],lwd=1.5)
}
legend.text.2 <- "quantile"
## plot legend
legend(x="topright",legend=c(st.name,dt.name,paste(probs,legend.text.2,sep=" ")),
col=c(1,rgb(0,0,0,0.5),cols),lty=c(NA,NA,rep(1,length(cols))),pch=c(NA,16,rep(NA,length(cols))))
} ## end of function IDF.plot
###################################################################################
## plot legend
legend(x="topright",legend=c(st.name,dt.name,paste(probs,legend.text.2,sep=" ")),
col=c(1,rgb(0,0,0,0.5),cols),lty=c(NA,NA,rep(1,length(cols))),pch=c(NA,16,rep(NA,length(cols))))
} ## end of function IDF.plot
###################################################################################