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
# load required packages library(extRemes) library(xts) # get data from eHYD ehyd_url <- "http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=107540&file=2" 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)