IDF.R 8.07 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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 onther 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.)
Jana Ulrich's avatar
Jana Ulrich committed
19
#' @param na.accept numeric giving maximum number of missing values for which annual max. should still be calculated
20
21
22
23
24
25
26
27
28
29
#' @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 seperately for the different groups of stations. 
#' Afterwards he results can be joint together using `rbind`.
#'
Jana Ulrich's avatar
Jana Ulrich committed
30
#' @return data.frame containing the annual intensity maxima [mm/h] in `$xdat`, the corresponding duration in `$ds` 
31
32
33
34
35
#' and the station id or name in `$station`.
#' 
#' @seealso \code{\link{pgev.d}}
#' 
#' @export
Jana Ulrich's avatar
Jana Ulrich committed
36
#' @importFrom pbapply pbsapply 
37
38
#' @importFrom zoo rollapply 
#'
Christoph Ritschel's avatar
Christoph Ritschel committed
39
#' @examples
40
41
42
43
44
45
46
47
48
49
#' 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
Jana Ulrich's avatar
Jana Ulrich committed
50
51
52
IDF.agg <- function(data,ds,na.accept = 0,
                    which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){
  
53
  if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data))
Jana Ulrich's avatar
Jana Ulrich committed
54
  
55
56
57
58
  # 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 "
Jana Ulrich's avatar
Jana Ulrich committed
59
                                    ,station," contains: ", class(data.s))}
60
    if(sum(is.element(names[1:2],names(data.s)))!=2){stop('Dataframe of station ', station 
Jana Ulrich's avatar
Jana Ulrich committed
61
62
                                                          ,' does not contain $', names[1]
                                                          ,' or $', names[2], '.')}
63
64
65
66
67
68
69
70
    dtime<-as.numeric((data.s[,names[1]][2]-data.s[,names[1]][1]),units="hours")
    
    if(any(ds %% dtime != 0)){
      stop('Given aggragation 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){
Jana Ulrich's avatar
Jana Ulrich committed
71
      runmean <- rollapply(data.s[,names[2]],ds/dtime,FUN=sum,fill =NA,align='right')
72
73
      runmean <- runmean/ds #intensity per hour
      subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
Jana Ulrich's avatar
Jana Ulrich committed
74
75
76
77
78
      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)})
79
      return(max) # maxima for single durations
Christoph Ritschel's avatar
Christoph Ritschel committed
80
    }
81
    # call function 1 in lapply to aggregate over all durations at single station
Jana Ulrich's avatar
Jana Ulrich committed
82
    data.agg <- pbsapply(ds,agg.ts,simplify=TRUE,cl=cl)  
Christoph Ritschel's avatar
Christoph Ritschel committed
83
    
84
85
86
    df <- data.frame(xdat = as.vector(data.agg), ds = rep(ds,each=length(data.agg[,1])))
    df$station <- rep(station,length(df$xdat))
    return(df) # maxima for all durations at one station
Rust Henning's avatar
Rust Henning committed
87
  }
88
89
90
  # 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
Jana Ulrich's avatar
Jana Ulrich committed
91
92
  station.list <- lapply(which.stations,agg.station)
  
93
  return(do.call('rbind',station.list))
Christoph Ritschel's avatar
Christoph Ritschel committed
94
}
Rust Henning's avatar
Rust Henning committed
95

96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#### IDF.plot ####



#' Plotting of IDF curves at a chosen station
#'
#' @param data matrix or dataframe containing: first column maxima,
#' second column coresponding durations
#' @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}}.
#' @param probs vector of exeedance probabilities for which to plot IDF curves (p = 1-1/ReturnPeriod)
#' @param calc.dur vector of durations for which to calculate IDF curves. If `NULL` (the default),
#' durations from data are taken
#' @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
#' @param xlim,ylim vectors of x- / y-plot-range
#' @param legend logical indicating if legend should be plotted
#' @param st.name string containing legend title (for example station name)
#' @param dt.name string containing name for data points in legend
#' @param ... additional parameters passed on to the \code{plot} function
#'
118
#' @export
119
120
121
122
123
124
#' @importFrom grDevices rgb
#' @importFrom graphics axis box lines plot points 
#' @examples
#' data('example',package = 'IDF')
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
#'                  ,mul = c(1,2),sigl = 1)
Jana Ulrich's avatar
Jana Ulrich committed
125
#' par <- gev.d.params(fit = fit, ydat = cbind(1,1))
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#' IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par),
#'          calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example')
IDF.plot <- function(data,fitparams,probs=c(0.5,0.9,0.99),calc.dur=NULL,
                     cols=4:2,add=FALSE,ylim=NULL,xlim=NULL,
                     legend=TRUE, st.name='Station',dt.name="observations",...){
  
  # if cols is to short, make longer
  if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
  
  # if no calc.dur given, take durations from data
  if(is.null(calc.dur))calc.dur <- data[,2]

  ## calculate IDF values for given probability and durations
  idf.array <- simplify2array(qgev.d(probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3],
                                     theta=fitparams[4],eta=fitparams[5],d=calc.dur)) # array[probs,durs]
  if(!add){ #new plot
    ## initialize plot window
    # check if limits where passed
    if(is.null(ylim)){ylim <- range(idf.array,data[,1],na.rm=TRUE)}
    if(is.null(xlim)){xlim <- range(data[,2],calc.dur,na.rm=TRUE)}

    # empty plot
    plot(NA,axes=F,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy",...)
    box('plot')
Jana Ulrich's avatar
Jana Ulrich committed
150
    axis(1,at=data[,2],labels=round(data[,2],digits = 2))
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
    axis(2)

    # check if `pch` was passed in ...
    if(is.element('pch',names(list(...))))pch <- list(...)[['pch']] else pch <- 16

    ## plot data points
    points(data[,2],data[,1],pch=pch,col=rgb(0,0,0,0.5))

  }else points(data[,2],data[,1],pch=16,col=rgb(0,0,0,0.5)) # add to existing plot

  # check if `lty,ltd` were passed in ...
  if(is.element('lty',names(list(...))))lty <- list(...)[['lty']] else lty <- 1
  if(is.element('lwd',names(list(...))))lwd <- list(...)[['lwd']] else lwd <- 1
    
  ## plot IDF curves
  for(i in 1:length(probs))
    lines(calc.dur,idf.array[i,],col=cols[i],lty=lty,lwd=lwd)

  if(legend){## plot legend
    legend(x="topright",title = st.name,legend=c(dt.name,paste(probs,"quantile",sep=" ")),
           col=c(rgb(0,0,0,0.5),cols),lty=c(NA,rep(lty,length(probs))),
Jana Ulrich's avatar
Jana Ulrich committed
172
           pch=c(pch,rep(NA,length(probs))),lwd=c(NA,rep(lwd,length(probs))))}
173
}