NDVI double logistic curve fitting in Python

Question

I am trying to do a NDVI double logistic curve fit in Python. This double logistic curve fitting is published by Beck et al. 2006. There is a R package greenbrown which does this over a year, so far I don't get really good or no curve fits. I would like to do the same curve fitting as this code in R: https://github.com/r-forge/greenbrown/blob/master/pkg/greenbrown/R/FitDoubleLogBeck.R

Beck et al. 2006

I have 9 NDVI observations per plot per year. This is my try out for one plot but so far I don't get really the expected results.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def double_logistic_function(t, wNDVI, mNDVI, S, A, mS, mA):
    return wNDVI + (mNDVI / (1 + np.exp(-mS * (t - S)))) * (1 - 1 / (1 + np.exp(-mA * (t - A))))

def weight_function(t, S, A, r):
    tr = 100 * (t - S) / (A - S)
    tr = np.clip(tr, 0, 100)
    return np.exp(-np.abs(r) / (1 + tr / 10))

def fit_curve(t, ndvi_observed):
    initial_guess = [np.min(ndvi_observed), np.max(ndvi_observed), np.mean(t), np.mean(t), 1, 1]
    params, _ = curve_fit(double_logistic_function, t, ndvi_observed, p0=initial_guess)
    residuals = ndvi_observed - double_logistic_function(t, *params)
    weights = weight_function(t, params[2], params[3], residuals)
    params, _ = curve_fit(double_logistic_function, t, ndvi_observed, p0=initial_guess, sigma=weights)
    return params

##Example usage with 9 observations over the year
t_9_obs = np.array([30, 60, 90, 120, 150, 180, 210, 240, 270])  

##ndvi_9_observed = np.array([0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2])
ndvi_9_observed = np.array([0.58, 0.583, 0.713, 0.807, 0.832, 0.878, 0.886, 0.863, 0.717])


##Fit the curve
parameters = fit_curve(t_9_obs, ndvi_9_observed)

##Plot the observed NDVI values
plt.scatter(t_9_obs, ndvi_9_observed, label='Observed NDVI')

##Generate points for the fitted curve
t_fit = np.linspace(min(t_9_obs), max(t_9_obs), 1000)
ndvi_fit = double_logistic_function(t_fit, *parameters)

##Plot the fitted curve
plt.plot(t_fit, ndvi_fit, label='Fitted Curve', color='red')

plt.xlabel('Day of the Year')
plt.ylabel('NDVI')
plt.legend()
plt.title('Double Logistic Curve Fitting for 9 NDVI Observations')
plt.show()

Outcome of this code:

Answer 1

I think there are two issues happening here:

1. Your mathematical model and your code do not match.
2. The optimizer is taking absurdly large steps when optimizing this function.

Your mathematical model and your code do not match, because the model is adding together the two sigmoid curves, but you are multiplying them. Also, the model uses mNDVI - wNDVI, but the code uses mNDVI.

I re-implemented your curve, using the mathematical model, and also breaking it onto several lines, so that it is easier to understand.

def double_logistic_function(t, wNDVI, mNDVI, S, A, mS, mA):
    sigmoid1 = 1 / (1 + np.exp(-mS * (t - S)))
    sigmoid2 = 1 / (1 + np.exp( mA * (t - A)))
    seasonal_term = sigmoid1 + sigmoid2 - 1
    return wNDVI + (mNDVI - wNDVI) * seasonal_term

You should always test the curve function to make sure that you are optimizing the function that you think you are optimizing.

The other issue I found was that the optimizer takes absurdly large steps for the S and A parameters.

These are the first 8 values of S, A that the optimizer tries.

S, A=(150.0, 150.0)
S, A=(150.0, 150.0)
S, A=(150.0, 150.0)
S, A=(150.00000223517418, 150.0)
S, A=(150.0, 150.00000223517418)
S, A=(150.0, 150.0)
S, A=(150.0, 150.0)
S, A=(15721.812769000722, 15722.840546778361)

When this happens, I find it is often helpful to make sure that X and Y are of similar magnitudes. To do this, I re-parameterized your X variable to vary from [0, 1] rather than [0, 365].

t_9_obs = np.array([30, 60, 90, 120, 150, 180, 210, 240, 270]) / 365

After doing this, I get a more reasonable fit.

Here is the full code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def double_logistic_function(t, wNDVI, mNDVI, S, A, mS, mA):
    sigmoid1 = 1 / (1 + np.exp(-mS * (t - S)))
    sigmoid2 = 1 / (1 + np.exp( mA * (t - A)))
    seasonal_term = sigmoid1 + sigmoid2 - 1
    return wNDVI + (mNDVI - wNDVI) * seasonal_term

def weight_function(t, S, A, r):
    tr = 100 * (t - S) / (A - S)
    tr = np.clip(tr, 0, 100)
    return np.exp(-np.abs(r) / (1 + tr / 10))

def fit_curve(t, ndvi_observed):
    initial_guess = [np.min(ndvi_observed), np.max(ndvi_observed), np.mean(t), np.mean(t), 1, 1]
    params, _ = curve_fit(double_logistic_function, t, ndvi_observed, p0=initial_guess)
    residuals = ndvi_observed - double_logistic_function(t, *params)
    weights = weight_function(t, params[2], params[3], residuals)
    params, _ = curve_fit(double_logistic_function, t, ndvi_observed, p0=initial_guess, sigma=weights)
    return params

##Example usage with 9 observations over the year
t_9_obs = np.array([30, 60, 90, 120, 150, 180, 210, 240, 270]) / 365

##ndvi_9_observed = np.array([0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2])
ndvi_9_observed = np.array([0.58, 0.583, 0.713, 0.807, 0.832, 0.878, 0.886, 0.863, 0.717])


##Fit the curve
parameters = fit_curve(t_9_obs, ndvi_9_observed)

##Plot the observed NDVI values
plt.scatter(t_9_obs * 365, ndvi_9_observed, label='Observed NDVI')

##Generate points for the fitted curve
t_fit = np.linspace(min(t_9_obs), max(t_9_obs), 1000)
ndvi_fit = double_logistic_function(t_fit, *parameters)

##Plot the fitted curve
plt.plot(t_fit * 365, ndvi_fit, label='Fitted Curve', color='red')

plt.xlabel('Day of the Year')
plt.ylabel('NDVI')
plt.legend()
plt.title('Double Logistic Curve Fitting for 9 NDVI Observations')
plt.show()

Warning: if you are interpreting the coefficients of this regression, do so with caution. There are multiple equivalent fits for this curve. For example, if you have the coefficients [wNDVI, mNDVI, S, A, mS, mA] then [wNDVI, mNDVI, A, S, -mA, -mS] is an equivalent solution that gives the same curve.