Skip to content

GitLab

Commit beaf7cb8 authored by Christoph Ritschel's avatar Christoph Ritschel
Browse files

included IDF.short in version 1.2

parent 9838eed6
Hide whitespace changes
Inline Side-by-side
Showing with 250 additions and 77 deletions
R/IDF.R
View file @ beaf7cb8
......@@ -529,7 +529,7 @@ IDF.init <- function(int.vec,durs,n.y,method="Nelder-Mead") {
#################################################################################
#' @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}
#' @description The function \code{IDF} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
#' to a data.frame of observations \code{data} with temporal inforamtion (at least years) and values of precipitation
#' at a given temporal resoultion. This precipitation time series gets aggregated at given aggregation levels.
#' \code{agg.lev} and yearly maxima of intensity are caluclated for a specific month or the whole year/dataset.
......@@ -605,84 +605,177 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
cat("Warning: Optimization not carried out due to invalid initial values. \n")
fit.min <- NA
}
######################################################
### 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")
}
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))
######################################################
### 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")
}
return(list("ints"=int.vec,"durs"=durs,"min"=fit$min,"par"=fit$par))
} ## End of function IDF.fit
######################################################################################################################
####
#' @title Fitting IDF model parameters to annual maximum intensity time series
#' @description The function \code{IDF.short} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
#' to vectors of annnual maximum intensities \code{int.vec} at different durations \code{durs}.
#' The starting values of the IDF model parameters can be determined by the user as well as specific options to use
#' during optimization. Logartihmic transformation, debugging, the optimization method, and an option to plot the
#' IDF curves.
#' @param ints.vec a \code{vector} of yearly maxima of intensity sorted by year and aggregatin level
#' @param durs a vector of aggregation levels used to fit the IDF curves. One value for each year. Has to have same lenght as \code{int.vec}
#' @param n.y \code{integer} value specifying the number of years of data
#' @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
# 'log-likelihood are printed on console for each iteration during optimization.
#' @param method \code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e
#' @param lower \code{vector} specifying the lower boundary of parameters for "L-BFGS-B" method
#' @param upper \code{vector} specifying the upper boundary of parameters for "L-BFGS-B" method
#' @param plot \code{logical} option of creating a plot of IDF curves with estimated parameters.
#' @param probs a vector of probabilities for which the IDF curves are calculated
#' @param cols a vector of colors for the seperate IDF curves, needs same length as \code{probs}
#' @param station.name \code{character} overall naming of the IDF plot, e.g. name of location or model name
#' @param data.name \code{character} naming the data points, e.g. obs or model name
#' @return $ints vector of sorted intensities for selected aggregation levels
#' @return $durs vector of sorted aggregation levels
#' @return $min minimum value of negative log-likelihood during optimization
#' @return $par vector of estimated IDF model parameters mu,sigma,xi,theta,eta at minimum value of negative log-likelihood.
#' @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,3,6,12,24))
#' int.vec <- data.agg$int.vec
#' durs <- data.agg$durs
#' n.y <- data.agg$n.y
#' fit <- IDF.short(int.vec,durs,n.y)
#' pars <- fit$par
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
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,
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") {
###################################################################################
### 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)) {
## 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)
}
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
legend.text.2 <- "quantile"
}
######################################################
### Estimate parameters for duration-dependent GEV ###
######################################################
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
}
######################################################
### success? Than plot! ###
######################################################
if(plot&& !is.na(fit$min)) {
d.all <- unique(durs)
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")
}
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)) {
## 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))))
## 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 function IDF.plot
###################################################################################
} ## 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"
\ No newline at end of file
## 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
###################################################################################
man/IDF.Rd
View file @ beaf7cb8
......@@ -62,7 +62,7 @@ $min minimum value of negative log-likelihood during optimization
$par vector of estimated IDF model parameters mu,sigma,xi,theta,eta at minimum value of negative log-likelihood.
}
\description{
The function \code{IDF.fit} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
The function \code{IDF} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
to a data.frame of observations \code{data} with temporal inforamtion (at least years) and values of precipitation
at a given temporal resoultion. This precipitation time series gets aggregated at given aggregation levels.
\code{agg.lev} and yearly maxima of intensity are caluclated for a specific month or the whole year/dataset.
......
man/IDF.short.Rd 0 → 100644
View file @ beaf7cb8
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/IDF.R
\name{IDF.short}
\alias{IDF.short}
\title{Fitting IDF model parameters to annual maximum intensity time series}
\usage{
IDF.short(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,
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")
}
\arguments{
\item{durs}{a vector of aggregation levels used to fit the IDF curves. One value for each year. Has to have same lenght as \code{int.vec}}
\item{n.y}{\code{integer} value specifying the number of years of data}
\item{mu.init}{initial estimation of location parameter, default is NA. Initial value estimated by fitting individual gev parameters}
\item{sigma.init}{initial estimation of scale parameter,default is NA. Initial value estimated by fitting individual gev parameters}
\item{xi.init}{inital estimation of shape parameter, default is NA. Initial value estimated by fitting individual gev parameters}
\item{theta.init}{inital value defining the curvature of the IDF, default is zero, it is not recommended to change it}
\item{eta.init}{intial estimation of slope parameter for sigma-power law, default is NA. Initial value estimated by fitting individual gev parameters}
\item{use.log}{\code{logical} value for usage of logarithmic values, default is \code{FALSE}}
\item{DEBUG}{\code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative}
\item{method}{\code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e}
\item{upper}{\code{vector} specifying the upper boundary of parameters for "L-BFGS-B" method}
\item{lower}{\code{vector} specifying the lower boundary of parameters for "L-BFGS-B" method}
\item{plot}{\code{logical} option of creating a plot of IDF curves with estimated parameters.}
\item{probs}{a vector of probabilities for which the IDF curves are calculated}
\item{cols}{a vector of colors for the seperate IDF curves, needs same length as \code{probs}}
\item{station.name}{\code{character} overall naming of the IDF plot, e.g. name of location or model name}
\item{data.name}{\code{character} naming the data points, e.g. obs or model name}
\item{ints.vec}{a \code{vector} of yearly maxima of intensity sorted by year and aggregatin level}
}
\value{
$ints vector of sorted intensities for selected aggregation levels
$durs vector of sorted aggregation levels
$min minimum value of negative log-likelihood during optimization
$par vector of estimated IDF model parameters mu,sigma,xi,theta,eta at minimum value of negative log-likelihood.
}
\description{
The function \code{IDF.short} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
to vectors of annnual maximum intensities \code{int.vec} at different durations \code{durs}.
The starting values of the IDF model parameters can be determined by the user as well as specific options to use
during optimization. Logartihmic transformation, debugging, the optimization method, and an option to plot the
IDF curves.
}
\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,3,6,12,24))
int.vec <- data.agg$int.vec
durs <- data.agg$durs
n.y <- data.agg$n.y
fit <- IDF.short(int.vec,durs,n.y)
pars <- fit$par
}
\author{
Christoph Ritschel \email{christoph.ritschel@met.fu-berlin.de}
}
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