Commit 096be90e authored by Henning Rust's avatar Henning Rust
Browse files

HR suggests some modifications

parent d4ef04ee
......@@ -394,6 +394,49 @@ 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
##################################################################################
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) {
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))
}
#' @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 data.frame of observations \code{data} with temporal inforamtion (at least years) and values of precipitation
......@@ -439,45 +482,11 @@ IDF.fit <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=
#########################################################################
### Calculate extreme values for each year and each aggregation level ###
#########################################################################
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) {
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
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)
###############################################
### Estimate Parameters for single duration ###
###############################################
......@@ -533,6 +542,10 @@ IDF.fit <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=
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)
......
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