README.md 3.43 KB
Newer Older
Jana Ulrich's avatar
Jana Ulrich committed
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

<!-- 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
analysis-tool in hydrology to assessthe characteristics of extreme
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")
Jana Ulrich's avatar
Jana Ulrich committed
35
36
37
38
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
```

## 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
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]')
```

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

  - 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
77
#> [1] 62.01441
Jana Ulrich's avatar
Jana Ulrich committed
78
79
#> 
#> $mle
80
#> [1]  5.983501e+00  4.844650e-01 -1.860657e-02  2.126704e-08  7.908172e-01
Jana Ulrich's avatar
Jana Ulrich committed
81
82
#> 
#> $se
83
#> [1] 3.974995e-01 8.085273e-02 8.330635e-02 2.000063e-06 1.230878e-02
Jana Ulrich's avatar
Jana Ulrich committed
84
85
86
87
88
89
90
91
92
93
# 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)
94
95
#>        mut   sigma0          xi        theta       eta
#> 1 5.983501 1.623306 -0.01860657 2.126704e-08 0.7908172
Jana Ulrich's avatar
Jana Ulrich committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

# 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:
  lines(q,dgev.d(q,params$mut,params$sigma0,params$xi,params$theta,params$eta,d = d),col=d) 
}
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%" />