In remymoine/PTS-R-package: Reshape and Parameters for Time Series analyses

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

Le parc national des Ecrins est un espace protégé des Alpes Françaises créé en 1973. Il concerne les vallées du Briançonnais, de la Vallouise, de l'Embrunais, du Champsaur, du Valgaudemar, du Valbonnais et de l'Oisans.

Particulièrement concerné par les enjeux climatiques du fait des activités pastorales et agricoles qu'il abrite et de sa fonction de conservation de la biodiversité, il est impliqué dans plusieurs protocoles de suivi de paramètres physiques à haute fréquence.

Le package RPTS vise à faciliter l'analyse de ces données.

Mise en forme des séries temporelles avec `time_rawshape()`

Afin de comparer des séries temporelles, cette première fonction permet de formater les données dans un même format avec un même fuseau horaire. Comme le montre les fichiers exemples contenus dans le package, les données ne présentent pas un format intuitif ni de fuseau hoaraire visible. Cette information doit être recherchée en amont de l'utilisation de la fonction.

library(RPTS)
path1<-system.file("data","hobodata1.csv",package="RPTS")
head(read.table(path1,sep=",",header=T))[,2]
path2<-system.file("data","hobodata2.csv",package="RPTS")
head(read.table(path2,sep=",",header=T,skip=1))[,2]

La fonction time_rawshape() va mettre en forme la colonne contenant les données temporelles afin d'obtenir des dates de format YYYY-MM-DD HH:MM:SS au fuseau horaire UTC. Pour cela, le chemin d'accès aux fichiers doit lui être précisé via le paramètre file.

L'utilisation de cette fonction réclame une exploration préalable du fichier brut en entrée. Il est nécessaire de rassembler ou vérifier les informations suivantes:

les décimales sont marquées par des points "." (fixe),
le caractère séparateur de colonnes (sep=),
le numéro ou nom de colonne contenant les données de date (Dcol=),
le format de date original (order=),
le fuseau horaire original (tz=) et
le numéro de colonne contenant les données du paramètre mesuré (Pcol=).

La fonction OlsonNames() rassemble les chaines de caractères intelligibles par R pour les fuseaux horaires. Utilisez l'une de ces chaines pour indiquer le bon fuseau.

## Première série
data1<-time_rawshape(file=path1,Pcol=5,Dcol=2,tz="Etc/GMT+2",order="mdy_HMS_p!",sep=",",skip=1)
str(data1)
lubridate::tz(data1$Dates)

## Seconde série
data2<-time_rawshape(file=path2,Pcol=5,Dcol=2,tz="Etc/GMT+2",order="mdy_HMS_p!",sep=",",skip=1)
str(data2)
lubridate::tz(data2$Dates)

data2<-na.omit(data2)

La seconde série contient des données NA. l'utilisation de na.omit(data2) permet de les éliminer.

Métadonnées temporelles basiques avec `time_summary`

Avant de réaliser des analyses plus poussées sur les différentes séries temporelles, il est pertinent de vérifier leur compatibilité. Il est notamment nécessaire de visualiser si:

les séries temporelles sont continues (sans interruption),
leur pas de temps est identique,
leur emprise présente un recouvrement et
la présence d'éventuelles lacunes.

La fonction time_summary() analyse des données temporelles rassemblées dans une liste (les données des paramètres ne sont pas acceptées et doivent être analysées séparément). Chaque élément de la liste peut être nommé par l'utilisateur.

Si une même série présente deux pas de temps différents, la fonction renvoie un message d'erreur avec le numéro de la série concernée et l'exécution s'arrête.

Si une série temporelle présente une lacune, deux pas de temps différents sont détectés. La fonction renvoie alors le même message d'erreur avec le numéro de la série et arrête son exécution.

ex<-list(run=c("2018-11-02 12:00:00","2018-11-02 13:00:00","2018-11-02 14:00:00"),notrun=c("2018-11-02 12:00:00","2018-11-02 13:00:00","2018-11-03 12:00:00"))
try(time_summary(ex))

Les différentes métadonnées temporelles peuvent être visualisée sous forme de lignes de commande graph = FALSE ou de résultat graphique graph = TRUE.

Résultat en lignes de commande

ok<-list(Serie_Temporelle_No.1=data1[,1],
         Serie_Temporelle_No.2=data2[,1])
time_summary(ok,graph=F)

Il renvoie pour chaque série temporelle:

Steptime la valeur du pas de temps en seconde avec une approximation en unité commune,
First la date de la première mesure réalisée,
Last la date de la dernière mesure,
Period la longueur de la période mesurée en secondes avec une approximation en unités communes,
N.Data le nombre de données contenues dans la série.

Résultat graphique

ok<-list(Serie_Temporelle_No.1=data1[,1],
         Serie_Temporelle_No.2=data2[,1])
time_summary(ok,graph=T)

Il permet de visualiser:

la chronologie des séries temporelles,
leur(s) recouvrement(s) éventuel(s) et
la taille de ce(s) recouvrement(s).

Aggrégation temporelle avec `time_agpar`

La comparaison inter-périodes constitues une entrée des analyses temporelles. Cette fonction permet de réaliser une aggrégation journalière, mensuelle ou annuelle des paramètres. Elle renvoie un tableau permettant de réaliser des traitements graphiques ou statistiques ultérieurs pour chaque pas de temps.

Chaque donnée est aggrégée suivant 8 méthodes:

Mi extrait la valeur maximale pour le niveau d'aggrégation demandé,
Q1 la valeur du premier quartile,
Me la médiane,
Av la moyenne,
Q3 la valeur du troisième quartile,
Ma la valeur minimale,
Su la somme des valeurs (cumul) et
Amp l'amplitude du paramètre.

ag1m<-time_agpar(data=data1,Dcol=1,Pcol=2,resolution="month")
head(ag1m)

Ces données sont formatées pour être représentées à l'aide du package ggplot2.

library(ggplot2)

ggplot(ag1m,aes(group=Year,col=Year))+
  geom_line(aes(x=Month,y=Me))+
  geom_line(aes(x=Month,y=Mi),linetype="dotted",size=0.75)+
  geom_line(aes(x=Month,y=Ma),linetype="dotted",size=0.75)+
  geom_line(aes(x=Month,y=Q3),linetype="dashed")+
  geom_line(aes(x=Month,y=Q1),linetype="dashed")+
  geom_point(aes(x=Month,y=Me))+
  scale_fill_brewer(palette="Dark2")+
  scale_color_brewer(palette="Dark2")+
  theme(axis.text.x = element_text(angle=45,hjust=1),
          axis.text.y = element_text(face="bold"),
          panel.background = element_rect(fill="white"),
          panel.grid.major.x = element_line(color="grey"),
          panel.grid.major.y = element_line(color="grey"))

Cette figure compare l'évolution des température médiane, minimale/ maximale (pointillés) et Q1/Q3 (tirets) entre les différentes années.

library(ggplot2)

ggplot(ag1m,aes(group=Year,col=Year))+
  geom_line(aes(x=Month,y=Av))+
  geom_line(aes(x=Month,y=Mi),linetype="dotted",size=0.75)+
  geom_line(aes(x=Month,y=Ma),linetype="dotted",size=0.75)+
  geom_point(aes(x=Month,y=Av))+
  scale_fill_brewer(palette="Dark2")+
  scale_color_brewer(palette="Dark2")+
  theme(axis.text.x = element_text(angle=45,hjust=1),
          axis.text.y = element_text(face="bold"),
          panel.background = element_rect(fill="white"),
          panel.grid.major.x = element_line(color="grey"),
          panel.grid.major.y = element_line(color="grey"))

Cette seconde figure compare l'évolution des température moyenne et minimale/ maximale (pointillés) entre les différentes années.

Extraction des GDD et FDD avec `temp_gfdd()`.

Les GDD et FFD sont deux indices permettant d'appréhender la période de végétation des végétaux et l'influence de la température sur les processus biologiques. Ils sont calculés à un pas journalier à l'aide de la température moyenne journalière. Celle-ci peut être obtenu à l'aide de time_agpar() avec resolution="day".-

ag1d<-time_agpar(data=data1,resolution="day",Dcol=1,Pcol=2)
ag2d<-time_agpar(data=data2,resolution="day",Dcol=1,Pcol=2)

L'extraction des GDD et FDD doit être réalisée sur une période spécifique. La colonne Date fournit par time_agpar() est au format date accepté par la fonction temp_gfdd.

gdd<-temp_gfdd(ag2d,Pcol="Av",Dcol=1,start="2016-04-01",end="2016-09-30",igt=1,sft=-1)
head(gdd)
fdd<-temp_gfdd(ag2d,Pcol="Av",Dcol=1,start="2015-11-01",end="2016-05-15",igt=1,sft=-1)
head(fdd)

Les deux courbes peuvent également être représentées à l'aide du package ggplot2.

library(ggplot2)
library(cowplot)

pfdd=ggplot(fdd,aes(x=Date,y=FDD))+
  geom_line(col="blue")+
  scale_x_date(date_breaks="1 month")+
  xlab("Dates (pas mensuel)")+
  theme(axis.text.x = element_text(angle=45,hjust=1),
          axis.text.y = element_text(face="bold"),
          panel.background = element_rect(fill="white"),
          panel.grid.major.x = element_line(color="grey"),
          panel.grid.major.y = element_line(color="grey"))

pgdd=ggplot(gdd,aes(x=Date,y=GDD))+
  geom_line(col="red")+
  scale_x_date(date_breaks="1 month")+
  xlab("Dates (pas mensuel)")+
  theme(axis.text.x = element_text(angle=45,hjust=1),
          axis.text.y = element_text(face="bold"),
          panel.background = element_rect(fill="white"),
          panel.grid.major.x = element_line(color="grey"),
          panel.grid.major.y = element_line(color="grey"))

plot_grid(pfdd,pgdd,ncol=2,nrow=1)

Couverture neigeuse en fonction de la température avec `temp_neige()`.

La couverture neigeuse agit comme un isolant sur le sol. L'amplitude thermique est donc très réduite en cas de couverture conséquente. Cette fonction permet de détecter la présence d'une couverture neigeuse conséquente en se basant sur un seuil d'amplitude thermique définit par l'utilisateur. Ainsi dans le cas présenté ci-dessous, si la différence de température entre la température minimale journalmière et maximale journalière est inférieure à 1°C, il est probable que le sol soit recouvert de neige sous nos latitudes (45°N).

neige1<-temp_neige(ag1d,Dcol="Date",tmax="Ma",tmin="Mi",threshold = 1)

ggplot(neige1,aes(x=Dates,y=Neige))+geom_point()+
  scale_x_date(date_breaks="1 month")+
  xlab("Dates (pas mensuel)")+
  theme(axis.text.x = element_text(angle=45,hjust=1),
          axis.text.y = element_text(face="bold"),
          panel.background = element_rect(fill="white"),
          panel.grid.major.x = element_line(color="grey"),
          panel.grid.major.y = element_line(color="grey"))