#ERI function
ERI <- function(proj, ref, index = "temp", period = 1, quantile = 0.9, threshold_upper = NULL, threshold_lower = 1, season_length = 365) {
if (!is.null(threshold_upper) && !is.null(threshold_lower)) {
stop('Only one of "threshold_upper" and "threshold_lower" can be specified')
}
if (!is.null(threshold_lower) && index == "drought") {
threshold <- -threshold_lower
ref <- -ref
proj <- -proj
}
#Calculate the number of seasons
ref.years <- length(ref) / season_length
proj.years <- length(proj) / season_length
#Calculate the season_length after the rolling average has been calculated
season_length <- season_length - (period - 1)
#Calculate the rolling sum for moving window of length 'period'
#for the reference and projections
if (index == "temp") {
ref.period <- rollsum(ref, k = period)
proj.period <- rollsum(proj, k = period)
#Calculate the 90th percentile from the reference
ref90.percentile <- quantile(ref.period, quantile)
#Calculate the annual index (percentage of days exceeding the 90th percentile of the reference period each year)
ref.annual.index <- proj.annual.index <- c()
for (yr in 1 : ref.years) {
ref.annual <- ref.period[(yr * season_length - (season_length - 1)) : (yr * season_length)]
ref.annual.index[yr] <- length(ref.annual[ref.annual > ref90.percentile]) / length(ref.annual) * 100
}
for (yr in 1 : proj.years) {
proj.annual <- proj.period[(yr * season_length - (season_length - 1)) : (yr * season_length)]
proj.annual.index[yr] <- length(proj.annual[proj.annual > ref90.percentile]) / length(proj.annual) * 100
}
#The interannual standard deviation and the reference value for the reference period
index.sd <- sd(ref.annual.index)
reference <- (1 - quantile) * 100
} else if (index == "flood") {
ref.period <- rollsum(ref, k = period)
proj.period <- rollsum(proj, k = period)
#Calculate the annual maximum
ref.annual.index <- proj.annual.index <- c()
for (yr in 1 : ref.years) {
ref.annual.index[yr] <- max(ref.period[(yr * season_length - (season_length - 1)) : (yr * season_length)])
}
for (yr in 1 : proj.years) {
proj.annual.index[yr] <- max(proj.period[(yr * season_length - (season_length - 1)) : (yr * season_length)])
}
reference <- mean(ref.annual.index)
index.sd <- sd(ref.annual.index)
}
if (index == "drought") {
ref.annual.index <- proj.annual.index <- c()
for (yr in 1 : ref.years) {
ref.annual.index[yr] <- max(c(unlist(lapply(clusters(ref[(yr * season_length - (season_length - 1)) : (yr * season_length)], u = threshold, keep.names = F), function(x) length(x)))))
}
for (yr in 1 : proj.years) {
proj.annual.index[yr] <- max(c(unlist(lapply(clusters(proj[(yr * season_length - (season_length - 1)) : (yr * season_length)], u = threshold, keep.names = F), function(x) length(x)))))
}
reference <- mean(ref.annual.index)
index.sd <- sd(ref.annual.index)
}
#The standardised index
invisible(result <- list(reference.index = ((ref.annual.index - reference) / index.sd), projected.index = ((proj.annual.index - reference) / index.sd), ref.index.sd = index.sd))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.