8

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

1

4 Answers 4

6

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.

4
  • 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
    Commented Dec 16, 2015 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
    Commented Dec 16, 2015 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
    Commented Dec 16, 2015 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
    Commented Dec 16, 2015 at 8:08
5

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:
            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
1

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.

1
  • 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
    Commented Mar 1, 2016 at 0:34
0

The python code above answered by Pateheo is quite useful, but the W does not seem to be updated as it iterates. This is how it works according to the paper paper_w

Also, the fitting score is calculated before the "smoother_ts" is assigned, which results in negative values to be summed to the fitting score. How about changing the code as follows:

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
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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