Reconstruction de séries chronologiques MODIS en appliquant le filtre Savitzky-Golay avec Python/Numpy

Vous devez interpoler les données manquantes avant de pouvoir appliquer le filtre Savitzky-Golay. TIMESAT est l'outil le plus largement utilisé pour ce travail et ils traitent les données manquantes avec une interpolation linéaire avant d'appliquer le filtre Savitzky-Golay. En supposant que vous avez déjà masqué les observations nuageuses et autres mauvaises comme np.nan voici comment vous pouvez interpoler une série chronologique avec pandas.interpolate() puis appliquer le filtre Savitzky-Golay scipy.signal.savgol_filter() .

import numpy as np
import pandas as pd
from scipy.signal import savgol_filter

#create a random time series
time_series = np.random.random(50)
time_series[time_series < 0.1] = np.nan
time_series = pd.Series(time_series)

# interpolate missing data
time_series_interp = time_series.interpolate(method="linear")

# apply SavGol filter
time_series_savgol = savgol_filter(time_series_interp, window_length=7, polyorder=2)

Il existe bien sûr d'autres moyens d'interpoler les données manquantes, mais pandas est l'un des moyens les plus pratiques de le faire, surtout si vous souhaitez tester les effets de différents algorithmes d'interpolation.

Basé sur le filtre SG de scipy.signal J'ai construit l'algorithme de lissage des séries temporelles NDVI proposé dans :

Une méthode simple pour reconstruire un ensemble de données de séries chronologiques NDVI de haute qualité basée sur le filtre Savitzky-Golay", Jin Chen et al. 2004

import pandas as pd
import numpy as np
from scipy.signal import savgol_filter
def savitzky_golay_filtering(timeseries, wnds=[11, 7], orders=[2, 4], debug=True):                                     
    interp_ts = pd.Series(timeseries)
    interp_ts = interp_ts.interpolate(method='linear', limit=14)
    smooth_ts = interp_ts                                                                                              
    wnd, order = wnds[0], orders[0]
    F = 1e8 
    W = None
    it = 0                                                                                                             
    while True:
        smoother_ts = savgol_filter(smooth_ts, window_length=wnd, polyorder=order)                                     
        diff = smoother_ts - interp_ts
        sign = diff > 0                                                                                                                       
        if W is None:
            W = 1 - np.abs(diff) / np.max(np.abs(diff)) * sign                                                         
            wnd, order = wnds[1], orders[1]                                                                            
        fitting_score = np.sum(np.abs(diff) * W)                                                                       
        print it, ' : ', fitting_score
        if fitting_score > F:
            break
        else:
            F = fitting_score
            it += 1        
        smooth_ts = smoother_ts * sign + interp_ts * (1 - sign)
    if debug:
        return smooth_ts, interp_ts
    return smooth_ts