time series: nettoyer et lisser

partons d’un jeu de données réel qui décrit des grandeurs géophysiques sur une période plus de 10 ans avec une fréquence de 1h

import pandas as pd

pd.read_excel¶

ce qui va nous donner une occasion de travailler sur un fichier .xlsx
pour cela le point d’entrée c’est pd.read_excel, mais il faut savoir que cela demande une dépendance supplémentaire

du coup si nécessaire, faites l’installation de openpyxl comme indiqué, par exemple:

# uncomment if needed
# %pip install openpyxl

chargement¶

le premier chargement “naif” avec read_excel nous permet de voir que la donnée temporelle ressemble à ceci

18/08/2006 21:00

aussi on va utiliser ce format pour charger “proprement” notre dataframe; cela va avoir deux avantages

• d’abord, on est sûr de notre affaire; on ne risque pas de laisser read_excel inventer des dates folkloriques

• ensuite, ça ne peut qu’accélérer le chargement

# c'est un peu long à charger
# aussi on va ranger ça dans un coin 
# et si nécessaire, on copiera ça au lieu de recharger

FORMAT = "%d/%m/%Y %H:%M"

try:
    df_original
    print("using preloaded data")
except NameError:
    df_original = pd.read_excel("data/DB_Galion.xlsx", date_format=FORMAT)

# repartir d'ici si on a besoin de repartir de 0:

df = df_original.copy()
df.head()

pour info, les données en question signifient ceci:

# quelques ordres de grandeur

df.info()

<class 'pandas.DataFrame'>
RangeIndex: 117188 entries, 0 to 117187
Data columns (total 5 columns):
 #   Column             Non-Null Count   Dtype         
---  ------             --------------   -----         
 0   Date_Heure_locale  117188 non-null  datetime64[us]
 1   NP                 116096 non-null  object        
 2   PA                 114921 non-null  float64       
 3   MG                 115843 non-null  float64       
 4   ET                 117188 non-null  float64       
dtypes: datetime64[us](1), float64(3), object(1)
memory usage: 4.5+ MB

ça fait tout de même 120.000 entrées de 4 mesures chacune; ça ne va pas être possible de faire une inspection visuelle de tout ça

voyons un peu le genre de chose qu’on peut faire automatiquement pour s’assurer de la complétude / cohérence de ces données

préparation¶

pour commencer, la todo list est souvent celle-ci:

• le plus urgent est presque toujours de transformer les données temporelles dans un type pertinent;
  ici on l’a déjà fait au moment du chargement, comme on le voit ci-dessus dans le résultat de df.info()

• avant de choisir un index; s’assurer de l’unicité

• adopter l’index

ici notre colonne candidat pour être l’index, c’est cette colonne temporelle, c’est vraiment naturel

unicité ?¶

# pour être propre, on veut s'assurer que la colonne en question
# est bien à valeurs uniques 

# normalement value_counts() ne devrait contenir que des 1
counts = df.Date_Heure_locale.value_counts()
dups = counts[counts != 1]

# en fait on a plusieurs entrées correspondant à ces moments-là:
dups

# pour combien de dates a-t-on 2 données

len(dups)

# voyons un peu ces multi-entrées

# on ne peut pas encore écrire directement ceci
# car justement on n'a pas encore mis la date en index
#df.loc[dups.index]

# on passe par une df temporaire qui a le bon index
# et on regarde les 3 premiers doublons

df.set_index('Date_Heure_locale').loc[dups.index[:3]]

mmh, c’est assez étrange, les 3 premières colonnes ont des données quasi identiques au moment du doublon, mais pas la 4ème...

compte tenu du faible nombre de lignes concernées, on va prendre une décision conservatoire, qui est de garder, pour chacun de ces 13 instants, la moyenne des deux mesures

sauf que, si on essaie de faire ça ce stade, on est confronté au fait que la colonne NP est de type object...

une colonne de type object¶

df.dtypes

pour bien trouver tous les soucis avec cette colonne:

# voyons ce que donne pd.to_numeric

try:
    pd.to_numeric(df.NP)
except Exception as exc:
    print(f"OOPS {type(exc)} {exc=}")

OOPS <class 'ValueError'> exc=ValueError('Unable to parse string " " at position 1695')

ce qui nous indique qu’il y a au moins un endroit où la colonne NP contient une chaine composée d’un espace; voyons combien il y en a de ce genre

(df.NP == ' ').value_counts()

fort heureusement on peut forcer la conversion comme ceci

pd.to_numeric(df.NP, errors='coerce')

cette fois la conversion se fait correctement; il ne faut pas oublier par contre de bien adopter la nouvelle colonne - car avec la cellule précédente on a calculé un nouvelle colonne mais elle ne fait pas partie de la dataframe

# comme on est satisfait on remplace la colonne dans la dataframe

df['NP'] = pd.to_numeric(df.NP, errors='coerce')

# et du coup maintenant on a bien des nombres partout

df.dtypes

la moyenne¶

on peut à présent faire la moyenne pour les doublons

df = df.groupby('Date_Heure_locale')[['NP', 'PA', 'MG', 'ET']].mean()
df.head()

# et pour être bien sûr

# maintenant value_counts() ne doit contenir que des 1
counts = df.reset_index().Date_Heure_locale.value_counts()
dups = counts[counts != 1]

if len(dups) == 0:
    print("OK !")

OK !

c’est terminé¶

et on n’a même pas besoin d’adopter l’index puisque ça a déjà été fait par le groupby..

df.head(2)

aperçu¶

notre objectif: passer le moins de temps possible pour voir une vague idée de ces données

# version simplissime: vraiment pas top

df.plot();

la première amélioration vient au prix d’une ligne qu’on mentionne généralement au début (autour de par exemple la ligne import pandas as pd)

%matplotlib ipympl pour des graphiques interactifs¶

%matplotlib ipympl

# une fois qu'on a choisi ce mode on obtient des visus interactives
# on peut agrandir la figure, zoomer, se déplacer, etc...

df.plot();

figsize¶

si vous voulez choisir une taille par défaut pour les figures

# il y a plein d'options pour faire ça
# j'aime bien celle-ci, mais bon...

from IPython.core.pylabtools import figsize
figsize(8, 6)

# et à partir de là...

df.plot();

# sachant qu'on peut toujours choisir la taille pour une figure donnée

df.plot(figsize=(3, 3));

subplots=True¶

comme les échelles ne sont pas forcément les mêmes, ou pour plein d’autres raisons, on peut avoir envie de voir les données indépendamment les unes des autres

df.plot(subplots=True);

alpha=0.1¶

si on insiste pour voir les 4 données dans la même vue, comme on avait fait pour commencer, il y a tellement de données en X qu’on ne voit que la dernière colonne !
dans ces cas-là le canal alpha (la transparence) est notre meilleur allié

df.plot(alpha=0.1);

les trous dans la donnée (aka nan)¶

à ce stade on a donc des nan - au moins dans NP, ça on en est sûr, et dans les autres colonnes peut-être aussi

# dans une colonne en particulier
df.NP.isna().sum()

# si on veut une information plus globale, par exemple
df.isna().sum(axis=0)

comment les calculer ?¶

ce qu’on aimerait savoir maintenant, c’est est-ce que les ‘trous’ sont plutôt éparpillés ou plutôt groupés; et pour ça on va

• fabriquer une série qui contient seulement les timestamps où on a une mesure (en enlevant ceux qui correspondent à un nan, donc)

• puis on fera la différence avec la valeur immédiatement voisine (pour cela on décale les valeurs de 1 cran, et on fait la différence)

si on avait une série parfaitement pleine, on n’obtiendrait que des valeurs de “1h” dans ce résultat
et donc en éliminant de cela les valeurs “1h”, on obtient la liste des trous dans la donnée (c’est-à-dire quand y a-t-il eu un trou, et combien de temps a-t-il duré)

# par exemple avec la colonne NP qui a les plus beaux trous
NP = df.NP.copy()

# on isole les lignes qui ont une valeur
NP_defined = NP[NP.notna()]

# ce qui nous intéresse ce sont les timestamps
NP_times = NP_defined.index

# combien de mesures en tout, combien de mesures pertinentes
NP.shape, NP_times.shape

# la grosse astuce consiste à mettre les timestamps comme valeurs
# et pour ça on utilise la fonction pd.to_series 
# souvenez-vous que NP_times est un objet Index

NP_times_series = NP_times.to_series()
NP_times_series

on décale¶

# et du coup maintenant on peut décaler (avec shift()) les valeurs de 1 cran

NP_times_series.shift()

# et surtout faire la différence entre les deux, qui normalement doit donner 1h

delta = NP_times_series - NP_times_series.shift()

delta.head(3)

# il faut enlever la première ligne, évidemment
# car le shift y a mis NaT

delta = delta.iloc[1:]

digression: comment fabriquer un timedelta¶

# c'est assez simple, on fait par exemple

pd.to_timedelta("01:00:00")

# ou encore, c'est assez flexible:

pd.to_timedelta(1, 'h')

reprenons¶

nous en sommes donc au stade où on veut enlever de delta les valeurs correspondant à 1h

holes = delta[delta != pd.to_timedelta(1, 'h')]

print(f"nous avons trouvé {len(holes)} trous dans la colonne NP")
holes

nous avons trouvé 76 trous dans la colonne NP

on les affiche¶

# si on veut les montrer
# c'est mieux de transformer en heures (sinon on a des nano-secondes apparemment)
# en plus comme on n'a pas une dataframe, avec le driver %matlotlib ipympl
# il faut créer une nouvelle figure

import matplotlib.pyplot as plt

plt.figure()
(holes / pd.to_timedelta(1, 'h')).plot(style=['r.']);

question¶

utilisez le zoom pour trouver à quelle époque a eu lieu la coupure la plus longue dans l’acquisition de NP

# on peut le retrouver comme ceci

when = holes.idxmax()
how_long = holes.loc[when]

print(f"the longest outage was {how_long} long and occurred on {when}")

the longest outage was 79 days 15:00:00 long and occurred on 2012-06-26 14:00:00

un peu de lissage¶

on regarde cette fois la courbe ‘PA’; utilisez le zoom pour regarder environ un mois

plt.figure()
PA = df.PA
PA.plot(alpha=0.4);

on voit que la fluctuation en tendance est obscurcie par des variations de plus haute fréquence; on va essayer de lisser ce signal pour ne conserver que la tendance de fond

normalement si vous zoomez encore davantage, vous devriez percevoir que la haute fréquence est de l’ordre de la demie journée

rolling¶

on va calculer le rolling de ce signal avec deux périodes (12h et 24h), et afficher tout cela

# rolling windowed: on essaie avec ces deux durées
# souvenez-vous que l'on a une mesure par heure

smooth24 = PA.rolling(24, center=True).mean()
smooth12 = PA.rolling(12, center=True).mean()

plt.figure()
PA.plot(alpha=0.2)

smooth12.plot(style=['g'], alpha=0.5);
smooth24.plot(style=['r'], alpha=0.5);

à nouveau, je vous invite à zoomer dans le graphique pour juger de la pertinence de chacune de ces deux approximations

laquelle garderiez-vous ?