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.

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


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 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)

enter image description here

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.

  • Hi Kersten, I tried the interpolate method in pandas but i cannot make it to work properly. The nan values in my array are still nan after interpolation. My NDVI values are in integer, so I am not sure if this is the issue. Here is the code that I am trying to run:\n ndvi=csvfile.irow(index)[4:50] \n ndvi[ndvi < 1900]=np.nan \n ndvi=pd.Series(ndvi) \n #x=np.arange(1,362,8) \n y=ndvi.interpolate(method="linear") – user32145 Dec 16 '15 at 6:13
  • I tried to convert it to float but the same result. The np.nan in my array is still np.nan after interpolation – user32145 Dec 16 '15 at 6:33
  • This code should work (and does at my PC) if ndvi is a Series (i.e. 1-dimensional). Maybe you can spot what is going on by consulting the pandas documentation? – Kersten Dec 16 '15 at 7:50
  • Hi Kersten, I found the culprit. So I am accessing the NDVI values from a CSV file just to test some algorithms before running it in large datasets. When accessing NDVI values from my CSV file the dtype is object and not numeric dtype I assume. Eventhough I converted it to numpy array, i.e., np.array(Series) the dtype is still an object. So I need to declare the data as np.float64(series) and it is working now. – user32145 Dec 16 '15 at 8:08

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[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:
            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) :

enter image description here

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.

  • Hi Nickves, I am thinking how can I perform step 3 (Determination of weights for each point in the NDVI time series). Is there a function already available in Python or I have to define a new one for this task? – user32145 Mar 1 '16 at 0:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.