The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
Intensity-duration-frequency (IDF) curves are a widely used analysis-tool in hydrology to assess the 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).
You can install the released version of IDF from CRAN with:
install.packages("IDF")or from gitlab using:
devtools::install_git("https://gitlab.met.fu-berlin.de/Rpackages/idf_package")Here are a few examples to illustrate the order in which the functions are intended to be used.
set.seed(999)
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)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]')
fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log'))
#> $conv
#> [1] 0
#>
#> $nllh
#> [1] 59.34496
#>
#> $mle
#> [1] 6.478887e+00 3.817184e-01 -1.833254e-02 2.843666e-09 7.922310e-01
#>
#> $se
#> [1] 0.25207846 0.02370771 0.04861600 NaN 0.01008561
# checking the fit
gev.d.diag(fit,pch=1,ci=TRUE)
# parameter estimates
params <- gev.d.params(fit)
print(params)
#> mut sigma0 xi theta eta eta2 tau
#> 1 6.478887 1.4648 -0.01833254 2.843666e-09 0.792231 0 0
# 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,params$tau,d = d),col=d)
}
legend('topright',col=1:2,lwd=1,legend = paste('d=',1:2,'h'),title = 'Duration')
plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/h]')
IDF.plot(durations,params,add=TRUE)
This Example depicts the different features that can be used to model the IDF curves, see Fauer et. al (2021, https://doi.org/10.5194/hess-25-6479-2021). Here we assume, that the block maxima of each duration can be modeled with the GEV distribution ():
where the GEV parameters depend on duration according to:
[ \sigma(d)=\sigma_0(d+\theta)^{-\eta_2}+\tau, \\ \mu(d) = \tilde{\mu}\cdot\sigma_0(d+\theta)^{-\eta}+\tau, \\ \xi(d) = \text{const.} ](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 ” (d)=_0(d+)^{-_2}+, \ (d) = _0(d+)^{-}+, \ (d) = “)
The function gev.d.fit provides the options:
theta_zero = TRUE eta2_zero = TRUE (default) tau_zero = TRUE (default) resulting in the following features for IDF-curves:
Example:
### 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.2)
# 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)
# group fits
all.fits <- list(simple=fit.simple,curvature=fit.theta,multiscaling=fit.eta2,flattening=fit.tau)
# compare parameter estimates:
print(t(sapply(all.fits,gev.d.params)))
#> mut sigma0 xi theta eta eta2
#> simple 1.754974 2.897875 -0.004559432 0 0.4978471 0
#> curvature 1.782777 3.187756 -0.01039414 0.02949747 0.5376237 0
#> multiscaling 1.896368 2.826447 0.1358611 0.03034914 0.4971133 0.1577612
#> flattening 2.038391 2.712067 0.1469354 0.09235748 0.6237206 0.2642782
#> tau
#> simple 0
#> curvature 0
#> multiscaling 0
#> flattening 0.1287691
### compare resulting idf-curves
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)
# create 4 plots: one for each additional parameter
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)){
# 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])
}
mtext(fit.labels[i.fit],3,-1.25)
}
legend('topright',lty=rev(1:3),legend = rev(idf.probs),title = 'p-quantile')
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.