---
title: "Introduction to ClimateAnomaly: Computing Climate Anomaly Indices"
author: "Pushkar Bora"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Introduction to ClimateAnomaly: Computing Climate Anomaly Indices}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  fig.width  = 7,
  fig.height = 4
)
```

## Overview

**ClimateAnomaly** is a general-purpose R package for computing standardised
rainfall anomaly indices and temperature anomalies from any monthly (or
annual) climate dataset.
All indices are computed **group-wise** — typically per *location × calendar
month* — so that each value reflects a departure from that location's own
long-term climatology for the same time of year.

### Indices implemented

| Index | Full name | Reference |
|-------|-----------|-----------|
| **RAI**  | Rainfall Anomaly Index | Van Rooy (1965) |
| **MRAI** | Modified Rainfall Anomaly Index | Ekpetere (2025) |
| **SRAI** | Standardised Rainfall Anomaly Index | Ekpetere (2025) |
| **SSPI** | Simplified Standardised Precipitation Index | Raziei (2021) |
| **SPI**  | Standardised Precipitation Index | McKee et al. (1993) |
| **TEMP** | Temperature anomaly (z-score) | — |

---

## 1  Quick start with the built-in dataset

```{r load-pkg}
library(ClimateAnomaly)
```

The package ships a sample dataset `india_monthly` — monthly climate records
for **six Indian states** (1976–2026) spanning arid, tropical, Himalayan,
semi-arid, high-rainfall, and peninsular climate zones.

```{r inspect-data}
data(india_monthly)
dim(india_monthly)
head(india_monthly[, 1:8])
```

Compute all six indices in one call:

```{r compute-all}
result <- calc_all_indices(
  data       = india_monthly,
  group_cols = c("state", "month_num"),  # group = state × calendar month
  rain_col   = "rainfall_mm",
  tmax_col   = "tmax_mean_c",
  tmin_col   = "tmin_mean_c"
)

# New columns added
setdiff(names(result), names(india_monthly))
```

Inspect a single state:

```{r show-maharashtra}
mh <- subset(result, state == "Maharashtra" & year >= 2015,
             select = c("year", "month_num", "rainfall_mm",
                        "RAI", "RAI_class", "SPI", "SPI_class",
                        "tmax_anomaly", "tmax_anom_class"))
print(head(mh, 24))
```

---

## 2  Using your own data file

ClimateAnomaly imposes **no fixed column names**.  Supply the actual column
names from your file as arguments.

### 2.1  Loading a CSV

```r
my_data <- read.csv("my_climate_data.csv", stringsAsFactors = FALSE)
head(my_data)
```

Assume the file has columns `district`, `yr`, `mon`, `precip_mm`,
`tmax_deg`, `tmin_deg`.

```r
my_result <- calc_all_indices(
  data       = my_data,
  group_cols = c("district", "mon"),   # group by district AND calendar month
  rain_col   = "precip_mm",
  tmax_col   = "tmax_deg",
  tmin_col   = "tmin_deg"
)
```

### 2.2  Loading an Excel file

```r
# install.packages("readxl")   # if not yet installed
library(readxl)
my_xl <- as.data.frame(read_excel("my_data.xlsx", sheet = "Monthly"))

my_result <- calc_all_indices(
  data       = my_xl,
  group_cols = c("Station", "Month"),
  rain_col   = "Rainfall_mm",
  tmax_col   = "Tmax_C",
  tmin_col   = "Tmin_C"
)
```

### 2.3  Annual data (no month grouping)

For annual totals, group only by location:

```r
annual_result <- calc_all_indices(
  data       = my_annual_data,
  group_cols = "station",   # or "state", "grid_id", …
  rain_col   = "annual_rain_mm",
  tmax_col   = NULL,   # skip temperature
  tmin_col   = NULL,
  indices    = c("RAI", "MRAI", "SPI")
)
```

### 2.4  Using the bundled full India dataset

```r
src <- system.file("extdata", "state_weather_monthly.csv",
                   package = "ClimateAnomaly")
india_full <- read.csv(src, stringsAsFactors = FALSE)
full_result <- calc_all_indices(india_full,
                                group_cols = c("state","month_num"),
                                rain_col   = "rainfall_mm",
                                tmax_col   = "tmax_mean_c",
                                tmin_col   = "tmin_mean_c")
```

---

## 3  Running individual index functions

Each index can also be called directly on a **single numeric vector**
(one group at a time).

### 3.1  RAI — Rainfall Anomaly Index

```{r rai-direct}
# Rainfall series for Kerala July (all available years)
kl_jul <- subset(india_monthly, state == "Kerala" & month_num == 7)$rainfall_mm
rai_vals <- calc_rai(kl_jul, n_extreme = 10L)
summary(rai_vals)
classify_rai(rai_vals)
```

### 3.2  MRAI and SRAI

```{r mrai-srai}
mrai_vals <- calc_mrai(kl_jul)
srai_vals <- calc_srai(kl_jul)
data.frame(year = subset(india_monthly, state=="Kerala" & month_num==7)$year,
           RAI  = round(rai_vals, 2),
           MRAI = round(mrai_vals, 2),
           SRAI = round(srai_vals, 2))
```

### 3.3  SPI — Standardised Precipitation Index

```{r spi-direct}
spi_vals <- calc_spi(kl_jul)
summary(spi_vals)
```

### 3.4  SSPI — Simplified SPI (skewness-based)

```{r sspi-direct}
sspi_vals <- calc_sspi(kl_jul)
# Compare SPI vs SSPI
cor(spi_vals, sspi_vals, use = "complete.obs")
```

### 3.5  Temperature anomaly

```{r temp-anom}
kl_tmax <- subset(india_monthly, state == "Kerala" & month_num == 6)$tmax_mean_c
z <- calc_temp_anomaly(kl_tmax)
summary(z)
classify_temp_anomaly(z)
```

---

## 4  Selecting a subset of indices

Use the `indices` argument to compute only what you need:

```{r subset-indices}
rain_only <- calc_all_indices(
  data       = india_monthly,
  group_cols = c("state", "month_num"),
  rain_col   = "rainfall_mm",
  tmax_col   = NULL,
  tmin_col   = NULL,
  indices    = c("RAI", "SPI"),
  classify   = FALSE   # skip classification labels
)
names(rain_only)
```

---

## 5  Drought classification

```{r class-summary}
# Frequency table of SPI classes across all states/months
tab <- table(result$SPI_class)
sort(tab, decreasing = TRUE)
```

```{r rai-class-by-state}
# RAI class distribution per state
with(result, tapply(RAI_class, state, function(x) {
  sort(table(x), decreasing = TRUE)[1:3]
}))
```

---

## 6  Simple time-series plot

```{r plot-spi, fig.cap = "SPI for Kerala June (1976–2026)"}
kl_jun <- subset(result, state == "Kerala" & month_num == 6)
kl_jun <- kl_jun[order(kl_jun$year), ]

plot(kl_jun$year, kl_jun$SPI,
     type = "h", col = ifelse(kl_jun$SPI >= 0, "steelblue", "firebrick"),
     lwd  = 2,
     xlab = "Year", ylab = "SPI",
     main = "Kerala – June SPI (1976–2026)")
abline(h = c(-1, 1), lty = 2, col = "grey50")
abline(h = 0, lwd = 1)
```

---

## 7  Minimum data requirements

| Index | Min non-NA observations per group | Notes |
|-------|------------------------------------|-------|
| RAI   | 4 | `n_extreme` capped at floor(n/2) |
| MRAI  | 4 | Uses full-series min/max |
| SRAI  | 4 | Uses full-series SD; output clamped to [−3, +3] |
| SSPI  | 4 | Uses 5th/95th percentiles; PZP-branch for dry months |
| SPI   | 4 | Gamma MLE; falls back to z-score if < 4 non-zero |
| TEMP  | 4 | z-score; returns 0 if SD = 0 |

Groups with fewer than 4 non-missing values return `NA`.

---

## References

Van Rooy, M.P. (1965). A rainfall anomaly index independent of time and
space. *Notos*, 14, 97–104.

Ekpetere, K.O. (2025). Development and evaluation of the modified and
standardized rainfall anomaly indices for extreme variability analysis.
*Journal of Environmental Management*, 375, 124160.
https://doi.org/10.1016/j.jenvman.2025.124160

Raziei, T. (2021). Revisiting the Rainfall Anomaly Index to serve as a
Simplified Standardized Precipitation Index. *Journal of Hydrology*, 602,
126761. https://doi.org/10.1016/j.jhydrol.2021.126761

Raziei, T. (2022). Improving the normalization procedure of the simplified
standardized precipitation index (SSPI) using Box-Cox transformation.
*Stochastic Environmental Research and Risk Assessment*, 37, 925–951.
https://doi.org/10.1007/s00477-022-02317-9

McKee, T.B., Doesken, N.J., & Kleist, J. (1993). The relationship of
drought frequency and duration to time scales. *Proceedings of the 8th
Conference on Applied Climatology*, AMS, pp. 179–184.