# This file contains the functions: # -IDF.agg for the preparing the data # -IDF.plot for plotting of IDF curves at a chosen station #### IDF.agg #### #' Aggregation and annual maxima for choosen durations #' @description Aggregates several time series for chosen durations and finds annual maxima #' (either for the whole year or chosen months). Returns data.frame that can be used for #' the function \code{\link{gev.d.fit}}. #' #' @param data list of data.frames containing time series for every station. #' The data.frame must have the columns 'date' and 'RR' unless other names #' are specified in the parameter `names`. The column 'date' must contain strings with #' standard date format. #' @param ds numeric vector of aggregation durations. #' (Must be multiples of time resolution at all stations.) #' @param na.accept numeric in [0,1] giving maximum percentage of missing values #' for which block max. should still be calculated #' @param which.stations optional, subset of stations. Either numeric vector or character vector #' containing names of elements in data. If not given, all elements in `data` will be used. #' @param which.mon optional, subset of months of which to calculate the annual maxima from. #' @param names optional, character vector of length 2, containing the names of the columns to be used. #' @param cl optional, number of cores to be used from \code{\link[pbapply]{pblapply}} for parallelization. #' #' @details If data contains stations with different time resolutions that need to be aggregated at #' different durations, IDF.agg needs to be run separately for the different groups of stations. #' Afterwards the results can be joint together using `rbind`. #' #' @return data.frame containing the annual intensity maxima [mm/h] in `$xdat`, the corresponding duration in `$ds` #' and the station id or name in `$station`. #' #' @seealso \code{\link{pgev.d}} #' #' @export #' @importFrom pbapply pbsapply #' @importFrom RcppRoll roll_sum #' #' @examples #' dates <- as.Date("2019-01-01")+0:729 #' x <- rgamma(n = 730, shape = 0.4, rate = 0.5) #' df <- data.frame(date=dates,RR=x) #' IDF.agg(list(df),ds=c(24,48)) #' #'## xdat ds station #'## 1 0.3025660 24 1 #'## 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 = list(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('At least one of the given aggregation durations is not multiple of the time resolution = ' ,dtime,' of station ',station,'.')} # function 1: aggregate over single durations and find annual maxima: agg.ts <- function(ds){ runsum = RcppRoll::roll_sum(data.s[,names[2]],ds/dtime,fill=NA) #runmean <- rollapplyr(as.zoo(data.s[,names[2]]),ds/dtime,FUN=sum,fill =NA,align='right') runsum <- runsum/ds #intensity per hour max.subset <- lapply(1:length(which.mon),function(m.i){ subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon[[m.i]]) max <- tapply(runsum[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset], function(vec){ n.na <- sum(is.na(vec)) max <- ifelse(n.na <= na.accept*length(vec),max(vec,na.rm = TRUE),NA) return(max)}) df <- data.frame(xdat=max,ds=ds,year=as.numeric(names(max)),mon=deparse(which.mon[[m.i]]), stringsAsFactors = FALSE) return(df)}) df <- do.call(rbind,max.subset) return(df) # maxima for single durations } # call function 1 in lapply to aggregate over all durations at single station data.agg <- pbapply::pblapply(ds,agg.ts,cl=cl) # df <- do.call(rbind,data.agg) 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 #### #' Plotting of IDF curves at a chosen station #' #' @param durations vector of durations for which to calculate the quantiles. #' @param fitparams vector containing parameters mut, sigma0, xi, theta, eta #' (modified location, scale, shape, duration offset, duration exponent) for chosen station #' as obtained from \code{\link{gev.d.fit}} #' (or \code{\link{gev.d.params}} for model with covariates). #' @param probs vector of exeedance probabilities for which to plot IDF curves (p = 1-1/ReturnPeriod) #' @param cols vector of colors for IDF curves. Should have same length as \code{probs} #' @param add logical indicating if plot should be added to existing plot, default is FALSE #' @param legend logical indicating if legend should be plotted (TRUE, the default) #' @param ... additional parameters passed on to the \code{plot} function #' #' @export #' @importFrom grDevices rgb #' @importFrom graphics axis box lines plot points #' @examples #' data('example',package = 'IDF') #' # fit d-gev #' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")]) #' ,mul = c(1,2),sigl = 1) #' # get parameters for cov1 = 1, cov2 = 1 #' par <- gev.d.params(fit = fit, ydat = matrix(1,1,2)) #' # plot quantiles #' IDF.plot(durations = seq(0.5,35,0.2),fitparams = par) #' # add data points #' points(example[example$cov1==1,]$d,example[example$cov1==1,]$dat) 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(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs)) ## calculate IDF values for given probability and durations qs <- lapply(durations,qgev.d,p=probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3], theta=fitparams[4],eta=fitparams[5]) idf.array <- simplify2array(qs) # array[probs,durs] if(!add){ #new plot ## initialize plot window # check if limits were passed if(is.element('ylim',names(list(...)))){ ylim <- list(...)[['ylim']] }else{ylim <- range(idf.array,na.rm=TRUE)} if(is.element('ylim',names(list(...)))){ xlim <- list(...)[['xlim']] }else{xlim <- range(durations,na.rm=TRUE)} # empty plot plot(NA,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy") } ## plot IDF curves for(i in 1:length(probs)){ lines(durations,idf.array[i,],col=cols[i],...) } if(legend){## plot legend # check if lwd, lty were passed if(is.element('lwd',names(list(...)))){ lwd <- list(...)[['lwd']] }else{lwd <- 1} if(is.element('lty',names(list(...)))){ lty <- list(...)[['lty']] }else{lty <- 1} legend(x="topright",title = 'p-quantile',legend=probs, col=cols,lty=lty,lwd=lwd) } }