Commit d98f1874 authored by Rust Henning's avatar Rust Henning 🐿
Browse files

initial commit

parents
Package: IDF
Type: Package
Title: Estimating intensity-duration-frequency parameters
Version: 1.01
Date: 2016-01-28
Author: Christoph Ritschel, Sarah Joedicke
Maintainer: Christoph Ritschel <christoph.ritschel@met.fu-berlin.de>
Description: The package 'IDF' reads precipitation data from webwerdis and Berlin Stadtmessnetz.
It aggregates these data and calculates yearly intensity maxima. From those maxima intensity-duration-frequency
parameters are estimated. IDF curves based on these estimated parameters can be plotted.
Depends: stats4, evd
License: GPL (>=2)
exportPattern(IDF.read,IDF.fit,IDF.plot)
# Import all packages listed as Imports or Depends
import(
stats4,
evd
)
##################################################
## IDF package
## Authors: Sarah Joedicke, Christoph Ritschel
## Update: 29.01.2016
###################################################
###############################################
############# Read Data function ##############
###############################################
##' @title Reading precipitation data
##' @description The function \code{IDF.read} reads a file in table format and creates a \code{list} from it, writing the data frame as first entry
##' and adds some attributes (station information, aggregation time, data source). The only data values used are:
##' date, precipitation
##' The \code{data.frame} will have the following format:
##' | year | mon | day | hour | min | RR |
##' |------+-----+-----+------+-----+----+
##' | | | | | | |
##' @usage IDF.read(file, type)
##' @param file a \code{character string} naming the file from which the data is to be read.
##' @param type a \code{character string} defining the type of data to be read: either "stadtmessnetz" or "webwerdis", depending on if the data comes from the Stadtmessnetz Berlin
##' or WebWerdis. If type = "webwerdis", the data will be read, then sorted, formatted and missing lines added,
##' while if type = "stadtmessnetz", the data will just be read and formatted.
##' Both source types have a different layout in the original file.
##' @return Liste a \code{list} containing a \code{data.frame} of date and time information and precipitation values for each time step
##' @details This function is designed to prepare a data file for aggregation by functions \code{IDF.agg} or \code{IDF.magg}.
##' The time given in the data is the end time, so the precipitation was measured up to that time.
##' @note The first entry of the list will be a \code{data.frame} named 't1'.
##' @seealso read.table, IDF.agg, IDF.magg
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
IDF.read <- function(file,type){
if(type != "stadtmessnetz" && type != "webwerdis") {
cat("Warning: wrong type declared for input file")
stop()
}
if (type == "stadtmessnetz") {
Tab_MN <- read.csv2(file) #STADTMESSNETZ
new_time <- strptime(Tab_MN$Zeitstempel,format="%d.%m.%Y %H:%M") #STADTMESSNETZ Datumsvektor
}
# Da die Stadtmessnetzdaten (bisher) konstistent aussehen, wird auf das Erstellen einer neuen Tabelle mit sicher allen
# Zeiten verzichtet, da die Minutendaten sehr gross sind. Sollte es inkonsistente Tabellen geben, sollte man diese seperat behandeln,
# sonst wird viel Rechenzeit fuer die kompletten Tabellen verschwendet.
if (type == "webwerdis") {
Tab <- read.table(file,header=TRUE,sep=";") #WEBWERDIS
Tab_kurz <- Tab[,c("Date","precipitation")]
## Tabelle in gewuenschtes Format sortieren; nur WEBWERDIS, da Stadtmessnetz sortiert und keine fehlenden Zeilen
time <- strptime(Tab_kurz$Date,format="%Y-%m-%d T %H:%M:%S")
Tab_sort <- Tab_kurz[order(as.character(time)),]
time_sort <- strptime(Tab_sort$Date,format="%Y-%m-%d T %H:%M:%S")
Tab_sort$Date <- as.character(time_sort)
## Fehlende Zeilen in Daten hinzufuegen, Niederschlagswert = NA; nur WEBWERDIS
h_diff <- as.numeric(difftime(format(time_sort[length(time_sort)],"%Y-%m-%d T %H:%M:%S") ,
format(time_sort[1],"%Y-%m-%d T %H:%M:%S"),units="hours")) #h_diff ist die gesamte Zeitdifferenz
new_time <- seq(time_sort[1], length = h_diff+1, by = "hour")
new_tab <- data.frame(Date=as.character(new_time), precipitation=NA) # Tabelle mit NAs und richtigen Zeiten vordefinieren
Tab_na <- (merge(Tab_sort, new_tab, "Date", all.y=TRUE))[,1:2]
} #Daten sortieren und in gewuenschte Form bringen
new_timect <- as.POSIXct(new_time)
J <- as.numeric(format(new_timect,'%Y'))
M <- as.numeric(format(new_timect,'%m'))
d <- as.numeric(format(new_timect,'%d'))
h <- as.numeric(format(new_timect,'%H'))
m <- as.numeric(format(new_timect,'%M'))
if (type == "webwerdis") Tab_end <- data.frame(J,M,d,h,m,Tab_na$precipitation.x) #WEBWERDIS
if (type == "stadtmessnetz") Tab_end <- data.frame(J,M,d,h,m,Tab_MN[,2]) #STADTMESSNETZ
## Tabellen-Attribute benennen:
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)
if (type == "webwerdis"){
# WEBWERDIS:
attr(Liste,"StationName") <- as.character(Tab$Stationname[1])
attr(Liste,"StationID") <- "NA"
attr(Liste,"Long (deg N)") <- Tab$Longitude[1]
attr(Liste,"Lat (deg E)") <- Tab$Latitude[1]
attr(Liste,"Heigth (m)") <- Tab$StationHeight[1]
attr(Liste,"Source") <- "Web-WERDIS"
} #Listen-Attribute benennen
if (type == "stadtmessnetz"){
# STADTMESSNETZ:
attr(Liste,"StationName") <- colnames(Tab_MN)[2]
attr(Liste,"StationID") <- "NA"
attr(Liste,"Long (deg N)") <- "NA"
attr(Liste,"Lat (deg E)") <- "NA"
attr(Liste,"Height (m)") <- "NA"
attr(Liste,"Source") <- "Stadtmessnetz"
} #Listen-Attribute benennen
cat(paste("read.data of", file , "done \n"))
str(Liste) # optional; so sieht man beim Einlesen, womit man es zu tun hat und ob alles geklappt hat
return(Liste)
}
############################################
########## Moving aggragation #############
############################################
##'@title Moving aggregation of a time series
##'@description The function \code{IDF.magg} aggregates the precipitation values of a \code{data.frame} in a \code{list} by moving an aggregation
##'window over the time steps. The list should be created by funcntion \code{IDF.read} in a previous step.
##'@usage IDF.magg(x, ks=c(2,4,8,16,32,64,128))
##'@param x a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param ks a numeric \code{vector} of aggregatation levels. Each value represents an aggregation level and will create
##'a new \code{data.frame} in the output list.
##'@details This function is designed to aggregate data which has been read by function \code{IDF.read}.
##'The time given in the data is the end time, so the precipitation was added up to that time.
##'Every line of data will be used in this process. In the end, the data frames containing aggregated precipitation
##'will be exactly as long as the original one. In the end, a new \code{list} containing several \code{data.frames} will be created.
##'The first one will be the initial \code{data.frame}, now named agg0. The other data frames will be named 'agg' + the aggregation level number from ks.
##'@return Liste a \code{list} containing \code{data.frames} of original and aggregated precipitation time series with time information and
##'precipitation values.
##'@seealso IDF.read, IDF.agg.
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
IDF.magg <- function(x,ks=c(2,4,8,16,32,64,128)){
NS <- x$t1$Niederschlagssumme #Niederschlagsvektor
Liste <- vector("list", length(ks)+1) # Leere Liste vordefinieren, +1 fuer Ursprungstabelle
Liste[[1]] <- x$t1
j <- 2
for (y in ks){
agg <- vector(length=nrow(x$t1))
print(paste("Aggregation ueber Aggregationsstufe", y, "laeuft"))
for (i in y:nrow(x$t1)){
agg[i] <- sum(NS[(i-(y-1)):i])
}
agg[1:(y-1)] <- NA
Tab_agg <- x$t1[,1:5]
Tab_agg[,6] <- agg
colnames(Tab_agg) <- c("Jahr","Monat","Tag","Stunde","Minute","Niederschlagssumme (mm)")
attr(Tab_agg,"Aggregationszeit (min)") <- attr(x$t1, "Aggregationszeit (min)")*y
Liste[[j]] <- Tab_agg
j <- j+1
}
attributes(Liste) <- attributes(x)
names(Liste) <- paste("agg", c(0,ks), sep="")
print("Vorgang abgeschlossen")
return(Liste)
}
############################################
########## diskret aggregation #############
############################################
##'@title Aggregation of a time series for seperated time steps
##'@description The function \code{IDF.agg} aggregates the precipitation values of a \code{data.frame} in a \code{list} in groups.
##'The list should be created by function \code{IDF.read} in a previous step.
##'@usage IDF.agg(x, ks=c(2,4,8,16,32,64,128))
##'@param x a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param ks a numeric \code{vector} of aggregatation levels. Each value represents an aggregation level and will create
##'a new \code{data.frame} in the output list.
##'@return Liste a \code{list} containing \code{data.frames} of original and aggregated precipitation time series with time information and
##'precipitation values.
##'@details This function is designed to aggregate data which has been read by function \code{IDF.read}.
##'The time given in the data is the end time, so the precipitation was added up to that time.
##'In this process the data will be summed up in groups, depending on the chosen aggregation levels. In the end,
##'the aggregated \code{data.frames} will be much shorter than the original one.
##'In the end, a new \code{list} containing several \code{data.frames} will be created. The first one will be the
##'initial data frame, now named agg0. The other data frames will be named 'agg' + the aggregation level number from ks.
##'@seealso read.data, agg.f().
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
IDF.agg <- function(x,ks=c(2,4,8,16,32,64,128)){
NS <- x$t1$RR
date <- paste(as.character(x$t1$year), as.character(x$t1$mon), as.character(x$t1$day),
as.character(x$t1$hour), as.character(x$t1$min))
Liste <- vector("list", length(ks)+1) # Leere Liste vordefinieren, +1 fuer Ursprungstabelle
Liste[[1]] <- x$t1
j <- 2
for (y in ks){
agg_NS <- apply(matrix(data=NS[1:((length(NS)%/%y)*y)], nrow=(length(NS)%/%y), ncol=y, byrow=TRUE), 1, sum)
# %/% dividiert und rundet ganzzahlig ab
# Variable if NS%%y > 0 sind das die unbenutzen Resttage, die letzten "" Zeitschritte werden ignoriert
if (length(NS)%%y > 0) cat("Note: the last", length(NS)%%y,"timesteps are ignored \n")
agg_date <- matrix(data=date[1:((length(NS)%/%y)*y)], nrow=(length(NS)%/%y), ncol=y, byrow=TRUE)[,y]
new_time2 <- (strptime(agg_date, format="%Y %m %d %H %M"))
J <- as.numeric(format(new_time2,'%Y'))
M <- as.numeric(format(new_time2,'%m'))
d <- as.numeric(format(new_time2,'%d'))
h <- as.numeric(format(new_time2,'%H'))
m <- as.numeric(format(new_time2,'%M'))
Tab_agg <- data.frame(J,M,d,h,m,agg_NS)
colnames(Tab_agg) <- c("year","mon","day","hour","min","RR")
attr(Tab_agg,"accumulation time (min)") <- as.numeric(difftime(new_time2[2],new_time2[1]),units="mins")
Liste[[j]] <- Tab_agg
j <- j+1
}
attributes(Liste) <- attributes(x)
names(Liste) <- paste("agg", c(0,ks), sep="")
return(Liste)
} # End of function agg
################################################################################################################
##'@title Yearly precipitation maximum
##'@description The function \code{IDF.max} generates an \code{vector} of yearly precipitation maxima from a \code{data.frame}
##'containing time information and precipitation data at a single aggregation level generated by functions \code{IDF.agg} or \code{IDF.magg}. An option
##'for selection a single month or search the whole year exists.
##'@param df a \code{data.frame} containing date and time information and the corresponding precipitation values at a single aggregation level.
##'@param month a \code{value} defining the month to use. A single number between 1 and 12 selects the month. The \code{character string} "all"
##'can be used if the maximum for all months is needed.
##'@return RR.max a \code{vector} containing the precipitation maxima for each year of a defined month or all months.
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.max <- function(df,month=1) {
years <- unique(df$year)
y.len <- length(years)
RR.max <- array(NA,dim=c(y.len))
for(y in 1:y.len) {
if(month=="all") { index <- which(df$year==years[y])
}else index <- which(df$year==years[y] & df$mon==month)
if(sum(index)>0) RR.max[y] <- max(df$RR[index],na.rm=T)
}
RR.max <- RR.max[!is.na(RR.max)]
return(RR.max)
}
#######################
## Fitting Functions ##
#######################
##'@title Density function of modified generalized extreme value distribution
##'@description The function \code{dgev.d} is a modified version of the function \code{dgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param q Vector of quantiles
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@param log \code{logical} option to use logarithmic parameter values, default=FALSE
##'@seealso \code{\link[evd]{dgev}}
##'@return dgev.d gives the density function
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
dgev.d <- function(q,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,log=FALSE) {
sigma.d <- sigma/(d+theta)^eta
dgev(q,loc=mu*sigma.d,scale=sigma.d,shape=xi,log=log)
}
##'@title Quantile function of modified generalized extreme value distribution
##'@description The function \code{qgev.d} is a modified version of the function \code{qgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param p Vector of probabilities
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@param lower.tail \code{logical} if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
##'@seealso \code{\link[evd]{qgev}}
##'@return qgev.d gives the quantile function
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
qgev.d <- function(p,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,lower.tail=TRUE) {
sigma.d <- sigma/(d+theta)^eta
qgev(p,loc=mu*sigma.d,scale=sigma.d,shape=xi,lower.tail=lower.tail)
}
##'@title Random generation for the modified generalized extreme value distribution
##'@description The function \code{rgev.d} is a modified version of the function \code{rgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param n Number of observations
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@seealso \code{\link[evd]{rgev}}
##'@return rgev.d generates random derivates
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
rgev.d <- function(n,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1) {
## gumbel
sigma.d <- sigma/(d+theta)^eta
rgev(n, loc=mu*sigma.d,scale=sigma.d,shape=xi)
}
#######################################################################
##' @title Negativ log-likelihood of modified GEV
##' @description The function \code{IDF.nll} calculates the negative log-likelihood for a given set of model parameters
##' \code{mu,sigma,xi,theta,eta}, given observations \code{x} and given durations \code{d}. Options for the usage of
##' logartihmic values \code{use.log} and a debugging function \code{DEBUG} are available.
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param x vector of observations at different durations d
##'@param d vector of durations
##'@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.
##'@return retruns weightes negative log-likelihood by number of observatons uesd
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
## mu is the mu~ from Koutsoyiannis
if(use.log){
## ensure that critical parameters are positive
sigma <- exp(sigma)
theta <- exp(theta)
eta <- exp(eta)
}
sigma.d <- sigma/((d+theta)^eta)
if(DEBUG) debug.values <- c(mu,sigma,xi,theta,eta)
## Weibull und Frechet
if(xi!=0){
C <- 1 + xi * (x/sigma.d - mu )
nll <- switch((sum(C<0)>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),
Inf)
# + 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))-sum((Y))-sum(exp(-Y)))
}
if(DEBUG){
cat(debug.values,nll,"\n")
options(digits.secs=6)
## debug.values <- c(debug.values,nll,as.character(Sys.time()))
## write(debug.values,file="optim.log",append=TRUE,ncolumns=length(debug.values))
## cat(debug.values,nll,sum(A<0),"\n")
}
return(nll/length(x))
} # end of function IDF.nll
######################################################################################################
##' @title Fitting function to optimize IDF model parameters
##' @description The function \code{fit.fun} fits IDF model parameters \code{mu,sigma,xi,theta,eta} to a set of given observations \code{obs},
##' typically a series of yearly maxima at different durations \code{d}. Options for using logarithmic parameter values and debugging
##' are given. Also the \code{optim} parameters \code{method} and \code{upper,lower} can be defined.
##' @param obs vector of yearly intensity maxima at different durations. Order: Y1D1, Y2D1,...,YnD1,Y1D2,...YnD2,Y1D3,...,YnDk
##' @param dur vector of durations with same length as \code{obs}. Order: n x D1, n x D2, ... n x Dk
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@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
##' @return $min value of negative log-likelihood at optimization minimum
##' @return $par vector of IDF parameters at optimization minimum
##' @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) {
sigma <- log(sigma)
theta <- log(theta)
eta <- log(eta)
if(method=="L-BFG-S") {
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])
}
}
## check initial value of negative log-Likelihood function
nll <- IDF.nll(mu,sigma,xi,theta,eta,x=obs,d=dur,use.log=use.log,DEBUG=DEBUG)
## if initial value is acceptable...
if(!is.infinite(nll)&!is.na(nll)) {
if(method=="L-BFGS-B") {
## 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
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),
method=method,upper=upper,lower=lower), error=function(e) NULL)#,
#upper=upper,lower=lower)
}else{
## 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
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),
method=method), error=function(e) NULL)#,
#upper=upper,lower=lower)
}
## if there was no error
if(!is.null(fit)) {
fit.min <- fit@min
fit.par <- fit@coef
}else { ## else return NA
fit.min <- NA
fit.par <- rep(NA,5)
} ## end if error
}else { ## else retunr NA
fit.min <- NA
fit.par <- rep(NA,5)
} ## end if initial value..
if(use.log){
fit.par[2] <- exp(fit.par[2])
fit.par[4] <- exp(fit.par[4])
fit.par[5] <- exp(fit.par[5])
}
names(fit.par) <- c("mu","sigma","xi","theta","eta")
return(list("min"=fit.min,"par"=fit.par))
} ## end of function fit.fun
##################################################################################
##' @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}
##' to a list object of observations \code{data} in a data.frame with temporal inforamtion 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.
##' 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 data a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param agg.lev a vector of aggregation levels used to fit the IDF curves.
##'@param month a value specifying the month to be used for estimating the IDF parameters. Type "all" for all months.
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@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 station.name \code{character} naming of data input, e.g. a observational station 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.
##' @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",mu=1,sigma=1,xi=0,theta=0,eta=1,
use.log=FALSE,DEBUG=FALSE,method="Nelder-Mead",upper=Inf,lower=-Inf,plot=FALSE,station.name="Berlin") {
data.agg <- IDF.agg(data,ks=agg.lev)
max.agg <- do.call("cbind",lapply(data.agg,FUN=IDF.max,month=month))
d.all <- c(1,agg.lev)
int <- t(max.agg)/d.all
y.len <- dim(int)[2]
int.vec <- as.vector(t(int))
durs <- rep(d.all,each=y.len)
fit <- fit.fun(obs=int.vec,dur=durs,mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta,use.log=use.log,
DEBUG=DEBUG,method=method,upper=upper,lower=lower)
if(plot&& !is.na(fit$min)) {
ds <- sort(rep(d.all,length(int.vec)/length(d.all)))
IDF.plot(pars=fit$par,name=station.name,ints=int.vec,ds=ds)
}
if(plot && is.na(fit$min)) {
cat("Warning: optimization did not converge and no parameters where 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 name a character string defining the station name for which the IDF curves are calculated, will be printed in legend
##' @param ints \code{vector} of observational intensities (surted by durations)
##' @param ds \code{vector} of durations (same length as intensities)
##' @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)),name="Berlin-Dahlem",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),max(dur)),ylim=c(min(idf.array[,1]),max(idf.array[,3])),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)
}
## plot legend
legend(x="topright",legend=c(name,"obs","IDF 0.5 quantile","IDF 0.9 quantile","IDF 0.99 quantile"),
col=c(1,rgb(0,0,0,0.5),cols[1],cols[2],cols[3]),lty=c(NA,NA,1,1,1),pch=c(NA,16,NA,NA,NA))
} ## end of function IDF.plot
###################################################################################