Reconstrucción de series temporales MODIS aplicando el filtro Savitzky-Golay con Python/Numpy

Question

Quiero aplicar el filtro 'Savitzky-Golay' (savgol) a mi serie temporal, el conjunto de datos MODIS, para eliminar el ruido (es decir, los píxeles de las nubes, etc.) en mis datos. MODIS tiene banderas de calidad que indican la fiabilidad de los valores de cada píxel o si el píxel está posiblemente afectado por las nubes. Así que me gustaría incorporar estas banderas de calidad en mi filtro poniendo menos peso o ignorando esos valores de píxeles y dejar que el filtro savgol prediga el valor óptimo del píxel. Estoy probando np.NaN / np.nan / isnull pero parece que elimina el elemento de la matriz, y en consecuencia el filtro savgol también omite esos valores. Me gustaría que mis datos resultantes fueran como en la figura adjunta.

( https://matinbrandt.wordpress.com/2014/12/02/smoothingfiltering-a-ndvi-time-series-using-a-savitzky-golay-filter-and-r/ )

Answer 1

Es necesario interpolar los datos que faltan antes de aplicar el filtro Savitzky-Golay. TIMESAT es la herramienta más utilizada para este trabajo y manejan los datos faltantes con interpolación lineal antes de aplicar el filtro Savitzky-Golay. Suponiendo que ya ha enmascarado las observaciones nubladas y otras malas como np.nan así es como se puede interpolar una serie temporal con pandas.interpolate() y luego aplicar el filtro 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)

Por supuesto, hay otras formas de interpolar los datos que faltan, pero pandas es una de las formas más convenientes de hacerlo, especialmente si quiere probar los efectos de diferentes algoritmos de interpolación.

Answer 2

Basado en el filtro SG de scipy.signal Construí el algoritmo de suavización de series temporales de NDVI propuesto en:

A simple method for reconstructing a high quality NDVI time-series data set based on the Savitzky-Golay filter", 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

Answer 3

Vas por buen camino, por favor lee el documento que he puesto en los comentarios de tu pregunta.

Un diagrama de flujo con los pasos que daría ese papel (por si no tienes acceso) :

Como puede ver, en el paso 1 se eliminan las nubes de la serie temporal y se aplican técnicas de interpolación para rellenar los huecos.

Un método sería como sugieren Chen et al. la interpolación lineal, con los métodos descritos por Kersten, en su respuesta.

Answer 4

El código python de arriba contestado por Pateheo es bastante útil, pero la W no parece actualizarse a medida que itera. Así es como funciona según el documento

Además, la puntuación de ajuste se calcula antes de asignar el "smoother_ts", lo que hace que se sumen valores negativos a la puntuación de ajuste. ¿Qué tal si se cambia el código de la siguiente manera?

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).interpolate(method='linear', limit=9999)
    smooth_ts = interp_ts                                                                                              
    wnd, order = wnds[0], orders[0]
    F = 1e8 
    it = 0                                                                                                             
    while True:
        smoother_ts = savgol_filter(smooth_ts, window_length=wnd, polyorder=order)                                     
        diff = (smoother_ts - interp_ts).clip(0)
        sign = diff > 0                                                                                                                       
        W = 1 - diff / np.max(diff)                                                         
        wnd, order = wnds[1], orders[1]                                                                            
        fitting_score = np.sum(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