#' @title Compute the Standardized Precipitation Index.
#'
#' @description
#' \code{spi_week} computes the weekly spi index with gamlss models.
#' It can also with work with precipitation series containing NA, 0 or
#' only non-zero values.
#'
#' @details
#' details.
#'
#' @param rain Precipitation level.
#' @param tscale Time-scale to compute the SPI in week units.
#' @param period Period (e.g. 53 weeks) defined to model the seasonal effect.
#' @param package Package to use. \code{gamlss} for backfitting and \code{bamlss} for baysian gamlss models.
#' @param plot Logical value indicating if a plot should be draw or not.
#'
#' @return dataframe consisting of spi, mu, sigma and pzero.
#'
#' @author Erick A. Chacon-Montalvan
#'
#' @importFrom mgcv gam predict.gam
#' @importFrom gamlss pbc gamlss
#' @importFrom utils getFromNamespace
#' @importFrom graphics par points lines
#'
#' @export
spi_week <- function(rain, tscale = 1, period = 365 / 7, package = "gamlss", plot = FALSE) {
# require(mgcv)
# Compute moving average of rain based on time-scale.
rain <- runmean(rain, tscale)
# Create dataframe.
data <- data.frame(rain, weeks = 1:length(rain))
data <- transform(data, no_rain = (rain == 0) * 1, rain_level = rain)
data <- within(data, rain_level[rain == 0] <- NA)
data <- transform(data, weeks2 = weeks %% period) # for seasonal trend
# Modelling probability of no rain.
gam0 <- mgcv::gam(no_rain ~ s(weeks2, bs = "cc"), binomial("logit"), data)
data$pzero <- mgcv::predict.gam(gam0, data, type = "response")
if (package == "gamlss") {
# require(gamlss)
predict.gamlss <- getFromNamespace("predict.gamlss", "gamlss")
form_mu <- rain_level ~ gamlss::pbc(weeks2)
form_sg <- ~ gamlss::pbc(weeks2)
data0 <<- na.omit(data)
lss <- gamlss::gamlss(form_mu, sigma.formula = form_sg, data = data0, family = GA,
trace = FALSE)
mu <- predict.gamlss(lss, what = "mu", newdata = data, type = "response")
sigma <- predict.gamlss(lss, what = "sigma", newdata = data, type = "response")
data$mu <- mu
data$sigma <- 1 / sigma ^ 2
# Notes:
# 1) predict only works if data0 exists in the global env.
# 2) stepGAIC only works if form_mu and form_sg exits in the global env.
# Solution: use <<- instead of <- for global assignment.
} else {
# # require(bamlss)
# # Modelling meanrain to select harmonic terms.
# lm0 <- lm(paste0("rain_level ~ 1", terms), data)
# lm1 <- step(lm0, trace = 0)
# # lm1 <- step(lm0, trace = 0, k = log(nrow(data)))
#
# # Distributional model for rain_level.
# # form_lss <- list(update(formula(lm1), rain_level ~ .),
# # update(formula(lm1), sigma ~ .))
# form0 <- as.character(formula(lm1))
# form0 <- as.formula(paste(form0[2], form0[1], form0[3]))
# # form0 <- formula(lm1)
# form_lss <- list(update(form0, rain_level ~ .),
# update(form0, sigma ~ .))
#
# # form_lss <- list(rain_level ~ 1,
# # sigma ~ 1)
# lss <- bamlss(form_lss, family = "gamma", data = data)
# data[c("mu", "sigma")] <- predict(lss, data, type = "parameter")
}
# Compute the spi and additional variables.
data <- transform(data, cdf_gamma = pgamma(rain, sigma, sigma / mu))
data <- transform(data, cdf_mix = pzero + (1-pzero) * cdf_gamma)
data <- transform(data, spi = qnorm(data$cdf_mix),
q025 = qgamma((0.025 - pzero) / (1-pzero), sigma,
sigma / mu),
q975 = qgamma((0.975 - pzero) / (1-pzero), sigma,
sigma / mu))
data <- within(data, {q025[is.na(q025)] <- 0;
q975[is.na(q975)] <- 0})
if (plot) {
opar <- par(mfrow = c(1, 2))
plot(data$weeks, data$spi)
points(data$weeks[abs(data$spi) > qnorm(0.975)],
data$spi[abs(data$spi) > qnorm(0.975)], col = 2, lwd = 2)
plot(data$weeks, data$rain)
points(data$weeks[abs(data$spi) > qnorm(0.975)],
data$rain[abs(data$spi) > qnorm(0.975)], col = 2, lwd = 2)
lines(data$weeks, data$mu)
lines(data$weeks, data$q025)
lines(data$weeks, data$q975)
par(opar)
}
return(subset(data, select = c(spi, mu, sigma, pzero)))
}
#' @title Identify floods and droughts based on SPI.
#'
#' @description
#' \code{find_flood_drought} identifies floods and droughts for a vector series of
#' standardized precipitation values.
#'
#' @details
#' Floods (droughts) are identified by a subset of negative values where at least one value
#' exceeds -1 (1).
#'
#' @param spi A vector series of standardized precipitation values.
#'
#' @return return.
#'
#' @author Erick A. Chacon-Montalvan
#'
#' @export
find_flood_drought <- function (spi) {
# Convert spi to string to detect extreme events.
spi <- (spi > -1) + (spi > 0) + (spi > 1)
spi[is.na(spi)] <- 9
spichar <- spi %>% paste(collapse = "")
# Detect flood and drought.
flood_exp <- gregexpr("[2-3]*3+[2-3]*", spichar)
drought_exp <- gregexpr("[0-1]*0+[0-1]*", spichar)
flood <- function(n) strrep("8", n)
drought <- function(n) strrep("7", n)
regmatches(spichar, flood_exp) <-
Map(flood, lapply(regmatches(spichar, flood_exp), nchar))
regmatches(spichar, drought_exp) <-
Map(drought, lapply(regmatches(spichar, drought_exp), nchar))
# Convert to vector.
spichar <- substring(spichar, 1:nchar(spichar), 1:(nchar(spichar)))
spi <- as.numeric(spichar)
spi[spi == 9] <- NA
spi[!(spi %in% c(7:8, NA))] <- 0
spi <- factor(spi, c(0, 7, 8), c("normal", "drought", "flood"))
return(spi)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.