climatic-indices.Rmd

---
title: "Calculating basic climatic indices with data from easyclimate"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Calculating basic climatic indices with data from easyclimate}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r climatic-indices-1, include = FALSE}
library(knitr)
knitr::opts_chunk$set(collapse = TRUE, warning = FALSE, message = FALSE, 
                      fig.width = 6, fig.height = 4, cache = TRUE)
```

First, let's download daily climatic data for a specific location and store it in a dataframe.

```{r climatic-indices-2}

library(easyclimate)
library(tidyr)
library(dplyr)

coords <- matrix(c(-5.36, 37.40), ncol = 2)

daily <- get_daily_climate(coords,
                           period = 2001:2005,
                           climatic_var = c("Prcp","Tmin","Tmax"))  

```


Here we calculate the daily mean temperature:

```{r climatic-indices-3}

daily <- daily |> 
  mutate(Tmean = (Tmin + Tmax) / 2, #Daily mean temperature
         date = as.Date(date),
         month = format(date, format = "%m"),
         year = format(date, format = "%Y"))
```

<br>

## Average climatic values per site or time period

To calculate average temperatures by site or time period we can use `group_by` and `summarise` from `dplyr`, or `by` and `aggregate` from base R: 

```{r climatic-indices-4}

daily |> 
  group_by(ID_coords) |> 
  summarise(Tmin.site = mean(Tmin),
            Tmean.site = mean(Tmean),
            Tmax.site = mean(Tmax))

daily |> 
  group_by(year) |>  
  summarise(Tmin.year = mean(Tmin),
            Tmean.year = mean(Tmean),
            Tmax.year = mean(Tmax)) |> 
  kable(digits = 1)

daily |>  
  group_by(month) |> 
  summarise(Tmin.month = mean(Tmin),
            Tmean.month = mean(Tmean),
            Tmax.month = mean(Tmax)) |> 
  kable(digits = 1)


```


Similarly, you can calculate accumulated precipitation over different time periods: 

```{r climatic-indices-5}
daily |> 
  group_by(year) |>  
  summarise(sumPrec = sum(Prcp)) |> 
  kable(digits = 0)

daily |> 
  group_by(year, month) |> 
  summarise(sumPrec = sum(Prcp)) |>  
  pivot_wider(names_from = month, values_from = sumPrec) |> 
  kable(digits = 0)
```


<br>

## Climatic anomalies

Sometimes it can be useful to calculate the deviation of climate values from a reference value. For example, if we consider as the reference the mean value for the period 2001-2005, we can calculate the deviations per site and year as follows:

```{r climatic-indices-6}

daily |>    
  group_by(ID_coords, year) |>  
  summarise(avgTm = mean(Tmean)) |>  
  group_by(ID_coords) |>  
  mutate(refTm = mean(avgTm)) |>  #Reference mean temperature
  mutate(devTm = avgTm - refTm) |>  
  kable(digits = 1)

daily |>    
  group_by(ID_coords, year) |>  
  summarise(sumPrec = sum(Prcp)) |> 
  group_by(ID_coords) |> 
  mutate(refPrec = mean(sumPrec)) |>  #Reference precipitation
  mutate(devPrec =sumPrec - refPrec) |>  
  kable(digits = 0)

```

<br>

## Climatic events

Having daily values allows us to calculate climatic indices such as the number of consecutive days above a given temperature (e.g. heat waves or growing degree days), the number of spring frosts or the number of consecutive days without rain.

<br>

Let's start with an example on heat waves, or number of days with maximum temperatures above a given threshold:

```{r climatic-indices-7}

threshold <- 32  

daily |>  
  mutate(ID_coords = as.factor(ID_coords),
         year = factor(year)) |>  
  group_by(ID_coords, year) |>  
  mutate(event = ifelse(Tmax >= threshold, 1, 0), # NAs will be 0 
         start = c(1, diff(event) != 0),
         run = cumsum(start)) |>  
  filter(event == 1) |> 
  group_by(run, .add = TRUE) |>  
  mutate(length = n()) |>  
  group_by(ID_coords, year) |> 
  summarise(max_heatwave = max(length), # longest heatwave per year (in days)
            mean_heatwave = mean(length), # average length of heatwaves per year
            n_heatwave = n()) |>  # No. of events per year
  kable(digits = 1)

```

<br>
To calculate the number of days with spring frosts (minimum temperatures below 0) we can use the next code:

```{r climatic-indices-8}

spring.months <- c("03","04","05") # March to May

daily |>  
  mutate(ID_coords = as.factor(ID_coords),
         year = factor(year)) |>  
  filter(month %in% spring.months) |>  
  mutate(event = ifelse(Tmin < 0, 1, 0)) |>  # NAs will be 0 
  group_by(ID_coords, year) |>  
  mutate(n_frost = sum(event)) |>  # No. of days with minimum temperature < 0 per year
  filter(event == 1) |>  
  group_by(ID_coords, year, n_frost) |>  
  summarise(Tmin_frost_avg = mean(Tmin)/100) |>  # Mean minimum temperature of frost days
  ungroup() |>  
  complete(ID_coords, year, fill = list(n_frost = 0, Tmin_frost_avg = NA)) |>  
  kable()

```

<br>
A dry periods can be defined as a specific number of consecutive days with no precipitation. 

```{r climatic-indices-9}

threshold <- 30 #threshold in days to count extreme long events of no-rain

daily %>%  
  mutate(event = ifelse(Prcp < 0.01, 1, 0), # NAs will be 0 
         start = c(1, diff(event) != 0)) |>  
    group_by(ID_coords, year) |>  
    mutate(run = cumsum(start)) |> 
    filter(event == 1) |>  
    group_by(ID_coords, year, run) |>  
    summarise(length = n()) |>  
    ungroup(run) |> 
    summarise(max_days_norain = max(length), # Maximum length of periods of consecutive days without rain per year
              mean_days_norain = mean(length), 
              nextreme_events = sum(length > threshold)) %>% # No. of events per year
  ungroup() |>  
  complete(ID_coords, year, 
           fill = list(max_days_norain = NA, mean_days_norain = NA, nevents = 0))  |>  
  kable(digits = 1)

```


<br>

## Using ClimInd in combination with easyclimate

The package [`climInd`](https://cran.r-project.org/package=ClimInd) can be used to calculate multiple climatic indices from data obtained through `easyclimate`. For example:

<br>
Number of days with maximum temperature > 32ºC in summer (June-August):

```{r climatic-indices-10}
library(ClimInd)

Tmax <- as.vector(daily$Tmax)
names(Tmax) <- as.character(format(daily$date, format = "%m/%d/%Y"))

d32(Tmax)
```


<br>
Growing degree days:

```{r climatic-indices-11}
Tmean <- as.vector(daily$Tmean)
names(Tmean) <- as.character(format(daily$date, format = "%m/%d/%Y"))

gd4(data = Tmean, time.scale = "year") 
```

<br>
Growing season length:

```{r climatic-indices-12}
gsl(data = Tmean, time.scale = "year") 
```

<br>
Heavy precipitation days:

```{r climatic-indices-13}
Prec <- as.vector(daily$Prcp)
names(Prec) <- as.character(format(daily$date, format = "%m/%d/%Y"))

d50mm(data = Prec, time.scale = "month")
```

<br>

## Learn more

To learn more about downloading daily climatic data, please consult [this vignette](https://verughub.github.io/easyclimate/articles/points-df-mat-sf.html) for point coordinates, or [this other](https://verughub.github.io/easyclimate/articles/polygons-raster.html) if you need to extract data for several polygons or a region.