Commit d5f936e4 authored by Laura Mack's avatar Laura Mack
Browse files

computational accuracy in if-statement

parent 07c6ca9a
......@@ -48,52 +48,52 @@
#'## 2 0.4112304 24 1
#'## 3 0.1650978 48 1
#'## 4 0.2356849 48 1
IDF.agg <- function(data,ds,na.accept = 0,
which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){
if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data))
# function 2: aggregate station data over durations and find annual maxima:
agg.station <- function(station){
data.s <- data[[station]]
if(!is.data.frame(data.s)){stop("Elements of 'data' must be data.frames. But element "
,station," contains: ", class(data.s))}
if(sum(is.element(names[1:2],names(data.s)))!=2){stop('Dataframe of station ', station
,' does not contain $', names[1]
,' or $', names[2], '.')}
dtime<-as.numeric((data.s[,names[1]][2]-data.s[,names[1]][1]),units="hours")
if(any(ds %% dtime != 0)){
stop('Given aggregation durations are not multiples of the time resolution = '
,dtime,' of station ',station,'.')}
# function 1: aggregate over single durations and find annual maxima:
agg.ts <- function(ds){
runmean <- rollapplyr(as.zoo(data.s[,names[2]]),ds/dtime,FUN=sum,fill =NA,align='right')
runmean <- runmean/ds #intensity per hour
subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
max <- tapply(runmean[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],
function(vec){
n.na <- sum(is.na(vec))
max <- ifelse(n.na <= na.accept,max(vec,na.rm = TRUE),NA)
return(max)})
return(max) # maxima for single durations
IDF.agg <- function(data,ds,na.accept = 0,
which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){
if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data))
# function 2: aggregate station data over durations and find annual maxima:
agg.station <- function(station){
data.s <- data[[station]]
if(!is.data.frame(data.s)){stop("Elements of 'data' must be data.frames. But element "
,station," contains: ", class(data.s))}
if(sum(is.element(names[1:2],names(data.s)))!=2){stop('Dataframe of station ', station
,' does not contain $', names[1]
,' or $', names[2], '.')}
dtime<-as.numeric((data.s[,names[1]][2]-data.s[,names[1]][1]),units="hours")
if(any(ds %% dtime > 10e-16)){
stop('Given aggregation durations are not multiples of the time resolution = '
,dtime,' of station ',station,'.')}
# function 1: aggregate over single durations and find annual maxima:
agg.ts <- function(ds){
runmean <- rollapplyr(as.zoo(data.s[,names[2]]),ds/dtime,FUN=sum,fill =NA,align='right')
runmean <- runmean/ds #intensity per hour
subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
max <- tapply(runmean[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],
function(vec){
n.na <- sum(is.na(vec))
max <- ifelse(n.na <= na.accept,max(vec,na.rm = TRUE),NA)
return(max)})
return(max) # maxima for single durations
}
# call function 1 in lapply to aggregate over all durations at single station
data.agg <- pbsapply(ds,agg.ts,simplify=TRUE,cl=cl)
df <- data.frame(xdat = as.vector(data.agg), ds = rep(ds,each=length(data.agg[,1])))
df$station <- station
df$year <- rep(unique(as.POSIXlt(data.s[,names[1]])$year+1900),length(ds))
return(df) # maxima for all durations at one station
}
# which stations should be used?
if(is.null(which.stations))which.stations <- 1:length(data)
# call function 2 in lapply to aggregate over all durations at all stations
station.list <- lapply(which.stations,agg.station)
return(do.call('rbind',station.list))
}
# call function 1 in lapply to aggregate over all durations at single station
data.agg <- pbsapply(ds,agg.ts,simplify=TRUE,cl=cl)
df <- data.frame(xdat = as.vector(data.agg), ds = rep(ds,each=length(data.agg[,1])))
df$station <- station
df$year <- rep(unique(as.POSIXlt(data.s[,names[1]])$year+1900),length(ds))
return(df) # maxima for all durations at one station
}
# which stations should be used?
if(is.null(which.stations))which.stations <- 1:length(data)
# call function 2 in lapply to aggregate over all durations at all stations
station.list <- lapply(which.stations,agg.station)
return(do.call('rbind',station.list))
}
#### IDF.plot ####
......@@ -130,7 +130,7 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
cols=4:2,add=FALSE,
legend=TRUE,...){
# if cols is to short, make longer
# if cols is to short, make longer
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
## calculate IDF values for given probability and durations
......
......@@ -230,18 +230,18 @@ gev.d.fit<-
if(z$trans) { # for nonstationary fit
print(z[c(2, 4)]) # print model, link (3) , conv
# print names of link functions:
cat('$link$name\n')
cat(paste('mu:\t',z$link$mulink$name))
cat(paste('\nsigma:\t',z$link$siglink$name))
cat(paste('\nxi:\t',z$link$shlink$name))
cat(paste('\ntheta:\t',z$link$thetalink$name))
cat(paste('\neta:\t',z$link$etalink$name,'\n\n'))
#cat('$link$name\n')
#cat(paste('mu:\t',z$link$mulink$name))
#cat(paste('\nsigma:\t',z$link$siglink$name))
#cat(paste('\nxi:\t',z$link$shlink$name))
#cat(paste('\ntheta:\t',z$link$thetalink$name))
#cat(paste('\neta:\t',z$link$etalink$name,'\n\n'))
}else{print(z[4])} # for stationary fit print only conv
if(!z$conv){ # if fit converged
print(z[c(5, 7, 9)]) # print nll, mle, se
}
}
class( z) <- "gev.d.fit"
class(z) <- "gev.d.fit"
invisible(z)
}
......
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