I am attempting to optimize the parameters for a double-logistic function on an annual MODIS NDVI time series in Google Earth Engine. I have previously looked at Non Linear Regression in Google Earth Engine on FeatureCollection, which is a similar problem without much clarity, other than hinting that it is possible to optimize the parameters for nonlinear regression within GEE. Previously, I have been able to use Python to change each image of the collection into an array, stack the arrays, then subset each length 46 array from each pixel and fit my function using lmfit.

However, I run into memory issues when attempting to stack arrays when using it on more than the 1x1 degree block of images I initially got it working on. So, my thought is that I need to do this without manipulating the Image Collection to a stack of np.arrays as that then attempts to load the Image Collection into local memory.

-Is this only a problem that can be solved by cutting the arrays into chunks in Python, or can I leverage GEE to reduce my computing times, beginning with the script below?

//Date range for image collection
var date_start = ee.Date('2002-01-01');
var end_date = ee.Date('2002-12-31');

//Region of interest
var roi = ee.Geometry.Rectangle(-100.1, 40, -100, 40.1);

//Clipping function for image collection
var clipping = function(image){
  return image.clip(roi)};

// define the image collection
var collection = ee.ImageCollection('MODIS/006/MOD09Q1')
    .filterDate(date_start, end_date)

// NDVI function
var NDVIcalc = function(image){
  var red = image.select('sur_refl_b01');
  var nir = image.select('sur_refl_b02');
  var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
  return ndvi;

//apply NDVI function to collection
var series = collection.map(NDVIcalc);

//6 parameter double logistic fundtion from Beck et al. 2006 where; t = doy, wNDVI = Annual minimum NDVI, 
//wNDVI = Annual maximum NDVI, mS = Spring slope, S = SOS, mA = Fall slope, A = EOS
var dbl_logistic_model = function(t, wNDVI,mNDVI,mS,S,mA,A){
    return wNDVI + (mNDVI - wNDVI)* ( 1/(1+Math.exp(-mS*(t-S))) + 1/(1+Math.exp(mA*(t-A))) - 1)};

//Function to create range of numbers
var Arange= function(a, b, step){
    var A= [];
        A[0]= a;
        step= step || 1;
        while(a+step<= b){
            A[A.length]= a+= step;
    return A;

//Initial parameter guesses
var params = [0.1,0.7,0.1,100,0.1,250]
//create day of year variable for "t" in model
var xdata = Arange(1,361,8);

//In python here for a length 46 array within a single cell the normal process here I would take is to wrap dbl_logistic_function in lmfit.Model and use the .fit function from lmfit to optimize parameters
//Where somefunction uses an optimizer like Levenberg-Marquardt to optimize parameters for each pixel, and mod_fit is the 6 parameter output
var mod_fit = somefunction(series, x=xdata, params)
  • The example you pointed is suggesting that you can transform your inputs non-linearly, but you can still only do linear regression with them (ie: you can regress over exp(t) + t^2 + log(t^3), but that's still a linear regression). In order to be able to use some form of python-based non-linear regression, you'll have to get all the values into your python client (perhaps using reduceRegion with a toList() or toCollection reducer, and getInfo() the result), which isn't going to scale too well. – Noel Gorelick Dec 29 '18 at 18:25

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.