Combined Return Level Plots in R

Hey there,

I have received multiple requests over the last few months to provide information on how I created the combined return level plots when using both annual maxima series and partial duration series. I have put off the evil hour to provide these scripts since I intended to come up with a nice ggplot solution.

As luck would have it I just needed to create another one of these combined return level plots for a conference poster, and I really shuddered at the thought of digging up my old R scripts I had used for that occasion. Even though it immediately turned out that my fear was not clearly not unfounded, for lack of time I have not been able to rework this behemoth of code (which is based on plot.fevd() from the extRemes package)to make it as tidy as initially intended. However, I think I have at least been able to restructure the basic plotting function rlplot() in a way that it might be understandable to other users. I might completely rework this function by building it from scratch one day, since the current version is an extremely verbose and somewhat very hacky solution. But it works, and that should be sufficient for now.

Since the function is somewhat cumbersome, I have provided it via our GitLab repository.

Attached is a short example of application, based on the same precipitation data set I have already used for the other posts on this topic. The functions required to run the code are read_ehyd() and rlplot().

# load required packages

# get data from eHYD
ehyd_url <- ""
precipitation_xts <- read_ehyd(ehyd_url)
precipitation_xts <- na.omit(precipitation_xts)

# derive AMS for maximum precipitation
ams <- apply.yearly(precipitation_xts, max)
ams <- as.vector(ams)

# annual maxima series
gev_mle <- fevd(ams, method = "MLE", type="GEV")
gev_lmom <- fevd(ams, method = "Lmoments", type="GEV")
gev_gmle <- fevd(ams, method = "GMLE", type="GEV")
gev_bpe <- fevd(ams, method = "Bayesian", type="GEV")

# partial duration series
u <- 40
precipitation_xts <- as.vector(precipitation_xts)
gp_mle <- fevd(precipitation_xts, method = "MLE", type="GP", threshold=u)
gp_lmom <- fevd(precipitation_xts, method = "Lmoments", type="GP", threshold=u)
gp_gmle <- fevd(precipitation_xts, method = "GMLE", type="GP", threshold=u)
gp_bpe <- fevd(precipitation_xts, method = "Bayesian", type="GP", threshold=u)

main_title = paste("Return Level Plot of Daily Precipitation Totals at Bärnkopf")
rlplot(GEV_MLE = gev_mle,
       GEV_LMOM = gev_lmom,
       GEV_GMLE = gev_gmle,
       GEV_BPE = gev_bpe,
       GP_MLE = gp_mle,
       GP_LMOM = gp_lmom,
       GP_GMLE = gp_gmle,
       GP_BPE = gp_bpe,
       plottitle = main_title,
       unit = "[mm/d]",
       y_seq = seq(0, 180, 10),
       do_ci = FALSE)
About This Author

Matthias studied Environmental Information Management at the University of Natural Resources and Life Sciences Vienna and holds a PhD in environmental statistics. The focus of his thesis was on the statistical modelling of rare (extreme) events as a basis for vulnerability assessment of critical infrastructure. He is working at the Austrian national weather and geophysical service (ZAMG) and at the Institute of Mountain Risk Engineering at BOKU University. He currently focuses the (statistical) assessment of adverse weather events and natural hazards, and disaster risk reduction. His main interests are statistical modelling of environmental phenomena as well as open source tools for data science, geoinformation and remote sensing.


You can post comments in this post.

  • Hi Matthias,

    Nice write up. I am new to GIS and I have data which is unclassified Lidar points . I have to classify them into grounds and vegetation classes. How can I go for it ?

    Atul Shukla 3 years ago Reply

    • Hi,

      this is actually a totally different topic.
      A simple internet search reveals interesting tutorials and blog posts on this topic, e.g.:

      Matthias 3 years ago Reply

      • Hi

        Thank you for your interesting blog

        When I ran the code I get this error.

        > # get data from eHYD
        > ehyd_url precipitation_xts precipitation_xts
        > # derive AMS for maximum precipitation
        > ams ams
        > # annual maxima series
        > gev_mle gev_lmom gev_gmle gev_bpe
        > # partial duration series
        > u precipitation_xts gp_mle gp_lmom gp_gmle gp_bpe <- fevd(precipitation_xts, method = "Bayesian", type="GP", threshold=u)
        Error in is.element("formula", class(x)) :
        object 'precipitation_xts' not found
        could please tell me how to solve that?

        Mustafa 1 year ago Reply

        • Hi,

          apparently, there is no object called precipitation_xts. Have you correctly specified this object?


          Matthias 11 months ago Reply

  • Hi Matthias,
    I have tried the scripts in Rstudio but they are not working. Can you help on this regard.

    Temesgen 2 years ago Reply

    • Hi,

      to be honest, my code for the combined return level plot is quite a mess. Its basically just a really ugly hack of the extRemes function plot.fevd() to produce the desired output once. It might be easier to simply use ggplot to construct the plot from scratch.
      Apart from that, please provide a more specific description of the problem. Otherwise, it is difficult to provide any help.

      Matthias 1 year ago Reply

Post A Reply