I'm trying to fit a stack of NDVI values to a Gaussian model to allow for determining dates of certain NDVI values using Python and NumPy/SciPy. I've attempted to do this with scipy.optimize.curve_fit.

The following is the function I'm using when applying curve_fit to the stack

def func(x, *p):
    A, mu, sigma = p
    return A * np.exp(-(x-mu)**2/(2.*sigma**2))

The parameters that I used for the curve fitting are the following

p0 = [1., 0., 1.]  # initial guess for the fitting coefficients
newX = np.linspace(date[0], date[-1], (date[-1] - date[0]) + 1)  # range of days of the year between first image in stack and last

To get the stack of values for each pixel and apply the curve fit, I used the following nested loop

for x in range(imageHeight):
    for y in range(imageWidth):
        lai = imgStack[x, y]
        popt, pcov = curve_fit(func, date, lai, p0)
        yFit = func(newX, *popt)

When I run the script, the output of yFit is always just an array of 0.0. The values of the popt variable are always the same as p0 which leads me to believe that the curve was not fit properly. Is the model not suitable for a stack of 20 NDVI values?

Credit to this post for most of this code.

  • 1
    Why do you suppose a Gaussian function ought to fit these values? I rather suspect this fitting might be doing something quite a bit different from what you were hoping for, given that you are using a date of 0 as your initial estimate of the peak (mu) and assuming a width (sigma) of just one day! For some examples and solutions check out stats.stackexchange.com/questions/11546 and stats.stackexchange.com/questions/70870. – whuber Oct 20 '14 at 19:50
  • @whuber I was suggested the Gaussian function to try and reduce the noise that occurs. I don't have very much experience in statistics so I am basically trying all the possible approaches to see which yield accurate results. – Dzinic Oct 21 '14 at 13:13
  • Not geographic? Move to stats? – Martin F Oct 21 '14 at 15:58
  • @martin That's a good idea. However, to be on topic on stats this question would have to be substantially reformulated, so I don't recommend migrating it. It usually is not constructive to start by giving a procedure that may be irrelevant to the actual research objective, noting that it is failing, and asking how to fix it up. A good question on stats (or anywhere else, for that matter) will explain what the research goal is, describe the data, and ask for specific help in achieving a clearly-stated objective. – whuber Oct 21 '14 at 17:30
  • 1
    I am with @whuber on this one. You may "reduce noise" but at the cost of biasing the resulting distribution. I would look into nonlinear smoothing functions more common in temporal analysis such as local polynomial regression (LOWESS). If you are wanting a linear fit you can explore a mixed effects model with an ARIMA-II term. – Jeffrey Evans Oct 21 '14 at 19:45

Your Answer

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

Browse other questions tagged or ask your own question.