PACKAGES
PRECIPITATION, TEMPERATURE ET EPAISSEUR NEIGEUSE
NDVI
NAO
- Téléchargement des données depuis NOAA website
- Visualisation de la série temporelle NAO
INDICE PAR SAISONS
GENERATE DATAFRAME
- Long format
- Short format
CORRELATION DES INDICES
ACP :
CONCLUSION DES ANALYSES
CODE GOOGLE ENGINE

PACKAGES

library(data.table)  # to decompress .gz files
library(stringr)  # Manipulation de chaînes de caractères
library(future)   # Planification de calculs parallèles
library(furrr)    # Fonctions purrr compatibles avec future pour le parallélisme
library(httr)     # Pour les requêtes HTTP (téléchargement de fichiers)
library(terra)    # Gestion et traitement de données raster (NDVI par exemple)
library(mapview)  # Visualisation carte 
library(sf)       # Gestion et traitement de données raster
library(tidyverse)

PRECIPITATION, TEMPERATURE ET EPAISSEUR NEIGEUSE

Les données sont fournies par METEO FRANCE et sont téléchargeables directement depuis le site .gouv : https://www.data.gouv.fr/organizations/meteo-france/datasets. Les relevés quotidiens (Données climatologiques de base - quotidiennes) et mensuels (Données climatologiques de base - mensuelles) sont disponibles. Il faut ensuite chercher par numéro de département pour obtenir les données de l’ensemble des stations de météo de celui-ci. Les données de précipitations, températures et vents sont disponibles dans les fichiers dont le noms finissent par RR-T-Vent et les données de neige dans autres-paramètres.

Téléchargement des données depuis METEO FRANCE

Les données sont téléchargeables directement depuis R avec l’URL du jeu de données.

# URL
url_data_temperature_precipitation <- "https://www.data.gouv.fr/api/1/datasets/r/828ec510-e062-48e4-9370-0bd9a5f65a9c" 

url_data_temperature_precipitation_2025_2026 <- "https://www.data.gouv.fr/api/1/datasets/r/1fb641c2-fb6b-4fa1-8961-a41c349e8629"

url_data_snow <-
  "https://www.data.gouv.fr/api/1/datasets/r/94757a8e-fa4a-4608-b846-316502983c4a"

url_data_snow_2025_2026 <- "https://www.data.gouv.fr/api/1/datasets/r/98852d18-6580-47d4-9442-0bb3184896c8"

# TELECHARGEMENT ET DECOMPRESSION DES DONNEES
df_temperature_precipitation <- 
  rbind(data.table::fread(url_data_temperature_precipitation), data.table::fread(url_data_temperature_precipitation_2025_2026)) 

df_snow <-   rbind(data.table::fread(url_data_snow), data.table::fread(url_data_snow_2025_2026))

colnames(df_temperature_precipitation)

##  [1] "NUM_POSTE" "NOM_USUEL" "LAT"       "LON"       "ALTI"      "AAAAMMJJ"  "RR"        "QRR"       "TN"        "QTN"       "HTN"       "QHTN"      "TX"        "QTX"       "HTX"       "QHTX"     
## [17] "TM"        "QTM"       "TNTXM"     "QTNTXM"    "TAMPLI"    "QTAMPLI"   "TNSOL"     "QTNSOL"    "TN50"      "QTN50"     "DG"        "QDG"       "FFM"       "QFFM"      "FF2M"      "QFF2M"    
## [33] "FXY"       "QFXY"      "DXY"       "QDXY"      "HXY"       "QHXY"      "FXI"       "QFXI"      "DXI"       "QDXI"      "HXI"       "QHXI"      "FXI2"      "QFXI2"     "DXI2"      "QDXI2"    
## [49] "HXI2"      "QHXI2"     "FXI3S"     "QFXI3S"    "DXI3S"     "QDXI3S"    "HXI3S"     "QHXI3S"    "DRR"       "QDRR"

La descrpition détaillée des variables est disponibles sur le site : https://www.data.gouv.fr/organizations/meteo-france/datasets en bas de la page où les données ont été téléchargées sous le titre Documentation dans les fichiers Q_descriptif_champ_*.

Formattage des données

Ajout des variables temporelles et trie des séquences

# GESTION DES DATES
df_temperature_precipitation$AAAAMMJJ <- as.Date(paste0(as.character(df_temperature_precipitation$AAAAMM), "01"), format = "%Y%m%d")
df_temperature_precipitation$year <- format(df_temperature_precipitation$AAAAMMJJ, format = "%Y")
df_temperature_precipitation$month <- format(df_temperature_precipitation$AAAAMMJJ, format = "%m")
df_temperature_precipitation$day <- format(df_temperature_precipitation$AAAAMMJJ, format = "%d")

# AJOUT DE LA VARIABLES ANNEES / MOIS / JOURS
df_snow$AAAAMMJJ <- as.Date(paste0(as.character(df_snow$AAAAMM), "01"), format = "%Y%m%d") 
df_snow$year <- format(df_snow$AAAAMMJJ, format = "%Y") 
df_snow$month <- format(df_snow$AAAAMMJJ, format = "%m")
df_snow$day <- format(df_snow$AAAAMMJJ, format = "%d")

# TRIER LES ANNEES
df_temperature_precipitation <- 
  df_temperature_precipitation |>
  filter(year %in% 1988:2025) 

df_snow <- 
  df_snow |>
  filter(year %in% 1988:2025)

On prend les données à partir de 1988 car on a besoin de décembre 1988 pour mesurer l’épaisseur neigeuse de 1989 (de décembre à mars, voir plus bas)

Filtrer les données par station d’intérêt

On va conserver les données des stations de Tignes et de Val d’Isère, utilisées dans le papier de Marion Tafani (2013). Ces deux stations sont les plus proches du lieu d’étude.

# EXTRACT STATION NAMES
nom_station <- unique(df_temperature_precipitation$NOM_USUEL) 
nom_station_tignes <- nom_station[grepl(pattern = "tignes", ignore.case = TRUE, x = nom_station)]
nom_station_valisere <- nom_station[grepl(pattern = "val.*isere", ignore.case = TRUE, x = nom_station)]

# FILTER
tignes <- df_temperature_precipitation |>
  filter(NOM_USUEL %in% nom_station_tignes) |>
  mutate(location = "tignes", station = NOM_USUEL)

valisere <- df_temperature_precipitation |>
  filter(NOM_USUEL %in% nom_station_valisere) |>
  mutate(location = "valisere", station = NOM_USUEL)

tignes_snow <- df_snow |>
  filter(NOM_USUEL %in% nom_station_tignes) |>
  mutate(location = "tignes", station = NOM_USUEL)

valisere_snow <- df_snow |>
  filter(NOM_USUEL %in% nom_station_valisere) |>
  mutate(location = "valisere", station = NOM_USUEL)

nom_station_tignes

## [1] "TIGNES-BREVIERES" "TIGNES_SAPC"      "Tignes"

nom_station_valisere

## [1] "VAL D ISERE"         "VAL D'ISERE JOS"     "Val d Isere Joseray"

On a 3 stations par site. On va maintenant chercher à déterminer quelles stations sont les plus pertinentes.

CHOIX DE LA STATION METEOROLOGIQUE

Critère 1 : Altitude des stations

##             station  ALTI
##              <char> <int>
## 1: TIGNES-BREVIERES  1570
## 2:      TIGNES_SAPC  2084
## 3:           Tignes  2080

##                station  ALTI
##                 <char> <int>
## 1:         VAL D ISERE  1850
## 2:     VAL D'ISERE JOS  1850
## 3: Val d Isere Joseray  1854

Les stations de Tignes paraissent plus pertinentes car plus proches en altitude de la zone d’étude.

Critère 2 : Données manquantes des séries temporelles par station

Données de précipitation (variables `RR`)

La variable RR est la quantité de précipitations de la journée.

table(tignes$year, !is.na(tignes$RR), tignes$NOM_USUEL) 
table(valisere$year, !is.na(valisere$RR), valisere$NOM_USUEL)

Variable	Station	Couverture inter-annuelle	Couverture intra-annuelle
`RR`	`Tignes`	Continue de 1989 à 2024	~150–220 jours/an (~40–60 %)
`RR`	`TIGNES-BREVIERES`	Complète 1989–2015, absente après	~365 jours/an (quasi complète)
`RR`	`TIGNES_SAPC`	Manquante 1989–1995, complète après	~360–366 jours/an
`RR`	`Val d Isere Joseray`	Continue de 1989 à 2024	~120–170 jours/an (~35–50 %)
`RR`	`VAL D'ISERE JOS`	Manquante 1989–1991	~330–366 jours/an
`RR`	`VAL D ISERE`	Complète jusqu’en 2013, absente après	~250–365 jours/an

On pourrait donc avoir une série complète intra/inter annuelle si on combine les données des stations TIGNES-BREVIERES de 1989 à 1995 et de TIGNES_SAPC de 1996 à 2024.

On vérifie la corrélation entre les deux séries temporelles.

Seules les séries de TIGNES_SAPC et TIGNES-BREVIERES sont bien corrélées. Je complète les données en prenant les données issues de la station avec le plus de données non-NA pour le jour * mois * année considéré.

df_precipitation <- tignes |> 
  filter(!is.na(RR) & NOM_USUEL %in% c("TIGNES_SAPC", "TIGNES-BREVIERES")) |> 
  group_by(year, NOM_USUEL) |> 
  summarise(n_jours = n(), # Compter le nombre de jours disponibles non manquants
    .groups = "drop" # Supprimer le regroupement après le résumé
  ) |>
  group_by(year) |>
  slice_max(
    n_jours,  # Sélectionner la station ayant le plus de jours disponibles
    n = 1, # Une seule station par année
    with_ties = FALSE # En cas d’égalité, ne garder qu’une station
  ) |>
  ungroup() |> # Retirer tout regroupement résiduel
  inner_join( # Rejoindre avec les données journalières complètes
    tignes |> filter(!is.na(RR) & NOM_USUEL %in% c("TIGNES_SAPC", "TIGNES-BREVIERES")), # Table source filtrée sur les précipitations non manquantes
    by = c("year", "NOM_USUEL") # Jointure sur l’année et le nom de la station
  ) |>
  dplyr::select(AAAAMMJJ, year, month, day, NOM_USUEL, RR)

On vérifie qu’on a toutes les données par an par mois :

table(df_precipitation$year, df_precipitation$month, dnn = c("Year", "Month"))

##       Month
## Year   01 02 03 04 05 06 07 08 09 10 11 12
##   1988 31 29 31 30 31 30 31 31 30 31 30 31
##   1989 31 28 31 30 31 30 31 31 30 31 30 31
##   1990 31 28 31 30 31 30 31 31 30 31 30 31
##   1991 31 28 31 30 31 30 31 31 30 31 30 31
##   1992 31 29 31 30 31 30 31 31 30 31 30 31
##   1993 31 28 31 30 31 30 31 31 30 31 30 31
##   1994 31 28 31 30 31 30 31 31 30 31 30 31
##   1995 31 28 31 30 31 30 31 31 30 31 30 31
##   1996 31 29 31 30 31 30 31 31 30 31 30 31
##   1997 31 28 31 30 31 30 31 31 30 31 30 31
##   1998 31 28 31 30 31 30 31 31 30 31 30 31
##   1999 31 28 31 30 31 30 31 31 30 31 30 31
##   2000 31 29 31 30 31 30 31 31 30 31 30 31
##   2001 31 28 31 30 31 30 31 31 30 31 30 31
##   2002 31 28 31 30 31 30 31 31 30 31 30 31
##   2003 31 28 31 30 31 30 31 31 30 31 30 31
##   2004 31 29 31 30 31 30 31 31 30 31 30 31
##   2005 31 27 31 30 31 30 31 31 30 31 30 31
##   2006 31 28 31 30 31 30 31 31 30 31 30 31
##   2007 31 28 31 30 31 30 31 31 30 31 30 31
##   2008 31 29 31 30 31 30 31 31 30 31 30 31
##   2009 31 28 31 30 31 30 31 31 30 31 30 31
##   2010 31 28 31 30 31 30 31 31 30 31 30 31
##   2011 31 28 31 30 31 30 31 31 30 31 30 31
##   2012 31 29 31 30 31 30 31 31 30 31 30 31
##   2013 31 28 31 30 31 30 31 31 30 31 30 31
##   2014 31 28 31 30 31 30 31 31 30 31 30 31
##   2015 31 28 31 30 31 30 31 31 30 31 30 31
##   2016 31 29 31 30 31 30 31 31 30 31 30 31
##   2017 31 28 31 30 31 30 31 31 30 31 30 31
##   2018 31 28 31 16 31 30 31 31 30 31 30 31
##   2019 31 28 31 30 31 30 31 31 30 31 30 31
##   2020 31 29 31 30 31 30 31 31 30 31 30 31
##   2021 31 28 31 30 31 30 31 31 30 31 30 31
##   2022 31 28 31 30 31 30 31 31 30 31 30 31
##   2023 31 28 31 30 31 30 31 31 30 31 30 31
##   2024 31 29 31 30 31 30 31 31 30 31 30 31

Tout bon !

Données de température (variables `TNTXM`)

La variable TNTXM est la moyenne entre la température minimal et la maximale de la journée.

table(tignes$year, !is.na(tignes$TNTXM), tignes$NOM_USUEL)
table(valisere$year, !is.na(valisere$TNTXM), valisere$NOM_USUEL)

Variable	Station	Couverture inter-annuelle	Couverture intra-annuelle
`TNTXM`	`Tignes`	Continue de 1989 à 2024	~150–210 jours/an (~40–60 %)
`TNTXM`	`TIGNES-BREVIERES`	Aucune donnée exploitable	0 jour
`TNTXM`	`TIGNES_SAPC`	Manquante 1989–1991	~360–366 jours/an
`TNTXM`	`Val d Isere Joseray`	Continue de 1989 à 2024	~120–170 jours/an (~35–50 %)
`TNTXM`	`VAL D'ISERE JOS`	Manquante 1989–1991	~350–366 jours/an
`TNTXM`	`VAL D ISERE`	Complète jusqu’en 2013, absente après	~250–365 jours/an

On peut avoir une série complète si on combine les données des stations VAL D ISERE de 1989 à 1991 et VAL D'ISERE JOS de 1992 à 2024.

On vérifie la corrélation des séries temporelles de température entre les stations.

Les trois séries sont très bien corrélées, je complète les séries avec les données de la station avec le maximum de données.

df_temperature <-  valisere |> 
  filter(!is.na(TNTXM)) |> 
  group_by(year, NOM_USUEL) |> 
  summarise(n_jours = n(), # Compter le nombre de jours disponibles non manquants
    .groups = "drop" # Supprimer le regroupement après le résumé
  ) |>
  group_by(year) |>
  slice_max(
    n_jours,  # Sélectionner la station ayant le plus de jours disponibles
    n = 1, # Une seule station par année
    with_ties = FALSE # En cas d’égalité, ne garder qu’une station
  ) |>
  ungroup() |> # Retirer tout regroupement résiduel
  inner_join( # Rejoindre avec les données journalières complètes
    valisere |> filter(!is.na(TNTXM)), # Table source filtrée sur les température non manquantes
    by = c("year", "NOM_USUEL") # Jointure sur l’année et le nom de la station
  ) |>
  dplyr::select(AAAAMMJJ, year, month, day, NOM_USUEL, TNTXM)

On vérifie qu’on a toutes les données par an par mois :

table(df_temperature$year, df_temperature$month, dnn = c("Year", "Month"))

##       Month
## Year   01 02 03 04 05 06 07 08 09 10 11 12
##   1988 31 29 31 30 31 30 31 31 30 31 30 31
##   1989 31 28 31 30 31  0  0  0 30 31 30 31
##   1990 31 28 31 30 31  0 31 31 30 31 30 31
##   1991 31 28 31 30 31 30  0  0 30 31 30 31
##   1992 31 29 31 30 31 30 31  0  0  0  0 31
##   1993 31 28 31 30 31 30 19 31 30 31 30 31
##   1994 31 28 31 30 31 30 31 31 30 31 30 31
##   1995 31 28 31 30 31 28 25 31 30 31 30 31
##   1996 31 29 31 30 31 30 31 31 30 31 30 31
##   1997 31 28 31 30 31 30 31 31 30 31 30 31
##   1998 31 28 31 30 31 30 31 31 30 31 30 31
##   1999 31 28 31 30 31 30 31 31 30 31 30 31
##   2000 31 29 31 30 31 30 31 31 30 31 30 31
##   2001 31 28 31 30 31 30 31 31 30 31 30 31
##   2002 31 28 31 30 31 30 31 31 30 31 30 31
##   2003 31 28 31 30 31 30 31 31 30 31 30 31
##   2004 31 29 31 30 31 30 31 31 30 31  0 31
##   2005 31 28 31 30 31 30 31 31 30 31 30 31
##   2006 31 28 31 30 31 30 31 31 30 31 30 31
##   2007 31 28 31 30 31 30 31 31 30 31 30 31
##   2008 31 29 31 30 31 30 31 31 30 31 30 31
##   2009 31 28 31 30 31 30 31 31 30 31 30 31
##   2010 31 28 31 30 31 30 31 31 30 31 30 31
##   2011 31 28 31 30 31 30 31 31 30 31 30 31
##   2012 31 29 31 30 31 30 31 31 30 31 30 31
##   2013 31 28 31 30 31 30 31 31 30 31 30 31
##   2014 31 28 31 30 31 30 31 31 30 31 30 31
##   2015 31 28 31 30 31 30 31 31 30 31 30 31
##   2016 31 29 31 30 31 30 31 31 30 31 30 31
##   2017 31 28 31 30 31 30 31 31 30 31 30 31
##   2018 21 28 13 19 31 30 31 31 30 31 30 31
##   2019 31 28 31 30 31 30 31 31 30 31 30 31
##   2020 31 29 31 30 31 30 31 31 30 31 30 31
##   2021 31 28 31 30 31 30 31 31 30 31 30 31
##   2022 31 28 31 30 31 30 31 31 30 31 30 31
##   2023 31 28 31 30 31 30 31 31 30 31 30 31
##   2024 31 29 31 30 31 30 31 31 30 31 30 31

Il manque quelques données, notamment fin d’été 1992, mais tout est globalement bon !

Données de neige (variables `NEIGETOT06`)

La variable NEIGETOT06 est l’épaisseur totale de neige au sol.

table(tignes_snow$year, !is.na(tignes_snow$NEIGETOT06), tignes_snow$NOM_USUEL)
table(valisere_snow$year, !is.na(valisere_snow$NEIGETOT06), valisere_snow$NOM_USUEL)

Pour l’épaisseur du manteau neigeux, seules les stations de Tignes et Val d Isere Joseray sont utilisables.

Variable	Station	Couverture inter-annuelle	Couverture intra-annuelle
`NEIGETOT06`	`Tignes`	Continue de 1989 à 2024	~150–220 jours/an
`NEIGETOT06`	`TIGNES-BREVIERES`	Données très discontinues	Variable selon les années
`NEIGETOT06`	`TIGNES_SAPC`	Majoritairement absente avant 2015	Données présentes après 2016
`NEIGETOT06`	`Val d Isere Joseray`	Continue de 1989 à 2024	~120–170 jours/an
`NEIGETOT06`	`VAL D'ISERE JOS`	Quasi totalement absente	Très rares valeurs
`NEIGETOT06`	`VAL D ISERE`	Absente après 2012	Très partielle

Pour l’épaisseur du manteau neigeux, les données sont assez variables en terme de disponiblité. Je regarde la corrélation entre les données pour voir dans quelle mesure on peut compléter.

Seules les données des stations Tigneset TIGNES_SAPC sont assez similaires pour mutuellement se compléter. J’utilise ces deux stations pour prendre la donnée avec le plus d’information par jour * mois * année.

df_snow <-  tignes_snow |> 
  filter(!is.na(NEIGETOT06) & NOM_USUEL %in% c("Tignes", "TIGNES_SAPC")) |> 
  group_by(year, NOM_USUEL) |> 
  summarise(n_jours = n(), # Compter le nombre de jours disponibles non manquants
    .groups = "drop" # Supprimer le regroupement après le résumé
  ) |>
  group_by(year) |>
  slice_max(
    n_jours,  # Sélectionner la station ayant le plus de jours disponibles
    n = 1, # Une seule station par année
    with_ties = FALSE # En cas d’égalité, ne garder qu’une station
  ) |>
  ungroup() |> # Retirer tout regroupement résiduel
  inner_join( # Rejoindre avec les données journalières complètes
    tignes_snow |> filter(!is.na(NEIGETOT06) & NOM_USUEL %in% c("Tignes", "TIGNES_SAPC")), # Table source filtrée sur les température non manquantes
    by = c("year", "NOM_USUEL") # Jointure sur l’année et le nom de la station
  ) |>
  dplyr::select(AAAAMMJJ, year, month, day, NOM_USUEL, NEIGETOT06)

table(df_snow$year, df_snow$month, dnn = c("Year", "Month"))

##       Month
## Year   01 02 03 04 05 06 07 08 09 10 11 12
##   1988 31 29 31 29  0  0  0  0  0  0  0 31
##   1989 31 28 31 30  0  0  0  0  0  0  0 31
##   1990 31 28 31 30  4  0  0  0  0  0  0 31
##   1991 31 28 31 30  7  0  0  0  0  0  0 31
##   1992 31 29 31 30 11  0  0  0  0  0  4 31
##   1993 31 28 31 30  7  0  0  0 17 26 25 31
##   1994 31 28 31 30 14  1 11  6 10 16 15 31
##   1995 31 28 31 30 30  3 23 26 17 25 20 31
##   1996 31 29 31 30 11 11 25 20 20 24 11 31
##   1997 31 28 31 30 14  5 17 20 13 17 14 31
##   1998 31 28 31 30  9  0  0  0 12 20 17 31
##   1999 31 28 31 30 13  4  0  0  0 15 13 31
##   2000 31 29 31 30 13  0  0  0  0  0 15 31
##   2001 31 28 31 30  6  0  0  0  0  0 14 31
##   2002 31 28 31 30  6  0  0  0  0  0 14 31
##   2003 31 28 31 30 22  9  8  0  0  1 25 31
##   2004 31 29 31 30 14  0  0  0  0  0 14 31
##   2005 31 28 31 30  7  0  0  0  0  0  6 28
##   2006 31 28 28 30  6  0  0  0  0  0  1 31
##   2007 31 28 31 30  8  1  1  0  0  0  0 31
##   2008 31 29 31 30  7  0  0  0  0  0  0 31
##   2009 31 28 31 30  7  0  0  0  0  0  1 31
##   2010 31 28 31 30 11  0  0  0  0  0  8 31
##   2011 31 28 31 30  8  0  0  0  0  0  1 31
##   2012 31 29 31 30  8  0  0  0  0  0  2 31
##   2013 31 28 31 30 10  0  0  0  0  0  1 31
##   2014 31 28 31 30  8  0  2  0  0  0  0 31
##   2015 31 28 31 30 10  0  0  0  1  0  0 31
##   2016 31 29 31 30 31 30 31 31 30 31 30 31
##   2017 31 28 31 30 31 30 31 31 30 31 30 31
##   2018  8  0  0 12 31 30 31 31 30 31 30 31
##   2019 31 28 31 30 31 30 31 31 30 31 30 31
##   2020 31 29 31 30 31 30 31 31 30 31 30 31
##   2021 31 28 31 30 31 30 31 31 30 31 30 31
##   2022 31 28 31 30 31 30 31 31 30 31 30 31
##   2023 31 28 31 30 31 30 31 31 30 31 30 31
##   2024 31 29 31 30 31 30 31 31 30 31 30 31

Pas mal de trou, mais la majorité de l’hiver est couvert.

Visualisation des séries temporelles

Précipitations

Température

Snow

NDVI

Les données NDVI sont disponibles à partir de 1981 sur le site du National Oceanic and Atmospheric Administration (NOA) : https://www.ncei.noaa.gov/data/land-normalized-difference-vegetation-index/access/. Ces données correspondent à des enregistrements journaliers sur une grille de 0.05° x 0.05°, soit environ une résolution de 4km x 4km. A noter que les données de 2025 génère une erreur, je pense qu’elles ne sont pas cleans. J’ai donc compute le NDVI jusqu’en 2024 uniquement.

Téléchargement du jeu de données

Il doit être possible d’importer les données directement dans R, mais l’ensemble des données est extrêmement lourd. De plus, le téléchargement est assez long, ce ne serait donc pas si pratique. J’ai donc décidé de les télécharger sur mon ordinateur en local directement depuis R, puis ensuite de les charger. J’utilise le package future pour paralléliser le téléchargement et le rendre plus rapide.

# SET UP PARALLELISATION
future::plan(future::multisession, workers = parallel::detectCores() - 1)  # utilise tous les coeurs dispo -1 pour ne pas bloquer l'ordinateur

# FUNCTION TO EXTRACT YEAR, MONTH AND DAY
extract_date <- function(file_names){
  date_raw <- stringr::str_extract(string = file_names, pattern = "_\\d{8}_")
  date_raw <- stringr::str_remove_all(date_raw, "_")
  date_raw <- as.Date(paste0(date_raw), format = "%Y%m%d")  # convertit en date
  year_temp <- (format(date_raw, "%Y"))  # récupère le mois
  month_temp <- (format(date_raw, "%m"))  # récupère le mois
  day_temp <- (format(date_raw, "%d"))    # récupère le jour
  
  return(list(date = date_raw, 
              year = year_temp,
              month = month_temp,
              day = day_temp))
}

Le code suivant a l’argument eval = F, il est à lancer “à la main” pour télécharger les fichiers. Penser à changer le chemin d’accès dest_dir.

# DOWNLOAD FILES FROM NOAA
for(year in 1989:2025){  # boucle sur chaque année du dataset

  # Spécifier l'URL où les fichiers sont stockés et extraire les noms de fichiers
  url <- paste0("https://www.ncei.noaa.gov/data/land-normalized-difference-vegetation-index/access/", year, "/")
  files <- readLines(url)  # lit la page web ligne par ligne
  file_names <- stringr::str_extract(files, "href=\"([^\"]+)\"")  # extrait les liens href
  file_names <- stringr::str_replace_all(file_names, "href=|\"", "")  # supprime 'href=' et les guillemets
  file_names <- file_names[file_names != ".."]  # enlève les références parent (..) dans la liste
  file_names <- as.character(na.omit(file_names))  # supprime les NA et convertit en vecteur caractère
  
  # Extraire les fichiers NDVI 
  month_temp <- extract_date(file_names)$month  # récupère le mois
  day_temp <- extract_date(file_names)$day    # récupère le jour
  file_names_filtered <- file_names[month_temp %in% c("03", "04", "05", "06", "07", "08", "09") & day_temp %in% c("01", "05", "10", "15", "20", "25")]  # Prendre quelques jours par mois

    # Créer le répertoire de destination si inexistant
  dest_dir <- paste0("C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/ENVIRONMENTAL_VARIABLES/NDVI/NOAA/", year)  # chemin du dossier de l'année
  
  # Checker si le fichier est déjà présent 
  existing_files <- list.files(dest_dir)
  
  # télécharger les fichiers non-présents uniquement 
  files_to_download <- setdiff(file_names_filtered, existing_files)
  
  dir.create(dest_dir, recursive = TRUE, showWarnings = FALSE)  # crée le dossier, ignore les avertissements si déjà présent

  # Fonction pour télécharger un fichier
  download_one <- function(f, base_url, dest_dir) {
    httr::GET(
      url = paste0(base_url, f),                  # URL complète du fichier
      httr::write_disk(file.path(dest_dir, f), overwrite = TRUE)  # chemin où sauvegarder le fichier
      # httr::progress()                            # montre la progression du téléchargement
    )
  }

 # Télécharger uniquement si nécessaire
if (length(files_to_download) > 0) {
  furrr::future_walk(
    files_to_download,
    download_one,
    base_url = url,
    dest_dir = dest_dir
  )
  message(year, " : ✔")
} else {
  message(year, " : ✔")
}
}

Computation du NDVI

On commence par tracer les frontière de la zone d’étude (fait avec Google Maps).

# BOUNDARIES POUR LA SASSIERE
longitude <- c(6.980740, 6.980891, 6.996077, 6.995761)  # NW, SW, SE, NE
latitude <- c(45.492250, 45.487339, 45.487896, 45.490588)  # NW, SW, SE, NE

lon <- c(min(longitude), max(longitude))
lat <- c(min(latitude), max(latitude))

On visualise les coordonnées extraites :

En cas de réutilisation du code, penser à changer le chemin d’accès de l’objet nc_dir vers le dossier où ont été téléchargé les données NDVI avec le script ci-dessus.

# CONVERT SPATIALLY THE BOUNDARIES
ext_ndvi <- terra::ext(min(lon), max(lon), min(lat), max(lat))


df_ndvi <- data.frame(file = as.character(), year = as.numeric(), mean_ndvi = as.numeric())

# COMPUTE NDVI
for(year in 1989:2024) {
  # boucle sur chaque année du dataset
  nc_dir <- paste0("C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/ENVIRONMENTAL_VARIABLES/NDVI/NOAA/", year)  # A MODIFIER SI FICHIERS NDVI PAS DANS DOCUMENT
  nc_files <- list.files(nc_dir, pattern = "\\.nc$", full.names = TRUE)
  ndvi_mean <- data.frame(file = basename(nc_files), mean_ndvi = NA_real_)
  
  mean_ndvi_from_nc <- function(nc_file, extent_obj) {
    r <- terra::rast(nc_file)  # create the raster of the zone
    # crs(r) <- "EPSG:4326"  # doit être défini sinon bug en 2025 ("crop do not overlap")
    r <- terra::crop(r, extent_obj)  # crop the zone to exxtract the specific NDVI
    
    # extract ndvi from layer
      terra::global(r, mean, na.rm = TRUE)[1, 1]  
  }
  
  # compute NDVI for each file of year t
  ndvi_mean["mean_ndvi"] <- sapply(nc_files, mean_ndvi_from_nc, extent_obj = ext_ndvi)
  df_ndvi <- rbind(df_ndvi, ndvi_mean)
}


# ADD YEAR, MONTH AND DAY
# Extraire les fichiers NDVI du 1er avril au 1er septembre
df_ndvi$date <- extract_date(df_ndvi$file)$date  # récupère la date
df_ndvi$year <- extract_date(df_ndvi$file)$year  # récupère l'année
df_ndvi$month <- extract_date(df_ndvi$file)$month  # récupère le mois
df_ndvi$day <- extract_date(df_ndvi$file)$day    # récupère le jour

# SAUVEGARDE DES DONNEES
save(df_ndvi, file = "C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/ENVIRONMENTAL_VARIABLES/NDVI/NOAA/df_ndvi.rds")

# LOAD
load("C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/ENVIRONMENTAL_VARIABLES/NDVI/NOAA/df_ndvi.rds")

NDVI issu de GOOGLE ENGINE

J’ai aussi compute les valeurs NDVI via la collection LANDSAT de Google Earth Engine qui me permettait d’avoir moins de points, mais une résolution spatiale beaucoup plus fine (~ 30m) pour comparer. Pour LANDSAT, j’ai également “réduit” la zone à la plaine de Milloz car la résolution le permettait, pour éviter d’inclure un pic montagneux dans la zone (et donc de la neige).

Le code est disponible dans le script rgee.R dont le code est disponible en fin de ce document. Les données générées ont été enregistrées en .csv car le package rgee permettant de se connecter à Google engine depuis R est assez compliqué à intégrer dans Rmarkdown.

# IMPORT DATA
df_gee <- read.csv("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/ENVIRONMENTAL_VARIABLES/NDVI/LANDSTAT/NDVI_LANDSAT_MILLOZ.csv", header = T)

# ADD DATE
df_gee <- df_gee |> select(-.geo)
df_gee$day <- (format(as.Date(df_gee$date), format = "%d"))    # récupère le jour
df_gee$month <- (format(as.Date(df_gee$date), format = "%m"))  # récupère le mois
df_gee$year <- (format(as.Date(df_gee$date), format = "%Y"))  # récupère le mois

Comparaison NDVI NOAA (~5km\(^2\)) et LANDSAT (~30m\(^2\))

Visualisation des séries temporelles NDVI

NAO

Le Northern Atlantic Oscillation index (NAO) mesure une différence de pression atmosphérique entre deux régions du globe, l’anticyclone des Açores et la dépression d’Islande.

Indice NAO	Europe du Nord / Ouest	Europe du Sud / Méditerranée	Impacts typiques
NAO positif (haut)	Hivers doux et humides, tempêtes fréquentes	Temps plus sec, risque de sécheresse	Inondations au nord, vents forts, douceur hivernale
NAO négatif (bas)	Hivers froids, neige fréquente, conditions plus sèches	Précipitations accrues, épisodes pluvieux intenses	Vagues de froid, inondations méditerranéennes
NAO positif persistant	Températures supérieures à la normale, pluies fréquentes	Déficit pluviométrique prolongé	Stress hydrique au sud, sols saturés au nord
NAO négatif persistant	Froid prolongé, épisodes neigeux durables	Cyclogenèse et pluies intenses	Gel prolongé, crues rapides

Les données de NAO peut-être téléchargé sur le site https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml et est disponible de 1950 à 2026. Comme c’est un site américain, j’ai téléchargé le fichier .csv par sécurité. Je mets néanmoins le code pour l’importer directement depuis R.

Téléchargement des données depuis NOAA website

# URL OF THE DATA
url_nao <- "https://ftp.cpc.ncep.noaa.gov/cwlinks/norm.daily.nao.cdas.z500.19500101_current.csv"

# READ .CSV AND FORMATING
df_nao <- read.csv(file = url_nao, header = T)
colnames(df_nao)[4] <- "nao"
df_nao <-
  df_nao |>
  filter(year %in% 1988:2026) |>
  mutate(year = as.character(year),
         month = as.character(month),
         day = as.character(day),
         month = factor(month,
                        levels = as.character(1:12),
                        labels = as.character(1:12)))

##       MONTH
## YEAR    1  2  3  4  5  6  7  8  9 10 11 12
##   1988 31 29 31 30 31 30 31 31 30 31 30 31
##   1989 31 28 31 30 31 30 31 31 30 31 30 31
##   1990 31 28 31 30 31 30 31 31 30 31 30 31
##   1991 31 28 31 30 31 30 31 31 30 31 30 31
##   1992 31 29 31 30 31 30 31 31 30 31 30 31
##   1993 31 28 31 30 31 30 31 31 30 31 30 31
##   1994 31 28 31 30 31 30 31 31 30 31 30 31
##   1995 31 28 31 30 31 30 31 31 30 31 30 31
##   1996 31 29 31 30 31 30 31 31 30 31 30 31
##   1997 31 28 31 30 31 30 31 31 30 31 30 31
##   1998 31 28 31 30 31 30 31 31 30 31 30 31
##   1999 31 28 31 30 31 30 31 31 30 31 30 31
##   2000 31 29 31 30 31 30 31 31 30 31 30 31
##   2001 31 28 31 30 31 30 31 31 30 31 30 31
##   2002 31 28 31 30 31 30 31 31 30 31 30 31
##   2003 31 28 31 30 31 30 31 31 30 31 30 31
##   2004 31 29 31 30 31 30 31 31 30 31 30 31
##   2005 31 28 31 30 31 30 31 31 30 31 30 31
##   2006 31 28 31 30 31 30 31 31 30 31 30 31
##   2007 31 28 31 30 31 30 31 31 30 31 30 31
##   2008 31 29 31 30 31 30 31 31 30 31 30 31
##   2009 31 28 31 30 31 30 31 31 30 31 30 31
##   2010 31 28 31 30 31 30 31 31 30 31 30 31
##   2011 31 28 31 30 31 30 31 31 30 31 30 31
##   2012 31 29 31 30 31 30 31 31 30 31 30 31
##   2013 31 28 31 30 31 30 31 31 30 31 30 31
##   2014 31 28 31 30 31 30 31 31 30 31 30 31
##   2015 31 28 31 30 31 30 31 31 30 31 30 31
##   2016 31 29 31 30 31 30 31 31 30 31 30 31
##   2017 31 28 31 30 31 30 31 31 30 31 30 31
##   2018 31 28 31 30 31 30 31 31 30 31 30 31
##   2019 31 28 31 30 31 30 31 31 30 31 30 31
##   2020 31 29 31 30 31 30 31 31 30 31 30 31
##   2021 31 28 31 30 31 30 31 31 30 31 30 31
##   2022 31 28 31 30 31 30 31 31 30 31 30 31
##   2023 31 28 31 30 31 30 31 31 30 31 30 31
##   2024 31 29 31 30 31 30 31 31 30 31 30 31
##   2025 31 28 31 30 31 30 31 31 30 31 30 31
##   2026 31 16  0  0  0  0  0  0  0  0  0  0

Visualisation de la série temporelle NAO

INDICE PAR SAISONS

On va caractériser les conditions climatiques des 3 périodes majeures des marmottes Alpines :

été à t-1 -> mi juin-août : accumulation des réserves pour l’hibernation [summer_condition]
hiver à t-1 -> décémbre-mars : condition d’hibernation [winter_condition]
printemps -> 15 avril-1er mai : condition à l’émergence et pour la reproduction [spring_condition]

On va opérationaliser les variables de la façon suivante :

summer_condition : (1) Précipitation, (2) date du NDVI maximum et (3) indice d’aridité Baynoul-Gaussen compute as \(RR -2 \cdot TNTXM\) et (4) NAO de l’été
winter_condition : (1) Température moyenne, (2) épaisseur neigeuse et (3) NAO de l’hiver
spring_condition : (1) Température moyenne, (2) NDVI, (3) NAO du printemps

Summer condition

# MEAN TEMPERATURE
summer_temperature <- 
  df_temperature |> 
  group_by(year) |> 
  filter(month %in% c("06", "07", "08")) |> 
  summarise(mean_temperature = mean(TNTXM),
            se = sd(TNTXM) / sqrt(n()),
            precision = 1/var(TNTXM, na.rm = T),
            ci_lower = mean_temperature - se,
            ci_upper = mean_temperature + se)

# MEAN PRECIPITATION
summer_precipitation <- 
  df_precipitation |> 
  group_by(year) |> 
  filter(month %in% c("06", "07", "08")) |> 
  summarise(rr = mean(RR),
            se = sd(RR) / sqrt(n()),
            precision = 1/var(RR, na.rm = T),
            ci_lower = rr - se,
            ci_upper = rr + se)


# DATE OF MAXIMUM NDVI
# summer_max_ndvi_date <- 
#   df_ndvi |> 
#     group_by(year) |> 
#     filter((month %in% "06" & day %in% c(15:31)) | (month %in% c("07", "08"))) |> 
#   filter(mean_ndvi == max(mean_ndvi, na.rm = T)) |> 
#   select(mean_ndvi, date, year, month, day) |> 
#   as.data.frame()
# 
# #' convert date to julian day
# summer_max_ndvi_date$julian_day_max_ndvi <- as.POSIXlt(as.Date(summer_max_ndvi_date$date))$yday + 1
# 
# summer_max_ndvi_date <- 
#   summer_max_ndvi_date |> 
#   mutate(se = NA,
#     ci_lower = NA, 
#          ci_upper = NA, 
#          precision = NA) |> 
#   select(year, julian_day_max_ndvi, se, precision, ci_lower, ci_upper) |> 
#   as_tibble()


# SUMMER NDVI 
summer_ndvi <- 
  df_ndvi |> 
  group_by(year) |> 
  filter(month %in% c("06", "07", "08")) |> 
    reframe(ndvi = mean(mean_ndvi, na.rm = T),
            precision = 1/var(mean_ndvi, na.rm = T),
            se = sd(mean_ndvi, na.rm = T) / sqrt(n()),
            ci_lower = mean(mean_ndvi, na.rm = T) - se,
            ci_upper = mean(mean_ndvi, na.rm = T) + se)


# SUMMER NAO
summer_nao <-
  df_nao |>
  group_by(year) |>
  filter(month %in% c("6", "7", "8")) |> 
  summarise(mean_nao = mean(nao),
            precision = 1/var(nao, na.rm = T),
            se = sd(nao) / sqrt(n()),
            ci_lower = mean(nao) - se,
            ci_upper = mean(nao) + se) |>
  mutate(year = as.character(year))

# BGI
df_merge <- merge(df_temperature, df_precipitation, by = c("year", "month", "day"))

summer_bgi <- 
  df_merge |> 
  group_by(year) |> 
  filter(month %in% c("06", "07", "08")) |> 
  reframe(bgi = mean(RR - 2*TNTXM, na.rm = T),
            precision = 1/var(RR - 2*TNTXM, na.rm = T),
            se = sd(RR - 2*TNTXM) / sqrt(n()),
            ci_lower = mean(RR - 2*TNTXM) - se,
            ci_upper = mean(RR - 2*TNTXM) + se) |> 
  as.data.frame()

Winter condition

# TEMPERATURE
winter_temperature <- 
  df_temperature |>
  filter(month %in% c("12", "01", "02", "03")) |>
  mutate(
    winter_year = if_else(month == "12",
                          as.numeric(as.character(year))+1, 
                          as.numeric(as.character(year)))
  ) |>
  group_by(winter_year) |>
  reframe(
    mean_temperature = mean(TNTXM, na.rm = TRUE),
    precision = 1 / var(TNTXM, na.rm = TRUE),
    se = sd(TNTXM, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean_temperature - se,
    ci_upper = mean_temperature + se
  ) |> 
  rename(year = winter_year) |> 
  mutate(year = as.character(year))


# WINTER PRECIPITATION
winter_precipitation <- 
  df_precipitation |>
  filter(month %in% c("12", "01", "02", "03")) |>
  mutate(
    winter_year = if_else(month == "12",
                          as.numeric(as.character(year))+1, 
                          as.numeric(as.character(year)))
  ) |>
  group_by(winter_year) |>
  reframe(
    mean_precipitation = mean(RR, na.rm = TRUE),
    precision = 1 / var(RR, na.rm = TRUE),
    se = sd(RR, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean_precipitation - se,
    ci_upper = mean_precipitation + se
  ) |> 
  rename(year = winter_year) |> 
  mutate(year = as.character(year))


# SNOW COVER DURING WINTER
winter_snow <- 
  df_snow |>
  filter(month %in% c("12", "01", "02", "03")) |>
  mutate(
    winter_year = if_else(month == "12",
                          as.numeric(as.character(year))+1, 
                          as.numeric(as.character(year)))
  ) |>
  group_by(winter_year) |>
  reframe(
    mean_snow = mean(NEIGETOT06, na.rm = TRUE),
    precision = 1 / var(NEIGETOT06, na.rm = TRUE),
    se = sd(NEIGETOT06, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean_snow - se,
    ci_upper = mean_snow + se
  ) |> 
  rename(year = winter_year) |> 
  mutate(year = as.character(year))
  

# WINTER NAO
winter_nao <- 
  df_nao |>
  filter(month %in% c("12", "1", "2", "3")) |>
  mutate(
    winter_year = if_else(month == "12",
                          as.numeric(as.character(year))+1, 
                          as.numeric(as.character(year)))
  ) |>
  group_by(winter_year) |>
  reframe(
    mean_nao = mean(nao, na.rm = TRUE),
    precision = 1 / var(nao, na.rm = TRUE),
    se = sd(nao, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean_nao - se,
    ci_upper = mean_nao + se
  ) |> 
  rename(year = winter_year) |> 
  mutate(year = as.character(year))

Spring condition

# TEMPERATURE
spring_temperature <- 
  df_temperature |> 
  group_by(year) |> 
    filter(month %in% c("04", "05")) |> 
  reframe(mean_temperature = mean(TNTXM),
            precision = 1/var(TNTXM, na.rm = T),
            se = sd(TNTXM) / sqrt(n()),
            ci_lower = mean_temperature - se,
            ci_upper = mean_temperature + se)


# PRECIPITATION
spring_precipitation <- 
  df_precipitation |> 
  group_by(year) |> 
    filter(month %in% c("04", "05")) |> 
  reframe(mean_precipitation = mean(RR),
            precision = 1/var(RR, na.rm = T),
            se = sd(RR) / sqrt(n()),
            ci_lower = mean_precipitation - se,
            ci_upper = mean_precipitation + se)


# BGI
# RAPPEL : 
# ... df_merge <- merge(df_temperature,
# ... df_precipitation,
# ... by = c("year", "month", "day"))

spring_bgi <-
  df_merge |>
  group_by(year) |>
   filter(month %in% c("04", "05")) |>
  reframe(bgi = mean(RR - 2*TNTXM, na.rm = T),
            precision = 1/var(RR - 2*TNTXM, na.rm = T),
            se = sd(RR - 2*TNTXM) / sqrt(n()),
            ci_lower = mean(RR - 2*TNTXM) - se,
            ci_upper = mean(RR - 2*TNTXM) + se) |>
  as.data.frame()


# SPRING NDVI
spring_ndvi <- 
  df_ndvi |> 
  group_by(year) |> 
    filter(month %in% c("04", "05")) |> 
    reframe(ndvi = mean(mean_ndvi, na.rm = T),
            precision = 1/var(mean_ndvi, na.rm = T),
            se = sd(mean_ndvi, na.rm = T) / sqrt(n()),
            ci_lower = mean(mean_ndvi, na.rm = T) - se,
            ci_upper = mean(mean_ndvi, na.rm = T) + se)

# SPRING NAO
spring_nao <-
  df_nao |>
  group_by(year) |>
    filter(month %in% c(4, 5)) |>
    summarise(mean_nao = mean(nao),
              precision = 1/var(nao, na.rm = T),
              se = sd(nao) / sqrt(n()),
            ci_lower = mean_nao - se,
            ci_upper = mean_nao + se) |>
  mutate(year = as.character(year))

GENERATE DATAFRAME

Long format

# Fonction générique pour standardiser un data.frame
standardize_df <- function(df,
                           value_col,
                           variable_name,
                           season = NA) {
  df |> 
    mutate(year = as.integer(year)) |> 
    transmute(  # transmute ne conserve que les colonnes spécifiées
      year,
      variable = variable_name,
      season = season,
      value = .data[[value_col]],
      ci_lower,
      ci_upper,
      precision
    )
}

# Liste des data.frames standardisés
df_abiotic_long <- bind_rows(
  # Été
  standardize_df(summer_precipitation, "rr", "precipitation", "summer"),
  standardize_df(summer_nao, "mean_nao", "nao", "summer"),
  standardize_df(summer_temperature, "mean_temperature", "temperature", "summer"),
  standardize_df(summer_bgi, "bgi", "bgi", "summer"),
  standardize_df(summer_ndvi, "ndvi", "ndvi", "summer"),
  # standardize_df(summer_max_ndvi_date, "julian_day_max_ndvi", "julian_day_max_ndvi", "summer"),

    # Hiver
  standardize_df(winter_temperature, "mean_temperature", "temperature", "winter"),
  standardize_df(winter_precipitation, "mean_precipitation", "precipitation", "winter"),
  standardize_df(winter_snow, "mean_snow", "snow", "winter"),
  standardize_df(winter_nao, "mean_nao", "nao", "winter"),
  
  # Printemps
  standardize_df(spring_temperature, "mean_temperature", "temperature", "spring"),
  standardize_df(spring_precipitation, "mean_precipitation", "precipitation", "spring"),
  standardize_df(spring_nao, "mean_nao", "nao", "spring"),
  standardize_df(spring_ndvi, "ndvi", "ndvi", "spring"),
  standardize_df(spring_bgi, "bgi", "bgi", "spring")
)

# aperçu
glimpse(df_abiotic_long)

## Rows: 521
## Columns: 7
## $ year      <int> 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, …
## $ variable  <chr> "precipitation", "precipitation", "precipitation", "precipitation", "precipitation", "precipitation", "precipitation", "precipitation", "precipitation", "precipitation", "precipita…
## $ season    <chr> "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", "summer", …
## $ value     <dbl> 2.15760870, 2.00434783, 2.65869565, 2.40108696, 2.90652174, 2.47391304, 3.09456522, 1.67065217, 4.26413043, 3.26413043, 2.80543478, 3.30869565, 2.49021739, 3.64456522, 3.83913043, …
## $ ci_lower  <dbl> 1.69831429, 1.56951171, 2.10683620, 1.91948196, 2.27674462, 1.94737673, 2.38211139, 1.31697199, 3.10817452, 2.62719855, 2.21751619, 2.69333515, 1.94370419, 2.80542851, 3.14701741, …
## $ ci_upper  <dbl> 2.6169030965, 2.4391839432, 3.2105551076, 2.8826919500, 3.5362988563, 3.0004493585, 3.8070190434, 2.0243323617, 5.4200863464, 3.9010623229, 3.3933533705, 3.9240561507, 3.0367305967…
## $ precision <dbl> 0.051526408, 0.057485850, 0.035690711, 0.046863013, 0.027405544, 0.039206269, 0.021414046, 0.086894184, 0.008134477, 0.026793299, 0.031446896, 0.028704695, 0.036392413, 0.015436425…

# save
save(df_abiotic_long, file = "C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/df_abiotic_long.rds")

Short format

df_abiotic_short <- df_abiotic_long %>%
  pivot_longer(
    cols = c(value, ci_lower, ci_upper, precision),
    names_to = "stat",
    values_to = "val"
  ) %>%
  mutate(
    stat = dplyr::recode(stat,
                  value = "y",
                  ci_lower = "ci_lower",
                  ci_upper = "ci_upper",
                  precision = "precision")
  ) %>%
  pivot_wider(
    id_cols = year,
    names_from = c(season, variable, stat),
    values_from = val,
    names_sep = "_"
  )

# Aperçu
glimpse(df_abiotic_short)

## Rows: 39
## Columns: 57
## $ year                           <int> 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 201…
## $ summer_precipitation_y         <dbl> 2.157609, 2.004348, 2.658696, 2.401087, 2.906522, 2.473913, 3.094565, 1.670652, 4.264130, 3.264130, 2.805435, 3.308696, 2.490217, 3.644565, 3.839130, 2.272826,…
## $ summer_precipitation_ci_lower  <dbl> 1.698314, 1.569512, 2.106836, 1.919482, 2.276745, 1.947377, 2.382111, 1.316972, 3.108175, 2.627199, 2.217516, 2.693335, 1.943704, 2.805429, 3.147017, 1.840850,…
## $ summer_precipitation_ci_upper  <dbl> 2.616903, 2.439184, 3.210555, 2.882692, 3.536299, 3.000449, 3.807019, 2.024332, 5.420086, 3.901062, 3.393353, 3.924056, 3.036731, 4.483702, 4.531243, 2.704802,…
## $ summer_precipitation_precision <dbl> 0.051526408, 0.057485850, 0.035690711, 0.046863013, 0.027405544, 0.039206269, 0.021414046, 0.086894184, 0.008134477, 0.026793299, 0.031446896, 0.028704695, 0.0…
## $ summer_nao_y                   <dbl> 0.12547826, 0.09914130, 0.21504348, 0.07701087, 0.17538044, -0.54158696, 0.49968478, 0.10159783, 0.33613044, -0.01680435, -0.53920652, 0.12592391, -0.24552174,…
## $ summer_nao_ci_lower            <dbl> 0.04636496, 0.02487485, 0.15490086, -0.01713398, 0.09804032, -0.63427103, 0.41750672, 0.03178352, 0.27436436, -0.10883872, -0.62746861, 0.04315418, -0.32593093…
## $ summer_nao_ci_upper            <dbl> 0.2045915660, 0.1734077588, 0.2751860935, 0.1711557119, 0.2527205464, -0.4489028915, 0.5818628422, 0.1714121369, 0.3978965075, 0.0752300252, -0.4509444310, 0.2…
## $ summer_nao_precision           <dbl> 1.736653, 1.970728, 3.005021, 1.226363, 1.817199, 1.265325, 1.609535, 2.230094, 2.849130, 1.283252, 1.395288, 1.586606, 1.681128, 2.731721, 2.973311, 1.344841,…
## $ summer_temperature_y           <dbl> 11.054348, NA, 12.367742, 8.816667, 10.336066, 11.347500, 12.617391, 11.122619, 10.473913, 10.773913, 11.846739, 11.207609, 11.275000, 11.235870, 11.468478, 14…
## $ summer_temperature_ci_lower    <dbl> 10.724526, NA, 12.071781, 8.176335, 9.889579, 11.018453, 12.262760, 10.676670, 10.175350, 10.456171, 11.503761, 10.944581, 10.954843, 10.894165, 11.184948, 14.…
## $ summer_temperature_ci_upper    <dbl> 11.384170, NA, 12.663703, 9.456999, 10.782552, 11.676547, 12.972023, 11.568568, 10.772476, 11.091655, 12.189717, 11.470636, 11.595157, 11.577574, 11.752008, 14…
## $ summer_temperature_precision   <dbl> 0.09992004, NA, 0.18413678, 0.08129587, 0.08223451, 0.11544990, 0.08642860, 0.05986180, 0.12193857, 0.10766206, 0.09240183, 0.15711199, 0.10604423, 0.09309162,…
## $ summer_bgi_y                   <dbl> -19.95109, NA, -23.23065, -13.55000, -17.97705, -20.71000, -22.14022, -20.65476, -16.68370, -18.28370, -20.88804, -19.10652, -20.05978, -18.82717, -19.09783, -…
## $ summer_bgi_ci_lower            <dbl> -20.74155, NA, -24.00118, -15.32341, -19.16798, -21.48195, -23.19536, -21.64999, -17.88865, -19.24119, -21.77784, -19.92144, -20.95087, -20.01183, -19.99908, -…
## $ summer_bgi_ci_upper            <dbl> -19.16062, NA, -22.46011, -11.77659, -16.78612, -19.93805, -21.08508, -19.65953, -15.47874, -17.32620, -19.99825, -18.29160, -19.16869, -17.64251, -18.19657, -…
## $ summer_bgi_precision           <dbl> 0.017395793, NA, 0.027166025, 0.010598892, 0.011558405, 0.020976317, 0.009763193, 0.012019217, 0.007486357, 0.011856020, 0.013728851, 0.016367448, 0.013688943,…
## $ summer_ndvi_y                  <dbl> NA, 0.14790588, 0.19090556, 0.15554444, 0.09657778, 0.15134706, 0.12044118, 0.07385000, 0.16038000, 0.13297647, 0.16758889, 0.13557647, 0.11644706, 0.13225000,…
## $ summer_ndvi_ci_lower           <dbl> NA, 0.11958334, 0.16137764, 0.12408764, 0.08155751, 0.11166693, 0.09215471, 0.05669550, 0.13084197, 0.10019575, 0.13992383, 0.11322230, 0.09018792, 0.09788998,…
## $ summer_ndvi_ci_upper           <dbl> NA, 0.1762284, 0.2204335, 0.1870013, 0.1115980, 0.1910272, 0.1487276, 0.0910045, 0.1899180, 0.1657572, 0.1952539, 0.1579306, 0.1427062, 0.1666100, 0.1489616, 0…
## $ summer_ndvi_precision          <dbl> NA, 69.25690, 63.71795, 56.14334, 246.24758, 35.28429, 69.43369, 188.78660, 63.67434, 51.69998, 72.58791, 111.17582, 80.56872, 47.05660, 35.68627, 30.40610, 77…
## $ winter_temperature_y           <dbl> -4.521978, -1.892562, -1.466942, -4.114876, -3.287705, -3.008264, -2.940496, -4.013223, -3.459836, -1.740496, -2.670248, -4.152066, -4.202459, -1.876033, -2.93…
## $ winter_temperature_ci_lower    <dbl> -4.921241, -2.218477, -1.779834, -4.568607, -3.590500, -3.372245, -3.335636, -4.386654, -3.827188, -2.075733, -3.035694, -4.570139, -4.550624, -2.292408, -3.30…
## $ winter_temperature_ci_upper    <dbl> -4.1227150, -1.5666469, -1.1540502, -3.6611455, -2.9849099, -2.6442843, -2.5453557, -3.6397926, -3.0924842, -1.4052591, -2.3048024, -3.7339928, -3.8542944, -1.…
## $ winter_temperature_precision   <dbl> 0.06893509, 0.07780470, 0.08441622, 0.04014381, 0.08940106, 0.06238200, 0.05293125, 0.05926456, 0.06074013, 0.07353790, 0.06188269, 0.04728352, 0.06761933, 0.0…
## $ winter_precipitation_y         <dbl> 5.824176, 2.311570, 4.294215, 2.638843, 2.427049, 2.004132, 4.139669, 6.087603, 2.831967, 2.030579, 3.393388, 3.242975, 3.736885, 4.141322, 2.604132, 2.601653,…
## $ winter_precipitation_ci_lower  <dbl> 4.9725949, 1.7249398, 3.3055154, 2.0660912, 1.6671031, 1.4741699, 3.4352359, 5.2279734, 2.2054885, 1.5051208, 2.6116797, 2.6251045, 3.0044814, 3.3627909, 1.960…
## $ winter_precipitation_ci_upper  <dbl> 6.675757, 2.898201, 5.282914, 3.211595, 3.186995, 2.534095, 4.844103, 6.947233, 3.458446, 2.556036, 4.175097, 3.860846, 4.469289, 4.919854, 3.247955, 3.257258,…
## $ winter_precipitation_precision <dbl> 0.015153284, 0.024015159, 0.008454464, 0.025193108, 0.014193014, 0.029425555, 0.016654616, 0.011183857, 0.020884663, 0.029932227, 0.013524609, 0.021648081, 0.0…
## $ winter_snow_y                  <dbl> 134.35165, 82.84298, 71.95868, 111.53719, 99.67213, 131.81818, 123.54545, 146.38017, 84.72131, 121.00000, 91.05785, 105.36364, 125.39344, 71.40496, 72.76033, 8…
## $ winter_snow_ci_lower           <dbl> 128.92042, 80.14424, 65.69492, 108.38971, 96.19960, 130.25025, 120.11622, 138.82673, 81.27204, 118.94524, 88.22212, 100.45550, 121.43584, 68.92911, 69.30188, 8…
## $ winter_snow_ci_upper           <dbl> 139.78288, 85.54172, 78.22243, 114.68467, 103.14466, 133.38611, 126.97469, 153.93360, 88.17058, 123.05476, 93.89358, 110.27177, 129.35104, 73.88081, 76.21878, …
## $ winter_snow_precision          <dbl> 0.0003725308, 0.0011347301, 0.0002106421, 0.0008342366, 0.0006797490, 0.0033617212, 0.0007027818, 0.0001448522, 0.0006889463, 0.0019574579, 0.0010277440, 0.000…
## $ winter_nao_y                   <dbl> 0.39061538, 0.76978512, 0.42530578, 0.36851240, 0.43532787, 0.51536364, 0.65049587, 0.78293388, -0.22462295, 0.23099174, 0.14077686, 0.33494215, 0.66205737, -0…
## $ winter_nao_ci_lower            <dbl> 0.349724118, 0.715275276, 0.348456976, 0.307558422, 0.371987152, 0.467455065, 0.595475460, 0.734670162, -0.292627604, 0.137878065, 0.072290083, 0.293835246, 0.…
## $ winter_nao_ci_upper            <dbl> 0.43150665, 0.82429497, 0.50215459, 0.42946637, 0.49866858, 0.56327222, 0.70551627, 0.83119760, -0.15661830, 0.32410541, 0.20926364, 0.37604905, 0.73101343, -0…
## $ winter_nao_precision           <dbl> 6.5719989, 2.7814082, 1.3993951, 2.2243880, 2.0430276, 3.6007098, 2.7300276, 3.5479141, 1.7724047, 0.9532081, 1.7619808, 4.8908588, 1.7238338, 2.1798293, 1.733…
## $ spring_temperature_y           <dbl> 4.413115, 3.791803, 3.731148, 1.386885, 4.883607, 4.413115, 3.331148, 3.144262, 4.403279, 4.408197, 4.219672, 4.603279, 5.031148, 4.260656, 4.070492, 5.259016,…
## $ spring_temperature_ci_lower    <dbl> 4.0090002, 3.1950072, 3.2059333, 0.8510111, 4.2675678, 3.8912409, 2.6624462, 2.7495048, 3.8990996, 3.8804979, 3.6127457, 3.9827027, 4.5463371, 3.6140747, 3.646…
## $ spring_temperature_ci_upper    <dbl> 4.817229, 4.388599, 4.256362, 1.922759, 5.499645, 4.934989, 3.999849, 3.539020, 4.907458, 4.935896, 4.826599, 5.223855, 5.515958, 4.907237, 4.494771, 5.887743,…
## $ spring_temperature_precision   <dbl> 0.10038321, 0.04602759, 0.05942884, 0.05708797, 0.04319705, 0.06019204, 0.03666113, 0.10519847, 0.06449121, 0.05887052, 0.04450390, 0.04256770, 0.06974710, 0.0…
## $ spring_precipitation_y         <dbl> 1.640984, 3.062295, 1.857377, 1.652459, 1.757377, 3.167213, 4.555738, 3.273770, 1.880328, 2.498361, 2.093443, 4.490164, 2.516393, 3.181967, 2.867213, 2.180328,…
## $ spring_precipitation_ci_lower  <dbl> 1.234497, 2.130550, 1.370245, 1.209675, 1.174317, 2.297254, 3.519763, 2.508054, 1.387384, 1.743998, 1.483462, 3.577602, 1.853373, 2.292971, 2.084917, 1.471647,…
## $ spring_precipitation_ci_upper  <dbl> 2.047470, 3.994040, 2.344509, 2.095243, 2.340437, 4.037172, 5.591712, 4.039487, 2.373272, 3.252723, 2.703423, 5.402726, 3.179414, 4.070963, 3.649509, 2.889009,…
## $ spring_precipitation_precision <dbl> 0.09921495, 0.01888323, 0.06908390, 0.08361529, 0.04822185, 0.02166069, 0.01527467, 0.02795980, 0.06746452, 0.02880779, 0.04405942, 0.01968546, 0.03729201, 0.0…
## $ spring_nao_y                   <dbl> -0.154770488, 0.448016391, 0.131081970, 0.053229514, 1.116196720, 0.096147547, 0.102688525, -0.579950818, -0.255032787, -0.223868854, -0.423147538, -0.02888524…
## $ spring_nao_ci_lower            <dbl> -0.27969615, 0.34730133, -0.01368763, -0.05748052, 1.02540240, -0.05238801, 0.00489010, -0.69552821, -0.30866417, -0.31998757, -0.49560767, -0.11794042, 0.2328…
## $ spring_nao_ci_upper            <dbl> -0.02984483, 0.54873145, 0.27585157, 0.16393955, 1.20699104, 0.24468311, 0.20048695, -0.46437343, -0.20140141, -0.12775014, -0.35068741, 0.06016993, 0.39532357…
## $ spring_nao_precision           <dbl> 1.0504294, 1.6161488, 0.7821957, 1.3375075, 1.9886247, 0.7430350, 1.7139827, 1.2272256, 5.6994402, 1.7744110, 3.1222808, 2.0670541, 2.4837804, 2.1167932, 2.302…
## $ spring_ndvi_y                  <dbl> NA, -0.021383333, 0.005627273, -0.005341667, 0.010716667, 0.019045455, 0.010109091, -0.003100000, 0.014880000, 0.023540000, 0.013336364, 0.009625000, 0.0172416…
## $ spring_ndvi_ci_lower           <dbl> NA, -0.0351820972, -0.0040520961, -0.0185522768, -0.0018520525, 0.0115189086, 0.0018340133, -0.0106152370, 0.0074845960, 0.0139012549, -0.0002581958, 0.0010792…
## $ spring_ndvi_ci_upper           <dbl> NA, -0.0075845695, 0.0153066416, 0.0078689435, 0.0232853859, 0.0265720005, 0.0183841685, 0.0044152370, 0.0222754040, 0.0331787451, 0.0269309230, 0.0181707426, …
## $ spring_ndvi_precision          <dbl> NA, 437.6615, 889.4564, 477.4996, 527.5173, 1471.0496, 1216.9552, 1475.4802, 1523.6841, 896.9697, 450.9086, 1141.0880, 1247.4386, 1745.9602, 600.3846, 863.5844…
## $ spring_bgi_y                   <dbl> -7.185246, -4.521311, -5.604918, -1.121311, -8.009836, -5.659016, -2.106557, -3.014754, -6.926230, -6.318033, -6.345902, -4.716393, -7.545902, -5.339344, -5.27…
## $ spring_bgi_ci_lower            <dbl> -7.992711, -6.144003, -6.856932, -2.333053, -9.370581, -7.222644, -3.940855, -4.008716, -8.005576, -7.584601, -7.880432, -6.194400, -8.726326, -7.021493, -6.42…
## $ spring_bgi_ci_upper            <dbl> -6.37778129, -2.89862040, -4.35290452, 0.09043007, -6.64909112, -4.09538887, -0.27225927, -2.02079190, -5.84688316, -5.05146497, -4.81137091, -3.23838677, -6.3…
## $ spring_bgi_precision           <dbl> 0.025143351, 0.006225847, 0.010458084, 0.011164780, 0.008853535, 0.006705074, 0.004872259, 0.016593211, 0.014071763, 0.010219115, 0.006961759, 0.007504424, 0.0…

# save
save(df_abiotic_short, file = "C:/Users/quittet/Documents/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/df_abiotic_short.rds")

CORRELATION DES INDICES

On utilise le package psych. Par souci de lisibilité, seules les corrélations significatives (p \(\leq\) 0.10) sont représentées.

ACP :

On compute une ACP par saison avec l’ensemble des variables climatiques générées. L’idée étant de casser la multicollinéarité entre les variables climatiques en utilisant un score (index) annuel par saison.

Préparation des données

df_abiotic_short <- df_abiotic_short |>
  filter(year %in% 1990:2023) |> 
  as.data.frame()

winter_data_pca <- df_abiotic_short[c(
  "year",
  "winter_snow_y",
  "winter_temperature_y",
  "winter_nao_y",
  "winter_precipitation_y")]

rownames(winter_data_pca) <- winter_data_pca$year
winter_data_pca$period <- ifelse(as.numeric(as.character(rownames(winter_data_pca))) <= 2005,
                                 "P1",
                                 "P2")

spring_data_pca <- df_abiotic_short[c(
  "year",
  "spring_ndvi_y",
  "spring_temperature_y",
  "spring_precipitation_y",
  "spring_bgi_y"
 #  "spring_nao_y"
)]

rownames(spring_data_pca) <- spring_data_pca$year
spring_data_pca$period <-
ifelse(as.numeric(as.character(rownames(spring_data_pca))) <= 2005,
                                 "P1",
                                 "P2")

summer_data_pca <- df_abiotic_short[c(
  "year",
  "summer_bgi_y",
  "summer_precipitation_y",
  "summer_ndvi_y",
  "summer_temperature_y"
  # "summer_nao_y"
)]

rownames(summer_data_pca) <- summer_data_pca$year
summer_data_pca$period <- 
ifelse(as.numeric(as.character(rownames(summer_data_pca))) <= 2005,
                                 "P1",
                                 "P2")

Winter

ACP

Variances expliquées par axe et contribution des variables

##                        Dim.1 Dim.2 Dim.3 Dim.4
## winter_snow_y           42.2   2.3  20.3  35.2
## winter_temperature_y     0.2  63.4  11.1  25.3
## winter_nao_y            26.4  28.0   6.7  38.9
## winter_precipitation_y  31.2   6.3  61.9   0.6

Eigenvalues des axes

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1       1.61                  40.27                             40.27
## comp 2       1.28                  31.98                             72.26
## comp 3       0.67                  16.72                             88.97
## comp 4       0.44                  11.03                            100.00

Rappel : une eigenvalue > 1 signifique que l’axe explique plus de variance qu’une variable du jeu de données original. En général on garde toutes les composantes principales avec une eigenvalue > 1. Je confirme avec le test de simulation/rééchantillonnage ci-dessous.

Significativité des axes

## Parallel analysis suggests that the number of factors =  NA  and the number of components =  0

Je conserve l’axe 1 et l’axe 2 uniquement.

Corrélation des variables aux axes

##                        Dim.1 Dim.2 Dim.3 Dim.4
## winter_snow_y            0.8  -0.2  -0.4   0.4
## winter_temperature_y    -0.1   0.9   0.3   0.3
## winter_nao_y             0.7   0.6  -0.2  -0.4
## winter_precipitation_y   0.7  -0.3   0.6  -0.1

L’axe 1 de l’indice d’hiver oppose donc les hivers neigeux et humides avec un fort NAO aux hivers avec peu de neige et secs avec un faible NAO. L’axe 2 lui oppose les hivers chauds aux hivers froids.

Contribution des variables

# GET COORDINATES AXIS 1
ind <- factoextra::get_pca_ind(pca_winter)
index_bool <- match(df_abiotic_short$year, rownames(ind$coord))

# SAVE SEASON INDEX TO DATAFRAME
df_abiotic_short$winter_index <- ind$coord[index_bool, 1]
df_abiotic_short$winter_index_2 <- ind$coord[index_bool, 2]

Summer

Variances expliquées par axe et contribution des variables

##                        Dim.1 Dim.2 Dim.3 Dim.4
## summer_bgi_y            34.5   0.4  12.1  53.0
## summer_precipitation_y  15.7  56.7  25.4   2.1
## summer_ndvi_y           17.3  41.4  41.2   0.1
## summer_temperature_y    32.5   1.5  21.2  44.8

Eigenvalues des axes

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1       2.70                  67.57                             67.57
## comp 2       0.79                  19.72                             87.28
## comp 3       0.50                  12.54                             99.82
## comp 4       0.01                   0.18                            100.00

Significativité des axes

## Parallel analysis suggests that the number of factors =  NA  and the number of components =  1

Je conserve l’axe 1 uniquement.

Corrélation des variables aux axes

##                        Dim.1 Dim.2 Dim.3 Dim.4
## summer_bgi_y            -1.0   0.1   0.2   0.1
## summer_precipitation_y  -0.7   0.7  -0.4   0.0
## summer_ndvi_y            0.7   0.6   0.5   0.0
## summer_temperature_y     0.9   0.1  -0.3   0.1

L’axe 1 de l’ACP des variables estivales opposent les étés chauds et secs aux étés plus froids et humides.

Contribution des variables

# GET COORDINATES AXIS 1 et 2
ind <- factoextra::get_pca_ind(pca_summer)
index_bool <- match(df_abiotic_short$year, rownames(ind$coord))

# SAVE SEASON INDEX TO DATAFRAME
df_abiotic_short$summer_index <- ind$coord[index_bool, 1]

Spring

Variances expliquées par axe et contribution des variables

##                        Dim.1 Dim.2 Dim.3 Dim.4
## spring_ndvi_y            3.4  68.3  28.3   0.0
## spring_temperature_y    39.1  10.4   7.5  43.0
## spring_precipitation_y  11.5  19.1  64.2   5.1
## spring_bgi_y            46.0   2.2   0.0  51.8

Eigenvalues des axes

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1       2.13                  53.13                             53.13
## comp 2       0.99                  24.85                             77.99
## comp 3       0.88                  21.99                             99.97
## comp 4       0.00                   0.03                            100.00

Significativité des axes

## Parallel analysis suggests that the number of factors =  NA  and the number of components =  1

Je ne conserve que l’axe 1.

Corrélation des variables aux axes

##                        Dim.1 Dim.2 Dim.3 Dim.4
## spring_ndvi_y            0.3   0.8  -0.5     0
## spring_temperature_y    -0.9   0.3   0.3     0
## spring_precipitation_y   0.5   0.4   0.8     0
## spring_bgi_y             1.0  -0.1   0.0     0

Sur l’axe 1, l’indice de printemps oppose les printemps secs et chauds aux printemps doux et humides.

Contribution des variables

# GET COORDINATES AXIS 1
ind <- factoextra::get_pca_ind(pca_spring)
index_bool <- match(df_abiotic_short$year, rownames(ind$coord))

# SAVE SEASON INDEX TO DATAFRAME
df_abiotic_short$spring_index <- ind$coord[index_bool, 1]

Correlation of the indices

Save

CONCLUSION DES ANALYSES

D’après les analyses ci-dessus, on va donc conserver les axes 1 des 3 ACPs correspondant aux 3 saisons. Comme les axes 2 de chaque ACP n’expliquaient pas plus de variance qu’une simple variable, je vais conserver en plus des coordonnées de l’axe 1 les variables winter_temperature (axe 2 de l’ACP d’hiver), summer_precipitation (axe 2 de l’ACP d’été) et spring_ndvi (axe 2 de l’ACP printemps).

Ci-dessous, la représentation de la variation temporelle de l’ensemble des variables qui seront incluses dans les analyses.

Coordonnées sur l’axe 1 des ACP saisonnières

Été : valeurs positives : étés chauds et secs

Hiver : valeurs positives : NAO +, neige +, précipitations +

Printemps : valeurs négatives : printemps chauds et secs

Autres variables brutes

NOTE - Ces variables ont déjà été représentées plus haut dans le document, mais par praticité, je les intègre de nouveau.

CODE GOOGLE ENGINE

Pour voir le code, cliquer sur le bouton show ->

# > sessionInfo()
# R version 4.5.2 (2025-10-31 ucrt)
# Platform: x86_64-w64-mingw32/x64
# Running under: Windows 11 x64 (build 26100)
# 
# Matrix products: default
# LAPACK version 3.12.1
# 
# locale:
#   [1] LC_COLLATE=French_France.utf8  LC_CTYPE=French_France.utf8    LC_MONETARY=French_France.utf8 LC_NUMERIC=C                  
# [5] LC_TIME=French_France.utf8    
# 
# time zone: Europe/Paris
# tzcode source: internal
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] rgee_1.1.8
# 
# loaded via a namespace (and not attached):
#   [1] Matrix_1.7-4            jsonlite_2.0.0          compiler_4.5.2          crayon_1.5.3            Rcpp_1.1.0             
# [6] leaflet_2.2.3           jquerylib_0.1.4         png_0.1-8               fastmap_1.2.0           reticulate_1.44.1      
# [11] lattice_0.22-7          R6_2.6.1                classInt_0.4-11         sf_1.0-23               htmlwidgets_1.6.4      
# [16] leaflet.providers_2.0.0 units_1.0-0             DBI_1.2.3               pillar_1.11.1           rlang_1.1.6            
# [21] sp_2.2-0                terra_1.8-86            cli_3.6.5               withr_3.0.2             magrittr_2.0.4         
# [26] class_7.3-23            crosstalk_1.2.2         ps_1.9.1                digest_0.6.39           grid_4.5.2             
# [31] processx_3.8.6          rstudioapi_0.17.1       leafem_0.2.5            base64enc_0.1-3         rappdirs_0.3.3         
# [36] lifecycle_1.0.4         vctrs_0.6.5             KernSmooth_2.23-26      glue_1.8.0              proxy_0.4-28           
# [41] raster_3.6-32           codetools_0.2-20        e1071_1.7-16            tools_4.5.2             htmltools_0.5.9    




#' [QUELQUES REMARQUES]
#' ... [Configurer le package rgee est extrêment compliqué, j'y ai passé beaucoup d'heures]
#' ... [Cela implique de créer un compte google earth egine, de créer un projet, de configurer]
#' ... [les accès au projet etc... Des tutoriels existent sur internet, mais de manièrte générale]
#' ... [beaucoup de bugs et de messages d'erreurs. Le code suivant a fonctionné et devrait faire]
#' ... [gagner beaucoup de temps. Il est nécessaire d'installer python, et de lancer les codes] 
#' ... [commenter **une unique fois**.]


##############################################################################

#' A OUVRIR DANS UNE NOUVELLE SESSION R
#' ... Ne pas charger library(rgee) avant de lancer le code ci-dessous
#' 
##############################################################################


################################################################################
#                                                                              #
#                        CONFIGURER RGEE ET REITUCULATE                        #
#                                                                              #
################################################################################
# Reticulate sert à utiliser Python depuis R
# Rgee sert à configurer la communication entre R et google earth engine


# CONFIGURE CONDA ENVIRONMENT
# reticulate::conda_create(envname = "conda_env_rgee")  # to run once
reticulate::use_condaenv("conda_env_rgee", required = TRUE)  # remove C:/Users/quittet/AppData/Local/R/cache/R/reticulate/uv/ cache if bug
reticulate::py_config()  
reticulate::py_run_string("import ee; print(ee.__version__)")
# reticulate::py_install(
#   "earthengine-api",
#   envname = "conda_env_rgee",
#   pip = TRUE
# )  # to run once to download earthengine-api
reticulate::py_module_available("ee")  # Si TRUE, gogole engine bien installé 
# reticulate::py_run_string()  #permet de lancer un code python


# rgee::ee$Authenticate()  # run once to s'identifier
rgee::ee_clean_user_credentials(user = "pa.quittet@gmail.com", 
                                earthengine = T,
                                drive = T)  # si problème de credentials

rgee::ee_Initialize(user = "pa.quittet@gmail.com",
              project = "projetndvi-484500",
              drive = TRUE)    #' METTRE ECHAPE SI LE FICHIER ROOT EST DEJA CREE

#' NOTE : une fois que tout est configuré, il faut simplement lancer le code : 
#' ... 1. reticulate::use_condaenv("conda_env_rgee", required = TRUE)  # remove C:/Users/quittet/AppData/Local/R/cache/R/reticulate/uv/ cache if bug
#  ... 2. rgee::ee_Initialize(user = "pa.quittet@gmail.com",
#         project = "projetndvi-484500",
#         drive = TRUE)    
#  ... 3. library(rgee)



# ee <- reticulate::import("ee")  # important
# ee$Initialize(project = "projetndvi-484500")  # nécessite d'enregistrer le projet sur https://console.cloud.google.com/earth-engine/configuration?project=projetndvi-484500
# ee$Authenticate()

# SANITY CHECK
rgee::ee_check()
rgee::ee_check_python()
rgee::ee_check_credentials()
rgee::ee_check_python_packages()

# Install python package
# reticulate::py_install("numpy", envname = "conda_env_rgee", pip = TRUE)
# reticulate::py_install("pandas", envname = "conda_env_rgee", pip = TRUE)

# load python package
np = reticulate::import("numpy")       # Import Numpy 
pd = reticulate::import("pandas")      # Import Pandas


# TEST RAPIDE : doit afficher une carte du monde
library(rgee)

srtm <- ee$Image("USGS/SRTMGL1_003")
viz <- list(
  max = 4000,
  min = 0,
  palette = c("#000000","#5AAD5A","#A9AD84","#FFFFFF")
)
Map$addLayer(
  eeObject = srtm,
  visParams =  viz,
  name = 'SRTM'
)



#################################################################################
#                                                                               #
#        EXTRACTION & COMPARAISON NDVI : NOAA AVHRR – MODIS – LANDSAT           #
#                                                                               #
#                                                                               #
#  Objectif :                                                                   #
#  - Extraire des séries temporelles de NDVI sur une même zone                  #
#  - Comparer NOAA, MODIS et LANDSAT                                      #
#                                                                               #
#################################################################################

###############################################################################
# 0) PRINCIPES GÉNÉRAUX (À LIRE AVANT)
###############################################################################
# Google Earth Engine (GEE) travaille avec des objets distants :
#   
#  -  ee$ImageCollection → une collection d’images satellites
# 
#  -  ee$Image → une image satellite
# 
#  -  ee$Geometry → une géométrie (zone d’étude)
#  
#  -  Les calculs sont exécutés sur les serveurs de Google, pas localement
# 
#  -  rgee est une interface R pour écrire du code Earth Engine


##############################################################################
# DEFINITION DE LA ZONE D’ETUDE (SASSIERE)
##############################################################################

# Redéfinition sur la plaine de Milloz (zone plus homogène : pas d’eau / neige)
longitude <- c(6.975554, 6.980080, 6.979503, 6.974579)
latitude  <- c(45.488200, 45.488380, 45.490544, 45.490145)


# # Coordonnées longitude des 4 coins du polygone (ordre quelconque)
# longitude <- c(6.980740, 6.980891, 6.996077, 6.995761)
# 
# # Coordonnées latitude des 4 coins du polygone
# latitude  <- c(45.492250, 45.487339, 45.487896, 45.490588)

# Définition des bornes min et max pour créer un rectangle
lon <- c(min(longitude), max(longitude))  # [lon_min, lon_max]
lat <- c(min(latitude),  max(latitude))   # [lat_min, lat_max]

# Création de la géométrie Earth Engine (rectangle englobant)
geometry <- ee$Geometry$Rectangle(
  coords = c(
    lon[1],  # longitude minimale (Ouest)
    lat[1],  # latitude minimale (Sud)
    lon[2],  # longitude maximale (Est)
    lat[2]   # latitude maximale (Nord)
  )
)
# Possibilité d’ajouter un buffer : geometry$buffer(1000)

################################################################################
# EXTRACTION NDVI NOAA AVHRR (CDR V5)
################################################################################

noaa <- ee$ImageCollection("NOAA/CDR/AVHRR/NDVI/V5")$ # Chargement de la collection AVHRR NDVI
  filterBounds(geometry)$                             # Ne garde que les images intersectant la zone
  filterDate("1989-01-01", "2013-12-31")$             # Filtre temporel
  select("NDVI")                                      # Sélection de la bande NDVI uniquement

# Calcul du NDVI moyen par image
noaa_fc <- noaa$map(function(img) {
  
  ndvi <- img$select("NDVI")$multiply(0.0001)        # Application du facteur d’échelle
  
  mean_ndvi <- ndvi$reduceRegion(
    reducer   = ee$Reducer$mean(),                   # Calcul de la moyenne spatiale
    geometry  = geometry,                            # Zone de calcul
    scale     = 1000,                                # Résolution spatiale (AVHRR ≈ 1 km)
    maxPixels = 1e8                                  # Limite maximale de pixels autorisée
  )
  
  ee$Feature(
    NULL,                                            # Pas de géométrie associée au point
    list(
      date = img$date()$format("YYYY-MM-dd"),        # Date de l’image
      NDVI = mean_ndvi$get("NDVI")                   # NDVI moyen extrait
    )
  )
})$filter(ee$Filter$notNull(list("NDVI")))               # Suppression des valeurs NA

# Export vers Google Drive
task_noaa <- ee$batch$Export$table$toDrive(
  collection  = noaa_fc,                             # FeatureCollection à exporter
  description = "NDVI_NOAA",                          # Nom de la tâche
  folder      = "EarthEngine",                        # Dossier Drive
  fileFormat  = "CSV"                                 # Format du fichier
)

task_noaa$start()                                     # Lancement de l’export
ee_monitoring(task_noaa)                              # Suivi de l’état de la tâche

################################################################################
# EXTRACTION NDVI MODIS (MOD13A2)
################################################################################

modis <- ee$ImageCollection("MODIS/061/MOD13A2")$     # Collection MODIS NDVI 16 jours
  filterBounds(geometry)                              # Filtre spatial uniquement

modis_ndvi <- modis$map(function(img) {
  
  qa <- img$select("SummaryQA")$lt(3)                 # Masque qualité (0,1,2 acceptés)
  
  ndvi <- img$select("NDVI")$multiply(0.0001)$        # Mise à l’échelle du NDVI
    updateMask(qa)                                    # Application du masque qualité
  
  ndvi$copyProperties(
    img,
    list("system:time_start")                         # Conservation de la date
  )
})

modis_fc <- modis_ndvi$map(function(img) {
  
  stats <- img$reduceRegion(
    reducer = ee$Reducer$max()$
      combine(
        reducer2     = ee$Reducer$mean(),
        sharedInputs = TRUE
      )$
      combine(
        reducer2     = ee$Reducer$count(),
        sharedInputs = TRUE
      ),
    geometry  = geometry,
    scale     = 250,
    maxPixels = 1e8
  )
  
  ee$Feature(
    NULL,
    list(
      date        = img$date()$format("YYYY-MM-dd"),
      NDVI_max    = stats$get("NDVI_max"),
      NDVI_mean   = stats$get("NDVI_mean"),
      NDVI_count  = stats$get("NDVI_count")
    )
  )
})$filter(ee$Filter$notNull(list("NDVI_max")))


task_modis <- ee$batch$Export$table$toDrive(
  collection  = modis_fc,
  description = "NDVI_MODIS",
  folder      = "EarthEngine",
  fileFormat  = "CSV"
)

task_modis$start()
ee_monitoring(task_modis)



################################################################################
# CHARGEMENT LANDSAT ET CALCUL NDVI
################################################################################
# mask_landsat_quality <- function(img) {
#   
#   qa_pixel <- img$select("QA_PIXEL")
#   qa_radsat <- img$select("QA_RADSAT")
#   
#   mask_pixel <- qa_pixel$bitwiseAnd(1)$eq(0)   # pixel fill
#     # And(qa_pixel$bitwiseAnd(8)$eq(0))$          # cloud
#     # And(qa_pixel$bitwiseAnd(16)$eq(0))         # cloud shadow
#     # And(qa_pixel$bitwiseAnd(32)$eq(0))          # snow
#   
#   mask_radsat <- qa_radsat$eq(0)               # saturation
#   
#   img$updateMask(mask_pixel)$updateMask(mask_radsat)
# }

landsat5 <- ee$ImageCollection("LANDSAT/LT05/C02/T1_L2")$
  filterDate("1989-01-01", "2012-12-31")$
  filterBounds(geometry)
# map(mask_landsat_quality)

landsat7 <- ee$ImageCollection("LANDSAT/LE07/C02/T1_L2")$
  filterDate("1999-01-01", "2021-12-31")$
  filterBounds(geometry)
# map(mask_landsat_quality)

landsat8 <- ee$ImageCollection("LANDSAT/LC08/C02/T1_L2")$
  filterDate("2013-01-01", "2025-12-31")$
  filterBounds(geometry)
# map(mask_landsat_quality)

# Fusionner toutes les collections en une seule
landsat_all <- landsat5$merge(landsat7)$merge(landsat8)


# Calcul NDVI (formules spécifiques selon capteur)
computeNDVI_L5L7 <- function(img)
  img$normalizedDifference(c("SR_B4", "SR_B3"))$rename("NDVI")

computeNDVI_L8 <- function(img)
  img$normalizedDifference(c("SR_B5", "SR_B4"))$rename("NDVI")

landsat_fc <- landsat_all$map(function(img) {
  
  # Calcul NDVI
  ndvi <- img$normalizedDifference(c("SR_B5", "SR_B4"))$rename("NDVI")   # Landsat 5/7
  ndvi <- ee$Image(ee$Algorithms$If(img$bandNames()$contains("SR_B5"), ndvi, 
                                    img$normalizedDifference(c("SR_B5", "SR_B4"))$rename("NDVI")))
  
  # Réduire la région avec max, mean, count
  stats <- ndvi$reduceRegion(
    reducer = ee$Reducer$max()$
      combine(ee$Reducer$mean(), sharedInputs = TRUE)$
      combine(ee$Reducer$count(), sharedInputs = TRUE),
    geometry  = geometry,
    scale     = 30,
    maxPixels = 1e8
  )
  
  # Retourner un Feature (R/EE friendly)
  ee$Feature(
    NULL,
    list(
      date       = img$date()$format("YYYY-MM-dd"),
      NDVI_max   = stats$get("NDVI_max"),
      NDVI_mean  = stats$get("NDVI_mean"),
      NDVI_count = stats$get("NDVI_count")
    )
  )
})



task_landsat <- ee$batch$Export$table$toDrive(
  collection  = landsat_fc,
  description = "NDVI_LANDSAT",
  folder      = "EarthEngine",
  fileFormat  = "CSV"
)

task_landsat$start()
ee_monitoring(task_landsat)




################################################################################
# CHARGEMENT DES CSV ET PREPARATION POUR COMPARAISON
################################################################################

# Remplacer par les chemins exacts téléchargés depuis Google Drive
df_noaa  <- read.csv("C:/Users/quittet/Downloads/NDVI_NOAA.csv")
df_modis <- read.csv("C:/Users/quittet/Downloads/NDVI_MODIS.csv")
df_landsat <- read.csv("C:/Users/quittet/Downloads/NDVI_LANDSAT.csv")

df_noaa$date  <- as.Date(df_noaa$date)
df_modis$date <- as.Date(df_modis$date)
df_landsat$date <- as.Date(df_landsat$date)

df_noaa  <- df_noaa  |> mutate(year = year(date), month = month(date))
df_modis <- df_modis |> mutate(year = year(date), month = month(date))
df_landsat <- df_landsat |> mutate(year = year(date), month = month(date))


# CHECK DATA
table(df_noaa$year, df_noaa$month)
table(df_modis$year, df_modis$month)
table(df_landsat$year, df_landsat$month)

################################################################################
# AGRÉGATION MENSUELLE
################################################################################
ndvi_noaa_month <- df_noaa |>
  filter(year >= 2000 & year <= 2013) |>
  group_by(year, month) |>
  summarise(NDVI_NOAA = mean(NDVI, na.rm = TRUE), .groups = "drop")

ndvi_modis_month <- df_modis |>
  filter(year >= 2000 & year <= 2013) |>
  group_by(year, month) |>
  summarise(NDVI_MODIS = mean(NDVI_mean, na.rm = TRUE), .groups = "drop")

ndvi_landsat_month <- df_landsat |>
  filter(year >= 2000 & year <= 2013) |>
  group_by(year, month) |>
  summarise(NDVI_LANDSAT = mean(NDVI_mean, na.rm = TRUE), .groups = "drop")

################################################################################
# FUSION ET ANALYSE DE CORRELATION
################################################################################
ndvi_merged <- ndvi_noaa_month |>
  inner_join(ndvi_modis_month, by = c("year", "month")) |>
  inner_join(ndvi_landsat_month, by = c("year", "month")) |> 
  as.data.frame()

cor_noaa_modis <- cor(ndvi_merged$NDVI_NOAA, ndvi_merged$NDVI_MODIS, use = "complete.obs")
cor_noaa_landsat <- cor(ndvi_merged$NDVI_NOAA, ndvi_merged$NDVI_LANDSAT, use = "complete.obs")
cor_modis_landsat <- cor(ndvi_merged$NDVI_MODIS, ndvi_merged$NDVI_LANDSAT, use = "complete.obs")

print(paste("Corrélation NOAA–MODIS =", round(cor_noaa_modis, 3)))
print(paste("Corrélation NOAA–LANDSAT =", round(cor_noaa_landsat, 3)))
print(paste("Corrélation MODIS–LANDSAT =", round(cor_modis_landsat, 3)))

################################################################################
# VISUALISATION
################################################################################
plot(
  ndvi_merged$NDVI_MODIS, ndvi_merged$NDVI_NOAA,
  xlab = "NDVI MODIS",
  ylab = "NDVI NOAA",
  main = "Comparaison NDVI NOAA vs MODIS (NDVI=0 conservé, QA permissif)"
)
abline(lm(NDVI_NOAA ~ NDVI_MODIS, data = ndvi_merged), col = "red")

plot(
  ndvi_merged$NDVI_LANDSAT, ndvi_merged$NDVI_NOAA,
  xlab = "NDVI LANDSAT",
  ylab = "NDVI NOAA",
  main = "Comparaison NDVI NOAA vs LANDSAT"
)
abline(lm(NDVI_NOAA ~ NDVI_LANDSAT, data = ndvi_merged), col = "blue")



ndvi_noaa_month |> 
  group_by(year) |> 
  summarise(y = mean(NDVI_NOAA, na.rm = T)) |> 
  plot(type = "b")

ndvi_modis_month |> 
  group_by(year) |> 
  summarise(y = mean(NDVI_MODIS, na.rm = T)) |> 
  plot(type = "b")

ndvi_landsat_month |> 
  group_by(year) |> 
  summarise(y = mean(NDVI_LANDSAT, na.rm = T)) |> 
  plot(type = "b")

Données météo de la Grande Sassière

Pierre-Alexandre QUITTET

2026-01-14

PACKAGES

PRECIPITATION, TEMPERATURE ET EPAISSEUR NEIGEUSE

Téléchargement des données depuis METEO FRANCE

Formattage des données

Ajout des variables temporelles et trie des séquences

Filtrer les données par station d’intérêt

CHOIX DE LA STATION METEOROLOGIQUE

Critère 1 : Altitude des stations

Critère 2 : Données manquantes des séries temporelles par station

Données de précipitation (variables RR)

Données de température (variables TNTXM)

Données de neige (variables NEIGETOT06)

Visualisation des séries temporelles

Précipitations

Température

Snow

NDVI

Téléchargement du jeu de données

Computation du NDVI

NDVI issu de GOOGLE ENGINE

Visualisation des séries temporelles NDVI

NAO

Téléchargement des données depuis NOAA website

Visualisation de la série temporelle NAO

INDICE PAR SAISONS

Summer condition

Winter condition

Spring condition

GENERATE DATAFRAME

Long format

Short format

CORRELATION DES INDICES

ACP :

Préparation des données

Winter

Summer

Spring

Correlation of the indices

Save

CONCLUSION DES ANALYSES

CODE GOOGLE ENGINE

Données de précipitation (variables `RR`)

Données de température (variables `TNTXM`)

Données de neige (variables `NEIGETOT06`)