Commit d7e36a8e authored by Henning Rust's avatar Henning Rust
Browse files

HR and KK did some small changes: a) add ... to plotting functions, b) add ......

HR and KK did some small changes: a) add ... to plotting functions, b) add ... to pass arguments to control of optim c) small others
parent ee66eadd
Package: IDF Package: IDF
Type: Package Type: Package
Title: Estimation and Plotting of IDF Curves Title: Estimation and Plotting of IDF Curves
Version: 1.2 Version: 1.3
Date: 2018-01-22 Date: 2018-02-06
Author: Christoph Ritschel, Carola Detring, Sarah Joedicke Author: Christoph Ritschel, Carola Detring, Sarah Joedicke
Maintainer: Christoph Ritschel <christoph.ritschel@met.fu-berlin.de> 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. 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.
......
...@@ -313,7 +313,7 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) { ...@@ -313,7 +313,7 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
#' @return $par vector of IDF parameters at optimization minimum #' @return $par vector of IDF parameters at optimization minimum
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de} #' @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) { 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 use.log=use.log
...@@ -349,7 +349,7 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F, ...@@ -349,7 +349,7 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations ## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations
## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization ## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization
fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG), fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),
control=list(trace=0,maxit=5000), control=list(...),
method=method,upper=upper,lower=lower), error=function(e) NULL)#, method=method,upper=upper,lower=lower), error=function(e) NULL)#,
#upper=upper,lower=lower) #upper=upper,lower=lower)
...@@ -358,7 +358,7 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F, ...@@ -358,7 +358,7 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations ## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations
## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization ## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization
fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG), fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),
control=list(trace=0,maxit=5000), control=list(...),
method=method), error=function(e) NULL)#, method=method), error=function(e) NULL)#,
#upper=upper,lower=lower) #upper=upper,lower=lower)
...@@ -572,7 +572,7 @@ IDF.init <- function(int.vec,durs,n.y,method="Nelder-Mead") { ...@@ -572,7 +572,7 @@ IDF.init <- function(int.vec,durs,n.y,method="Nelder-Mead") {
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, 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, 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") { 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 ### ### Calculate extreme values for each year and each aggregation level ###
...@@ -600,7 +600,7 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL ...@@ -600,7 +600,7 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
###################################################### ######################################################
if(!is.na(mu.init) | !is.na(sigma.init) | !is.na(xi.init) | !is.na(eta.init)) { 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, 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) DEBUG=DEBUG,method=method,upper=upper,lower=lower,...)
}else { }else {
cat("Warning: Optimization not carried out due to invalid initial values. \n") cat("Warning: Optimization not carried out due to invalid initial values. \n")
fit.min <- NA fit.min <- NA
...@@ -673,7 +673,8 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL ...@@ -673,7 +673,8 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta.init=0,eta.init=NA, IDF.short <- function(int.vec,durs,n.y,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, 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") { 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="Station",data.name="obs",...) {
################################################################################### ###################################################################################
### Estimate Parameters for single duration if not given initial values by user ### ### Estimate Parameters for single duration if not given initial values by user ###
...@@ -692,7 +693,7 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta ...@@ -692,7 +693,7 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta
###################################################### ######################################################
if(!is.na(mu.init) | !is.na(sigma.init) | !is.na(xi.init) | !is.na(eta.init)) { 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, 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) DEBUG=DEBUG,method=method,upper=upper,lower=lower,...)
}else { }else {
cat("Warning: Optimization not carried out due to invalid initial values. \n") cat("Warning: Optimization not carried out due to invalid initial values. \n")
fit.min <- NA fit.min <- NA
...@@ -744,7 +745,10 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta ...@@ -744,7 +745,10 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta
#' IDF.plot(pars=param,st.name="example",dt.name="rgamma") #' IDF.plot(pars=param,st.name="example",dt.name="rgamma")
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de} #' @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) { 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=rainbow(length(probs)),lty=1,
st.name="Station",dt.name="obs",ints=NA,ds=NA,ylim=NA,add=FALSE,...) {
## initialize array for IDF values at different durations and for different probabilities ## initialize array for IDF values at different durations and for different probabilities
idf.array <- array(NA,dim=c(length(dur),length(probs))) idf.array <- array(NA,dim=c(length(dur),length(probs)))
...@@ -756,24 +760,25 @@ 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 ...@@ -756,24 +760,25 @@ 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
idf.array[,i] <- qgev.d(probs[i],mu=pars[1],sigma=pars[2],xi=pars[3],theta=pars[4],eta=pars[5],d=dur) 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 } ## end of loop over probs
if(!add){
## initiialize plot window with limits of IDF values ## 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") y.range <- ifelse(is.na(ylim), c(min(idf.array[,1],na.rm=T),max(idf.array[,3],na.rm=T)),ylim)
plot(NA,axes=F,xlim=c(min(dur,na.rm=T),max(dur,na.rm=T)),ylim=y.range,xlab="duration [h]",ylab="intensity [mm/h]",log="xy",...)
axis(1,at=dur,labels=dur) axis(1,at=dur,labels=dur)
axis(2) axis(2)
points(ds,ints,pch=16,col=rgb(0,0,0,0.5)) points(ds,ints,pch=16,col=rgb(0,0,0,0.5))
## loop over probabilities ## loop over probabilities
## plot IDF curve ## 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" legend.text.2 <- "quantile"
## plot legend ## plot legend
legend(x="topright",legend=c(st.name,dt.name,paste(probs,legend.text.2,sep=" ")), 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)))) 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))))
}
for(i in 1:length(probs))
lines(dur,idf.array[,i],col=cols[i],lwd=1.5,lty=lty)
} ## end of function IDF.plot } ## end of function IDF.plot
################################################################################### ###################################################################################
......
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