I want to apply the 'Savitzky-Golay' (savgol) filter to my time series, MODIS dataset, to remove noise (i.e., cloud pixels, etc.) in my data. MODIS have quality flags that indicates the reliability of each pixel values or if the pixel is possibly affected by clouds. So I would like to incorporate these quality flags in my filter by putting less weight or ignoring those pixel values and let savgol filter predict the optimal pixel value. I am testing np.NaN/np.nan/isnull but it seems that it removes the element in the array, and consequently savgol filter also skip those values. I would like my resulting data to be like in the attached figure.
You need to interpolate missing data before you can apply the Savitzky-Golay filter. TIMESAT is the most widely used tool for this job and they handle missing data with linear interpolation prior to applying the Savitzky-Golay filter. Assuming that you already masked cloudy and other bad observations as
np.nan here is how you can interpolate a time-series with
pandas.interpolate() and then apply the Savitzky-Golay 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)
There are of course other ways to interpolate the missing data but pandas is one of the most convenient ways to do this, especially if you want to test the effects of different interpolation algorithms.
Based on the SG filter from
scipy.signal I built the NDVI timeseries smoothing algorithm proposed in:
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, orders 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, orders 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
You're on the right track, please read the paper I've put in your question's comments.
A Flowchart with the steps you'd take from that paper (in case you don't have access) :
As you can see, for step 1 you would remove the clouds from the time series; and then apply interpolation techniques to fill the gaps.
A method would be as Chen et al. suggest linear interpolation, with methods described by Kersten, in his answer.