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, date[-1], (date[-1] - date) + 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.