R/ERI.R

In alasdairhunter/MagicWP7: Package for computing climate extremes from the CMIP data.

#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))
}