tyears: Calculate Low-Flow Quantiles for given Return Periods
In lfstat: Calculation of Low Flow Statistics for Daily Stream Flow Data

Description Usage Arguments Details Value Author(s) References See Also Examples

Fits an extreme value distribution using L-moments to the minima of a time series of discharges and subsequently estimates quantiles (the so called T-years event) for given return periods. In the presence of zero flow observations a mixed distribution is fitted.

tyears(lfobj, event = 1/probs, probs = 0.01,
       dist = "wei", check = TRUE, zeta = zetawei, zetawei = NULL,
       plot = TRUE, col = 1, log = TRUE, legend = TRUE,
       rp.axis = "top", rp.lab = "Return period",
       freq.axis = TRUE,
       freq.lab = expression(paste("Frequency " *(italic(F)),
                                   " = Non-Exceedance Probability P ",
                                   (italic(X) <= italic(x)))),
       xlab = expression("Reduced variate,  " * -log(-log(italic(F)))),
       ylab = "Quantile",
       hyearstart = hyear_start(lfobj),
       n = NULL)

`lfobj`	An object of class lfobj or an object which can be coerced to class xts. Either with a single column or with a column named 'discharge'.
`event`	numeric vector specifying the return periods. E.g. `event = 100` will yield the 100 years extreme low flow event.
`probs`	Alternate way to specify the return period of the event.
`dist`	A character vector of distributions to fit. Basically all distributions provided by Hosking's `lmom-package` and their reversed counterparts can be chosen.
`check`	logical, should `check_distribution` get called?
`zeta`	numeric vector of length one for manually setting a lower bound. Only a few distributions allow for a lower bound, namely `'gpa'`, `'ln3'`, `'wak'` and `'wei'`. The default value of `NULL` results in not bounding the distribution, therefore the parameter `zeta` is estimated.
`zetawei`	same as `zeta`
`plot`	logical. If `TRUE`, sample observations as well as estimated quantile functions are plotted.
`col`	numeric or character vector of length one or as long as `dist`, specifying the color used for plotting.
`log`	logical. If `TRUE` probabilities will be plotted on a double logarithmic scale.
`legend`	logical, should a legend be added to the plot?
`rp.axis`	vector of length one, specifying if and how an additional scale bar for the return periods is drawn. Possible choices are `'bottom'`, `'top'` and `'none'`. Alternatively, the position of the scale bar can be specified as an real number between 0 and 1, indicating the y-position of the legend.
`rp.lab`	character vector, text above the scale bar for return periods
`freq.axis`	logical, should an additional abscissa showing the probabilities be drawn on top of the plot?
`freq.lab`	character vector, text above the probability axis
`xlab`	character vector, a label for the x axis
`ylab`	character vector, a label for the y axis
`hyearstart`	vector of length one, providing the start of the hydrological year. This is evaluated by `water_year`. The default is, to retrieve the values stored in the attributes of the `lfobj`.
`n`	Argument 'n' is deprecated and ignored. To apply a moving average, do it prior to calling 'tyears'. See section Examples.

This function is vectorised over dist and event.

According to paragraph 7.4.2 of the WMO manual, special care has to be taken in the presence of zero flow observations. A cdf called G(x) is fitted to the non-zero values of the original time series

If a distribution is fitted which allows for finite lower bound (zeta), and zeta is estimated being negative, estimation is repeated constraining zeta = 0. If this behavior is not desired, the parameter zeta has to be set explicitly.

An object of class evfit, see evfit.

Daniel Koffler and Gregor Laaha

Gustard, A. & Demuth, S. (2009) (Eds) Manual on Low-flow Estimation and Prediction. Operational Hydrology Report No. 50, WMO-No. 1029, 136p.

There are methods for printing summarizing objects of class evfit.

evfit

data("ngaruroro")
ng <- subset(ngaruroro, hyear %in% 1964:2000)

# vector of return periods
rp <- c(1.5, 5, 10, 100)

# Fitting some distributions for the low flows (annual minima)
# and estimating the quantile for arbitrary return periods
y <- tyears(ng, dist = c("gum", "wei", "ln3", "pe3"), event = rp,
            plot = FALSE)

# print()ing the object shows just the return periods
y

# but y is actually a list
str(y)

# there is a summary method, returning L-moments and estimated parameters
summary(y)

plot(y)


# fitting just one distribution, with annotated quantiles
z <- tyears(ng, dist = c("gevR"), event = rp)
rpline(y, return.period = rp, suffix = c("a", "m\u00B3"))


# applying a moving average before fitting
ng2 <- ng
ng2$flow <- ma(ng2$flow, n = 4)
tyears(ng2, dist = c("gum", "wei", "ln3", "pe3"), event = rp,
       plot = FALSE)

Loading required package: xts
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

Loading required package: lmom
Loading required package: lattice
Warning messages:
1: For fitting minima, a Weibull distribution with parameter 'zeta = 0' may be best. 
2: In evfit(x = minima, distribution = dist, zeta = zeta, check = check,  :
  Estimation of parameter zeta in the 'ln3' distribution resulted in a negative value (-2.91).  As this is not meaningful for discharges, parameter estimation was done with a forced lower bound of '0'. To override this behavior, consider setting the 'zeta' argument explicitly when calling the function.
             distribution
return period      gum      wei      ln3      pe3
          1.5 4.309692 4.424127 4.372547 4.407726
          5   3.355655 3.329617 3.351916 3.355020
          10  3.107806 3.008854 3.057580 3.028993
          100 2.628128 2.445523 2.458004 2.329147
List of 11
 $ freq.zeros   : num 0
 $ parameters   :List of 4
  ..$ gum: Named num [1:2] 3.685 0.692
  .. ..- attr(*, "names")= chr [1:2] "xi" "alpha"
  ..$ wei: Named num [1:3] 2.08 2.26 2.54
  .. ..- attr(*, "names")= chr [1:3] "zeta" "beta" "delta"
  ..$ ln3: Named num [1:3] 0 1.385 0.209
  .. ..- attr(*, "names")= chr [1:3] "zeta" "mu" "sigma"
  ..$ pe3: Named num [1:3] 4.084 0.854 0.364
  .. ..- attr(*, "names")= chr [1:3] "mu" "sigma" "gamma"
 $ lmom         : Named num [1:4] 4.0844 0.4797 0.0594 0.1565
  ..- attr(*, "names")= chr [1:4] "l_1" "l_2" "t_3" "t_4"
 $ values       : num [1:37(1d)] 3.34 4.8 5.02 4.82 3.21 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr [1:37] "1964" "1965" "1966" "1967" ...
 $ is.censored  : logi FALSE
 $ extreme      : chr "minimum"
 $ zeta         : NULL
 $ estimates    : num [1:37, 1:4] 3.36 4.71 5.27 4.82 3.24 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$     : NULL
  .. ..$ dist: Named chr [1:4] "gum" "wei" "ln3" "pe3"
  .. .. ..- attr(*, "names")= chr [1:4] "gum" "wei" "ln3" "pe3"
 $ p.value      : Named num [1:4] 1 1 1 1
  ..- attr(*, "names")= chr [1:4] "gum" "wei" "ln3" "pe3"
 $ r.squared    : Named num [1:4] 0.971 0.978 0.982 0.981
  ..- attr(*, "names")= chr [1:4] "gum" "wei" "ln3" "pe3"
 $ T_Years_Event: num [1:4, 1:4] 4.31 3.36 3.11 2.63 4.42 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ return period: chr [1:4] "1.5" "5" "10" "100"
  .. ..$ distribution : chr [1:4] "gum" "wei" "ln3" "pe3"
 - attr(*, "class")= chr [1:2] "evfit" "list"
Values: num [1:37(1d)] 3.34 4.8 5.02 4.82 3.21 ...


L-Moments:
       l_1        l_2        t_3        t_4 
4.08443243 0.47967868 0.05942792 0.15648535 

Fitted Parameters of the Distribution:
gum    xi:     3.68498,  alpha:  0.69203
wei    zeta:   2.07519,  beta:   2.26359,  delta:   2.5411
ln3    zeta:         0,  mu:     1.38536,  sigma: 0.208916
pe3    mu:     4.08443,  sigma: 0.853741,  gamma: 0.364278
             distribution
return period      gum      wei      ln3      pe3
          1.5 4.433022 4.549116 4.497004 4.531669
          5   3.428910 3.400051 3.424759 3.426952
          10  3.168052 3.067906 3.116942 3.088750
          100 2.663195 2.493522 2.492250 2.370888
Warning message:
In evfit(x = minima, distribution = dist, zeta = zeta, check = check,  :
  Estimation of parameter zeta in the 'ln3' distribution resulted in a negative value (-2.51).  As this is not meaningful for discharges, parameter estimation was done with a forced lower bound of '0'. To override this behavior, consider setting the 'zeta' argument explicitly when calling the function.