4  Standardized Precipitation Evapotranspiration Calculation

4.1 Objectif

Ce document vise à calculer l’indice Standardized Precipitation Evapotranspiration Index (SPEI), une des variables de contrôle intégrées dans l’analyse. L’ajout de cette variable permet de réduire la variabilité inexpliquée du modèle et estimer plus précisément l’effet du traitement. Nous calculons l’indice à partir d’une référence à long terme de 1981 à 2021, en utilisant les données mensuelles de précipitation, de température minimale et maximale provenant de Worlclim. Le calcul du SPEI est effectué en plusieurs étapes en utilisant la méthode de Hargreaves ( Modified-Hargreaves methods) définie par Droogers et Allen Droogers and Allen (2002) pour estimer l’évapotranspiration potentielle (PET), puis en intégrant ce calcul dans un indice global de sécheresse.

4.2 Méthode

Nous avons chargé le dataframe contenant les données sur la précipitation, la température minimale et la température maximale pour chaque clusters. A partir des coordonnées géométriques de chaque clusters, nous calculons la latitude, qui nous servira de base de calcul pour estimer l’évapotranspiration potentielle (PET). Le calcul de la latitude garantira que l’estimation du PET est spécifique aux conditions géographiques locales. Nous récupérons ensuite les séries temporelles mensuelles des variables climatiques de chaque cluster: température minimale (tmin), température maximale (tmax), précipitation (prec). A travers la formule de Hargreaves, nous estimons le PET et le water balance (c’est-à-dire la différence entre la précipitation et le PET). Le water balance est converti en série temporelle, à partir de laquelle l’indice SPEI est estimé sur une échelle de 12 mois, en utilisant une période de référence de 1981 à 2021. Pour une analyse plus globale, on agrège les résultats à l’échelle annuelle en calculant la moyenne annuelle des valeurs de SPEI par année pour chaque cluster.

4.3 Calcul de SPEI mensuel

On extrait la première ligne du dataframe pour calculer le SPEI mensuel sur un seul cluster. L’analyse d’un seul cluster permet de valider les étapes du calcul, mais aussi d’appréhender la structure des données et les résultats.

library(SPEI) # Calcul de l'indice SPEI
library(tidyverse) # Manipulation et visualisation des données
library(sf) # Analyse des données spatiales
library(tmap) # Analyse cartographique
library(labelled) # Manipulation des labels
library(tibbletime) # Manipulation des données temporelles 
library(zoo) # Manipulation des données temporelles 
library(readr) # Lecture des données de texte rectangulaires
library(ggplot2) # visualisation

# Load data 
spatial_covars <- read_rds("data/derived/spatial_covars_staggered.rds")

# Function to compute SPEI
compute_spei_annual <- function(tmin_tbl, tmax_tbl, prec_tbl, lat_deg) {
  # tmin_tbl/tmax_tbl/prec_tbl: tibbles avec colonnes `datetime` (Date) et `value` (num)
  d_tmin <- tibble(date = tmin_tbl$datetime, tmin = tmin_tbl$value)
  d_tmax <- tibble(date = tmax_tbl$datetime, tmax = tmax_tbl$value)
  d_prec <- tibble(date = prec_tbl$datetime, prec = prec_tbl$value)

  d_merged <- reduce(list(d_tmin, d_tmax, d_prec), left_join, by = "date") %>%
    arrange(date)

  d_clean <- drop_na(d_merged)  # supprime lignes avec NA

  # Si séries trop courtes, renvoyer squelette 1981:2021 en NA
  if (nrow(d_clean) < 12) {
    return(tibble(year = 1981:2021, spei_mean = NA_real_))
  }

  # PET (Hargreaves), bilan hydrique et SPEI mensuel
  pet <- hargreaves(Tmin = d_clean$tmin,
                    Tmax = d_clean$tmax,
                    Pre  = d_clean$prec,
                    lat  = lat_deg)

  wb <- d_clean$prec - pet

  wb_ts <- ts(wb,
              start = c(year(min(d_clean$date)), month(min(d_clean$date))),
              frequency = 12)

  spei_obj <- spei(wb_ts,
                   scale = 12,
                   ref.start = c(1981, 1),
                   ref.end   = c(2021, 12))

  tibble(datetime = d_clean$date,
         spei      = as.numeric(spei_obj$fitted)) %>%
    filter(datetime >= as.Date("1981-01-01"),
           datetime <= as.Date("2021-12-31")) %>%
    mutate(year = year(datetime)) %>%
    group_by(year) %>%
    summarise(spei_mean = mean(spei, na.rm = TRUE), .groups = "drop") %>%
    complete(year = 1981:2021, fill = list(spei_mean = NA_real_))
}



# Latitude géodésique depuis la géométrie (centroïde)
spatial_covars_spei <- spatial_covars %>%
  mutate(lat = st_coordinates(st_centroid(geometry))[, 2])

row1 <- spatial_covars_spei[1, ]

spei_tbl_one <- compute_spei_annual(
  tmin_tbl = row1$temperature_min_wc[[1]],
  tmax_tbl = row1$temperature_max_wc[[1]],
  prec_tbl = row1$precipitation_wc[[1]],
  lat_deg  = row1$lat
)

[1] "Calculating reference evapotranspiration using the Hargreaves method. Using latitude (`lat`) to estimate extraterrestrial radiation. Using precipitation data, following Droogers and Allen (2002). Checking for missing values (`NA`): all the data must be complete. Input type is vector. Assuming the data are monthly time series starting in January, all regular (non-leap) years."
[1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 12. Using kernel type 'rectangular', with 0 shift. Fitting the data to a log-Logistic distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using a user-specified reference period. Input type is tsvector. Time series spanning janv. 1980 to oct. 2021, with frequency = 12."

# Graphique 
spei_2017 <- spei_tbl_one %>% 
  filter(year == 2017)

ggplot(spei_tbl_one, aes(x = year, y = spei_mean)) +
  geom_col(aes(y = pmax(0, -spei_mean)), alpha = 0.25) +  # barres pour sécheresse (optionnel)
  geom_line(linewidth = 0.8) +
  geom_vline(xintercept = 2017, color = "red", linetype = "dotted", linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_label(
    data = spei_2017, 
    aes(x = 2017, y = spei_mean, 
        label = paste0("Année: 2017\nSPEI: ", round(spei_mean, 2))),
    nudge_x = 1, 
    nudge_y = 0.3,
    fill = "white",
    color = "black",
    label.size = 0.4, 
    label.padding = unit(0.2, "lines")
  ) +
  
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Annual SPEI (scale=12) – cluster de démonstration",
       x = "Année", y = "SPEI (moyenne annuelle)") +
  theme_minimal()

Le graphique illustre l’évolution annuelle de l’indice SPEI moyen pour un cluster géographique entre 1981 et 2021. La courbe bleue représente la tendance annuelle du SPEI moyen, mettant en exergue les variations interannuelles du climat dans la zone étudiée. Les points rouges indiquent les valeurs du SPEI moyen pour chaque année. La ligne horizontale noire en pointillées, placées à y = 0, sert de référence pour les conditions climatiques neutres.

On observe qu’entre 1980 et 2000, la majorité des années a un SPEI supérieur à 0, correspondant à des conditions normales ou légèrement humides, avec peu d’épisodes de sécheresse significative selon la classification de Mckee et al (1993?). En revanche, à partir des années 2004, on obtient des SPEI inférieur à 0, marquant une augmentation de la fréquence des épisodes de sécheresse. En 2017, le SPEI a une valeur particulièrement basse (SPEI = -2.40012658), qui d’après la classification de Vicente-Serrano et al. (2010?) caractérise une sécheresse très sévère.

Nous appliquons la même procédure pour l’ensemble des clusters étudiés.

# SPEI for all clusters 
quiet <- function(x) {
  tmp <- capture.output(res <- x, file = NULL)
  res
}

spatial_covars_spei <- spatial_covars_spei %>%
  mutate(
    spei_wc = pmap(
      list(temperature_min_wc, temperature_max_wc, precipitation_wc, lat),
      ~ quiet(compute_spei_annual(..1, ..2, ..3, ..4))
    )
  )

spei_df <- spatial_covars_spei %>%
  select(DHSCLUST, spei_wc) %>%
  unnest(spei_wc)

ggplot(spei_df, aes(x = year, y = spei_mean, group = DHSCLUST)) +
  geom_line(alpha = 0.35) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Annual SPEI by cluster (1981–2021)",
       x = "Année", y = "SPEI (moyenne annuelle)") +
  theme_minimal()

Chaque ligne bleue correspond à un groupe de ménages différent. On constate que ces lignes ne se superposent pas parfaitement, ce qui suggère que les conditions d’humidité varient d’une région à l’autre. Un indice SPEI positif indique des périodes plus humides, tandis qu’un indice SPEI négatif correspond à des périodes plus sèches. Cela montre que, même si les ménages sont regroupés, ils n’ont pas vécu les mêmes conditions climatiques. On observe néanmoins des baisses significatives de l’indice SPEI, notamment des épisodes de sécheresse très sévère en 1992 et entre 2017 et 2020.

4.4 Transformation de la colonne spei_wc en format wide

# Chaque élément spei_wc devient une table {datetime, variable, unit, value}
spatial_covars_spei <- spatial_covars_spei %>%
  mutate(
    spei_wc = map(spei_wc, ~ .x %>%
      mutate(
        datetime = as.Date(paste0(.data$year, "-01-01")),
        variable = "spei_scale12_mean",
        unit     = "annual",
        value    = .data$spei_mean
      ) %>%
      select(datetime, variable, unit, value))
  )

# Sauvegarde (cohérente avec le reste de tes scripts)
spatial_covars_spei_df <- as.data.frame(spatial_covars_spei)
write_rds(spatial_covars_spei_df, "data/derived/spatial_covars_spei_staggered.rds")
cat("SPEI (annuel, 1981–2021) enregistré dans data/derived/spatial_covars_spei_staggered.rds\n")

SPEI (annuel, 1981–2021) enregistré dans data/derived/spatial_covars_spei_staggered.rds

Droogers, Peter, and Richard G. Allen. 2002. “Estimating Reference Evapotranspiration Under Inaccurate Data Conditions.” Irrigation and Drainage Systems 16 (1): 3345. https://doi.org/10.1023/A:1015508322413.