README.md 8.17 KB
Newer Older
Jana Ulrich's avatar
Jana Ulrich committed
1 2 3 4 5 6 7 8 9 10

<!-- README.md is generated from README.Rmd. Please edit that file -->

# IDF

<!-- badges: start -->

<!-- badges: end -->

Intensity-duration-frequency (IDF) curves are a widely used
Jana Ulrich's avatar
Jana Ulrich committed
11
analysis-tool in hydrology to assess the characteristics of extreme
Jana Ulrich's avatar
Jana Ulrich committed
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
precipitation. The package ‘IDF’ functions to estimate IDF relations for
given precipitation time series on the basis of a duration-dependent
generalized extreme value (GEV) distribution. The central function is ,
which uses the method of maximum-likelihood estimation for the d-GEV
parameters, whereby it is possible to include generalized linear
modeling for each parameter. For more detailed information on the
methods and the application of the package for estimating IDF curves
with spatial covariates, see Ulrich et. al (2020,
<https://doi.org/10.3390/w12113119>).

## Installation

You can install the released version of IDF from
[CRAN](https://CRAN.R-project.org) with:

``` r
install.packages("IDF")
```

or from gitlab using:

``` r
34
devtools::install_git("https://gitlab.met.fu-berlin.de/Rpackages/idf_package")
Jana Ulrich's avatar
Jana Ulrich committed
35 36 37 38 39 40 41 42 43 44 45 46
```

## Example

Here are a few examples to illustrate the order in which the functions
are intended to be used.

  - Step 0: sample 20 years of example hourly ‘precipitation’ data

<!-- end list -->

``` r
47
set.seed(999)
Jana Ulrich's avatar
Jana Ulrich committed
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
dates <- seq(as.POSIXct("2000-01-01 00:00:00"),as.POSIXct("2019-12-31 23:00:00"),by = 'hour')
sample.precip <- rgamma(n = length(dates), shape = 0.05, rate = 0.4)
precip.df <- data.frame(date=dates,RR=sample.precip)
```

  - Step 1: get annual maxima

<!-- end list -->

``` r
library(IDF)

durations <- 2^(0:6) # accumulation durations [h] 
ann.max <- IDF.agg(list(precip.df),ds=durations,na.accept = 0.1)
# plotting the annual maxima in log-log representation
plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/h]')
```

66
<img src="man/figures/README-annualmax-1.png" width="100%" />
Jana Ulrich's avatar
Jana Ulrich committed
67 68 69 70 71 72 73 74 75 76 77

  - Step 2: fit d-GEV to annual maxima

<!-- end list -->

``` r
fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log'))
#> $conv
#> [1] 0
#> 
#> $nllh
78
#> [1] 59.34496
Jana Ulrich's avatar
Jana Ulrich committed
79 80
#> 
#> $mle
81
#> [1]  6.478887e+00  3.817184e-01 -1.833254e-02  2.843666e-09  7.922310e-01
Jana Ulrich's avatar
Jana Ulrich committed
82 83
#> 
#> $se
84
#> [1] 4.018974e-01 6.677404e-02 4.931317e-02 2.000063e-06 1.141522e-02
Jana Ulrich's avatar
Jana Ulrich committed
85 86 87 88 89 90 91 92 93 94
# checking the fit 
gev.d.diag(fit,pch=1,)
```

<img src="man/figures/README-fit-1.png" width="100%" />

``` r
# parameter estimates 
params <- gev.d.params(fit)
print(params)
95 96
#>        mut sigma0          xi        theta      eta     eta2 tau
#> 1 6.478887 1.4648 -0.01833254 2.843666e-09 0.792231 0.792231   0
Jana Ulrich's avatar
Jana Ulrich committed
97 98 99 100 101 102 103 104 105 106 107

# plotting the probability density for a single duration 
q.min <- floor(min(ann.max$xdat[ann.max$ds%in%1:2]))
q.max <- ceiling(max(ann.max$xdat[ann.max$ds%in%1:2]))
q <- seq(q.min,q.max,0.2)
plot(range(q),c(0,0.55),type = 'n',xlab = 'Intensity [mm/h]',ylab = 'Density')
for(d in 1:2){ # d=1h and d=2h
  # sampled data:
  hist(ann.max$xdat[ann.max$ds==d],main = paste('d=',d),q.min:q.max
       ,freq = FALSE,add=TRUE,border = d)   
  # etimated prob. density:
108
  lines(q,dgev.d(q,params$mut,params$sigma0,params$xi,params$theta,params$eta,params$tau,d = d),col=d) 
Jana Ulrich's avatar
Jana Ulrich committed
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
}
legend('topright',col=1:2,lwd=1,legend = paste('d=',1:2,'h'),title = 'Duration')
```

<img src="man/figures/README-fit-2.png" width="100%" />

  - Step 3: adding the IDF-curves to the data

<!-- end list -->

``` r
plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/h]')
IDF.plot(durations,params,add=TRUE)
```

<img src="man/figures/README-idf-1.png" width="100%" />
125 126 127 128

## IDF Features

This Example depicts the different features that can be used to model
129 130 131 132
the IDF curves (submitted for publication: Fauer et al. An Extended
Model in Estimating Consistent Quantiles for
Intensity-Duration-Frequency Curves, 2020). Here we assume, that the
block maxima of each duration can be modeled with the GEV distribution
Jana Ulrich's avatar
Jana Ulrich committed
133 134 135 136 137 138 139 140
(![\\xi\\neq0](https://latex.codecogs.com/png.latex?%5Cxi%5Cneq0
"\\xi\\neq0")):   
![G(z;\\mu,\\sigma,\\xi)=\\exp \\left\\lbrace -\\left\[ &#10;1+\\xi
\\left( \\frac{z-\\mu}{\\sigma} \\right)&#10;\\right\]^{-1/\\xi}
\\right\\rbrace,](https://latex.codecogs.com/png.latex?G%28z%3B%5Cmu%2C%5Csigma%2C%5Cxi%29%3D%5Cexp%20%5Cleft%5Clbrace%20-%5Cleft%5B%20%0A1%2B%5Cxi%20%5Cleft%28%20%5Cfrac%7Bz-%5Cmu%7D%7B%5Csigma%7D%20%5Cright%29%0A%5Cright%5D%5E%7B-1%2F%5Cxi%7D%20%5Cright%5Crbrace%2C
"G(z;\\mu,\\sigma,\\xi)=\\exp \\left\\lbrace -\\left[ 
1+\\xi \\left( \\frac{z-\\mu}{\\sigma} \\right)
\\right]^{-1/\\xi} \\right\\rbrace,")  
141 142 143

where the GEV parameters depend on duration according to:

Jana Ulrich's avatar
Jana Ulrich committed
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
  
![&#10;\\sigma(d)=\\sigma\_0(d+\\theta)^{-\\eta\_2}+\\tau,
\\\\&#10;\\mu(d) =
\\tilde{\\mu}\\cdot\\sigma\_0(d+\\theta)^{-\\eta}+\\tau,
\\\\&#10;\\xi(d) = \\text{const.}
&#10;](https://latex.codecogs.com/png.latex?%0A%5Csigma%28d%29%3D%5Csigma_0%28d%2B%5Ctheta%29%5E%7B-%5Ceta_2%7D%2B%5Ctau%2C%20%5C%5C%0A%5Cmu%28d%29%20%3D%20%5Ctilde%7B%5Cmu%7D%5Ccdot%5Csigma_0%28d%2B%5Ctheta%29%5E%7B-%5Ceta%7D%2B%5Ctau%2C%20%20%5C%5C%0A%5Cxi%28d%29%20%3D%20%5Ctext%7Bconst.%7D%20%0A
"
\\sigma(d)=\\sigma_0(d+\\theta)^{-\\eta_2}+\\tau, \\\\
\\mu(d) = \\tilde{\\mu}\\cdot\\sigma_0(d+\\theta)^{-\\eta}+\\tau,  \\\\
\\xi(d) = \\text{const.} 
")  

The function `gev.d.fit` provides the options:

  - `theta_zero = TRUE` ![\\theta
    = 0](https://latex.codecogs.com/png.latex?%5Ctheta%20%3D%200
    "\\theta = 0")
161
  - `eta2_zero = TRUE` (default) ![\\eta\_2 =
Jana Ulrich's avatar
Jana Ulrich committed
162 163
    \\eta](https://latex.codecogs.com/png.latex?%5Ceta_2%20%3D%20%5Ceta
    "\\eta_2 = \\eta")
164
  - `tau_zero = TRUE` (default) ![\\tau
Jana Ulrich's avatar
Jana Ulrich committed
165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
    = 0](https://latex.codecogs.com/png.latex?%5Ctau%20%3D%200
    "\\tau = 0")

resulting in the following features for IDF-curves:

  - simple scaling: using only parameters ![\\tilde{\\mu}, \\sigma\_0,
    \\xi,
    \\eta](https://latex.codecogs.com/png.latex?%5Ctilde%7B%5Cmu%7D%2C%20%5Csigma_0%2C%20%5Cxi%2C%20%5Ceta
    "\\tilde{\\mu}, \\sigma_0, \\xi, \\eta")
  - curvature for small durations: allowing ![\\theta
    \\neq 0](https://latex.codecogs.com/png.latex?%5Ctheta%20%5Cneq%200
    "\\theta \\neq 0") (default)
  - multi-scaling: allowing ![\\eta\_2 \\neq
    \\eta](https://latex.codecogs.com/png.latex?%5Ceta_2%20%5Cneq%20%5Ceta
    "\\eta_2 \\neq \\eta")
  - flattening for long durations: allowing ![\\tau
    \\neq 0](https://latex.codecogs.com/png.latex?%5Ctau%20%5Cneq%200
    "\\tau \\neq 0").
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208

Example:

``` r
### sampling example data
set.seed(42)
# durations
ds <- 1/60*2^(seq(0,13,1))
# random data for each duration
xdat <- sapply(ds,rgev.d,n = 20,mut = 2,sigma0 =3,xi = 0.2,theta = 0.1,eta = 0.6,tau = 0.1,eta2 = 0.8)
# transform to data.frame
example <- data.frame(xdat=as.numeric(xdat),ds=rep(ds,each=dim(xdat)[1]))

### different fit options
fit.simple <- gev.d.fit(xdat=example$xdat,ds = example$ds,theta_zero = TRUE,show=FALSE)
fit.theta <- gev.d.fit(xdat=example$xdat,ds = example$ds,show=FALSE)
fit.eta2 <- gev.d.fit(xdat=example$xdat,ds = example$ds,eta2_zero = FALSE,show=FALSE)
fit.tau <- gev.d.fit(xdat=example$xdat,ds = example$ds,eta2_zero = FALSE,tau_zero = FALSE,show=FALSE)

### compare resulting idf-curves
# group fits
all.fits <- list(fit.simple,fit.theta,fit.eta2,fit.tau)
fit.cols <- c('red','purple','blue','darkgreen')
fit.labels <- c('simple-scaling','theta!=0','eta2!=eta','tau!=0')
# plotting probabilities
idf.probs <- c(0.5,0.75,0.99) 
209
# create 4 plots: one for each additional parameter
210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
par(mfrow=c(2,2),mar=c(0.2,0.2,0.2,0.2),oma=c(3.5,4.5,0,0),mgp=c(2.5,0.6,0))
for(i.fit in 1:length(all.fits)){
  plot(example$ds,example$xdat,log='xy',type='n',axes=FALSE)
  box()
  boxplot(example$xdat~example$ds,at=ds,add = TRUE,
        boxwex=0.2,cex=0.4,axes=FALSE)
  if(i.fit %in% c(1,3)){
    axis(2,las=2)
    mtext('Intensity [mm/h]',2,3)
  }
  if(i.fit %in% 3:4){
    axis(1,at=c(0.1,1,24,120),labels = c(0.1,1,24,120))
    mtext('Duration [h]',1,2)
  }
  for(i.p in 1:length(idf.probs)){
Jana Ulrich's avatar
Jana Ulrich committed
225 226 227
    # plotting IDF curves for each model (different colors) and probability (different lty)
    IDF.plot(1/60*2^(seq(0,13,0.5)),gev.d.params(all.fits[[i.fit]]),probs = idf.probs[i.p]
             ,add = TRUE,legend = FALSE,lty = i.p,cols = fit.cols[i.fit])
228 229 230 231 232 233 234
  }
  mtext(fit.labels[i.fit],3,-1.25)
}
legend('topright',lty=rev(1:3),legend = rev(idf.probs),title = 'p-quantile')
```

<img src="man/figures/README-idf_features-1.png" width="100%" />