IDF.R 7.74 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# 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.)
#' @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`.
#'
#' @return data.frame containing the annual maxima in `$xdat`, the corresponding duration in `$ds` 
#' and the station id or name in `$station`.
#' 
#' @seealso \code{\link{pgev.d}}
#' 
#' @export
#' @importFrom pbapply pblapply 
#' @importFrom zoo rollapply 
#'
Christoph Ritschel's avatar
Christoph Ritschel committed
38
#' @examples
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#' 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,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 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){
      runmean <- rollapply(data.s[,names[2]],ds/dtime,FUN=sum,fill =NA,align='right',na.rm= TRUE)
      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],max,na.rm=T)
      return(max) # maxima for single durations
Christoph Ritschel's avatar
Christoph Ritschel committed
73
    }
74
75
    # call function 1 in lapply to aggregate over all durations at single station
    data.agg <- sapply(ds,agg.ts,simplify=TRUE)  
Christoph Ritschel's avatar
Christoph Ritschel committed
76
    
77
78
79
    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
80
  }
81
82
83
84
  # 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 <- pblapply(which.stations,agg.station,cl=cl)
Rust Henning's avatar
Rust Henning committed
85

86
  return(do.call('rbind',station.list))
Christoph Ritschel's avatar
Christoph Ritschel committed
87
}
Rust Henning's avatar
Rust Henning committed
88

89

90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#### 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
#'
112
#' @export
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
#' @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)
#' par <- gev.d.params(fit = fit,cov.list = list(1,1))
#' 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')
    axis(1,at=calc.dur,labels=calc.dur)
    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))),
           pch=c(pch,rep(NA,length(probs))),lwd=c(NA,rep(lwd,length(probs))))
167
168
  }
}