Importation sur R des données météo

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_snow <-
  "https://www.data.gouv.fr/api/1/datasets/r/94757a8e-fa4a-4608-b846-316502983c4a"

# TELECHARGEMENT ET DECOMPRESSION DES DONNEES
df_temperature_precipitation <- data.table::fread(url_data_temperature_precipitation) 
df_snow <- data.table::fread(url_data_snow)

colnames(df_temperature_precipitation)

##  [1] "NUM_POSTE" "NOM_USUEL" "LAT"       "LON"       "ALTI"      "AAAAMMJJ" 
##  [7] "RR"        "QRR"       "TN"        "QTN"       "HTN"       "QHTN"     
## [13] "TX"        "QTX"       "HTX"       "QHTX"      "TM"        "QTM"      
## [19] "TNTXM"     "QTNTXM"    "TAMPLI"    "QTAMPLI"   "TNSOL"     "QTNSOL"   
## [25] "TN50"      "QTN50"     "DG"        "QDG"       "FFM"       "QFFM"     
## [31] "FF2M"      "QFF2M"     "FXY"       "QFXY"      "DXY"       "QDXY"     
## [37] "HXY"       "QHXY"      "FXI"       "QFXI"      "DXI"       "QDXI"     
## [43] "HXI"       "QHXI"      "FXI2"      "QFXI2"     "DXI2"      "QDXI2"    
## [49] "HXI2"      "QHXI2"     "FXI3S"     "QFXI3S"    "DXI3S"     "QDXI3S"   
## [55] "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% 1989:2025) 

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

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)

FALSE signifie une donnée manquante et TRUE une donnée disponible. 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
##   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
##   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
##   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.

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:2024){  # 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(3:9) & day_temp %in% c(1, 5, 10, 15, 20, 25)]  # Prendre quelques jours par mois

  # 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)
  
  # 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
  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"  # ne change rien, avec ou sans ce sont les mêmes valeurs
    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

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 dans la zone.

Le code est disponible dans le script rgee.R. 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% 1989:2024) |>
  mutate(year = as.character(year),
         month = as.character(month),
         day = as.character(day))

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

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

Plot function

index_plot_function <- function(data, y, ylab, title){
  data |>
    ggplot(aes(
      x = year,
      y = !!(sym(y)),
      size = precision,
      weight = precision,
      group = 1
    )) +
    geom_vline(xintercept = "2005", col = "#CD5555", linetype = 1) +
    geom_point() +
    scale_x_discrete(breaks = c(seq(1989, 2024, 2))) +
    theme_classic(base_size = 15) +
    theme(axis.text.x = element_text(angle = 90),
          legend.position = "none",
          panel.grid = element_line(color = "#CFCFCF",
                                        size = 0.5,
                                        linetype = 2)) +
    geom_smooth() +
    ylab(ylab) +
    xlab("") +
    ggtitle(title)
}

Summer condition

# MEAN PRECIPITATION
summer_mean_precipitation <- 
  df_precipitation |> 
  group_by(year) |> 
  filter(month %in% c("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 NAO
summer_nao <- 
  df_nao |> 
  group_by(year) |> 
  filter((month %in% 6 & day %in% c(15:31)) | (month %in% c(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% "06" & day %in% c(15:31)) | (month %in% c("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 |> 
  group_by(year) |> 
  filter(month %in% c("12","01","02","03")) |> 
  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)

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

# WINTER NAO
winter_nao <- 
  df_nao |> 
  group_by(year) |> 
  filter(month %in% c("12", "1", "2", "3")) |> 
  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))

Spring condition

# TEMPERATURE
spring_temperature <- 
  df_temperature |> 
  group_by(year) |> 
    filter(((month %in% "04" & day %in% c(15:31)) | (month %in% "05" & day %in% 1))) |> 
  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)

# SPRING NDVI
spring_ndvi <- 
  df_ndvi |> 
  group_by(year) |> 
    filter(((month %in% "04" & day %in% c(15:31)) | (month %in% "05" & day %in% "01"))) |> 
    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% 4 & day %in% c(15:31)) | (month %in% 5 & day %in% 1))) |> 
    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(
      year,
      variable = variable_name,
      season = season,
      value = .data[[value_col]],
      ci_lower,
      ci_upper
    )
}

# Liste des data.frames standardisés
df_long <- bind_rows(
  # Été
  standardize_df(summer_mean_precipitation, "rr", "precipitation", "summer"),
  standardize_df(summer_nao, "mean_nao", "nao", "summer"),
  standardize_df(summer_bgi, "bgi", "bgi", "summer"),
  
  # Hiver
  standardize_df(winter_temperature, "mean_temperature", "temperature", "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_ndvi, "ndvi", "ndvi", "spring"),
  standardize_df(spring_nao, "mean_nao", "nao", "spring")
)

# Aperçu
glimpse(df_long)

## Rows: 323
## Columns: 6
## $ year     <int> 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1…
## $ variable <chr> "precipitation", "precipitation", "precipitation", "precipita…
## $ season   <chr> "summer", "summer", "summer", "summer", "summer", "summer", "…
## $ value    <dbl> 1.988710, 1.504839, 1.587097, 2.814516, 2.437097, 2.879032, 1…
## $ ci_lower <dbl> 1.442727, 1.004614, 1.162954, 1.943240, 1.786062, 2.093053, 1…
## $ ci_upper <dbl> 2.534692, 2.005063, 2.011239, 3.685792, 3.088131, 3.665012, 1…

Short format

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

# Aperçu
glimpse(df_wide)

## Rows: 36
## Columns: 28
## $ year                          <int> 1989, 1990, 1991, 1992, 1993, 1994, 1995…
## $ summer_precipitation_y        <dbl> 1.988710, 1.504839, 1.587097, 2.814516, …
## $ summer_precipitation_ci_lower <dbl> 1.442727, 1.004614, 1.162954, 1.943240, …
## $ summer_precipitation_ci_upper <dbl> 2.534692, 2.005063, 2.011239, 3.685792, …
## $ summer_nao_y                  <dbl> 0.203256410, 0.146615383, 0.225551280, 0…
## $ summer_nao_ci_lower           <dbl> 0.121796160, 0.080777643, 0.125541302, 0…
## $ summer_nao_ci_upper           <dbl> 0.28471666, 0.21245312, 0.32556126, 0.17…
## $ summer_bgi_y                  <dbl> NA, -23.23065, -15.06875, -20.27660, -21…
## $ summer_bgi_ci_lower           <dbl> NA, -24.00118, -17.56434, -21.53252, -22…
## $ summer_bgi_ci_upper           <dbl> NA, -22.46011, -12.57316, -19.02067, -21…
## $ winter_temperature_y          <dbl> -1.081818, -2.847107, -3.636364, -2.8114…
## $ winter_temperature_ci_lower   <dbl> -1.379180, -3.221801, -4.056276, -3.0972…
## $ winter_temperature_ci_upper   <dbl> -0.784456, -2.472414, -3.216451, -2.5257…
## $ winter_snow_y                 <dbl> 70.04132, 85.04959, 110.05785, 118.50000…
## $ winter_snow_ci_lower          <dbl> 66.04221, 79.58947, 106.16045, 116.79924…
## $ winter_snow_ci_upper          <dbl> 74.04043, 90.50971, 113.95525, 120.20076…
## $ winter_nao_y                  <dbl> 0.54954545, 0.62576860, 0.41604132, 0.39…
## $ winter_nao_ci_lower           <dbl> 0.46828727, 0.56300726, 0.35865792, 0.33…
## $ winter_nao_ci_upper           <dbl> 0.63080363, 0.68852994, 0.47342473, 0.44…
## $ spring_temperature_y          <dbl> -0.96250, -0.34375, -2.39375, 2.30000, 3…
## $ spring_temperature_ci_lower   <dbl> -1.5421461, -1.0591152, -3.2761458, 1.14…
## $ spring_temperature_ci_upper   <dbl> -0.3828539, 0.3716152, -1.5113542, 3.459…
## $ spring_ndvi_y                 <dbl> -0.03110000, -0.00847500, -0.01165000, 0…
## $ spring_ndvi_ci_lower          <dbl> -0.057198180, -0.028293610, -0.038318755…
## $ spring_ndvi_ci_upper          <dbl> -0.005001820, 0.011343610, 0.015018755, …
## $ spring_nao_y                  <dbl> -0.07494117, 1.27323529, -0.76770587, 1.…
## $ spring_nao_ci_lower           <dbl> -0.262992396, 1.155697107, -0.841751155,…
## $ spring_nao_ci_upper           <dbl> 0.11311005, 1.39077348, -0.69366059, 1.5…