Map number of agricultural seasons

one, two or three seasons of rice?

Step 1: Import the MODIS aqua 16 days EVI data product (MYD13Q1.006 Vegetation Indices 16-Day Global 250m).

Step 2:  Set the time period

// set dates
var startdate = ee.Date.fromYMD(2010,1,1);
var enddate = ee.Date.fromYMD(2012,12,1);

Step 3:  filter for date

// select EVI
var EVI = myd.filterDate(startdate,enddate).select("EVI");

Step 4: apply scaling

// scaling factor
var multiply = function(image){
  // multply image
  var img = image.multiply(0.0001);
  return img.set('system:time_start',image.get('system:time_start')) ;

// use map function to apply function to all maps
EVI =;

Step 5: add points for harmonic trend analysis

// Define and display a FeatureCollection of three known locations.
var forest = ee.Feature(ee.Geometry.Point(104.96200561523438,19.36805056768655), {'label': 'Forest'})
var rice1 = ee.Feature(ee.Geometry.Point(105.10120957958861,19.029753537734564), {'label': 'rice1'})
var rice2 = ee.Feature(ee.Geometry.Point(105.63341617584229,18.780886954479392), {'label': 'rice2'})

var points = ee.FeatureCollection([forest,rice1,rice2])

// zoom and display

Step 6: Set variables for harmonic trend analysis

// The dependent variable we are modeling.
var dependent = 'EVI';

// The number of cycles per year to model.
var harmonics = 3;

// Make a list of harmonic frequencies to model.
// These also serve as band name suffixes.
var harmonicFrequencies = ee.List.sequence(1, harmonics);

Step 7: Construct names for harmonic trend analysis

// Function to get a sequence of band names for harmonic terms.
var getNames = function(base, list) {
  return ee.List(list).map(function(i) {
    return ee.String(base).cat(ee.Number(i).int());

// Construct lists of names for the harmonic terms.
var cosNames = getNames('cos_', harmonicFrequencies);
var sinNames = getNames('sin_', harmonicFrequencies);

Step 8: Add constant to each image as a band

// Function to add a constant band.
var addConstant = function(image) {
  return image.addBands(ee.Image(1));

// add a constant as a band
var harmonics =

Step 9: Add a time band to each image

// Function to add a time band.
var addTime = function(image) {
  // Compute time in fractional years since the epoch.
  var date = ee.Date(image.get('system:time_start'));
  var years = date.difference(ee.Date('1970-01-01'), 'year');
  var timeRadians = ee.Image(years.multiply(2 * Math.PI));
  return image.addBands(timeRadians.rename('t').float());
harmonics =

Step 10: compute harmonics

// Function to compute the specified number of harmonics
// and add them as bands.  Assumes the time band is present.
var addHarmonics = function(freqs) {
  return function(image) {
    // Make an image of frequencies.
    var frequencies = ee.Image.constant(freqs);
    // This band should represent time in radians.
    var time = ee.Image(image).select('t');
    // Get the cosine terms.
    var cosines = time.multiply(frequencies).cos()
    // Get the sin terms.
    var sines = time.multiply(frequencies).sin()
    return image.addBands(cosines).addBands(sines);

harmonics =;

Step 11: Apply regression reduction.

// Independent variables.
var independents = ee.List(['constant', 't'])

// The output of the regression reduction is a 4x1 array image.
var harmonicTrend = harmonics
  .reduce(ee.Reducer.linearRegression(independents.length(), 1));

Step 12: Turn the array image into a mult0-band image of coefficients.

// Turn the array image into a multi-band image of coefficients.
var harmonicTrendCoefficients ='coefficients')

Step 13: compute the fitted values

// Compute fitted values.
var fittedHarmonic = {
  return image.addBands(

Step 14: compute the chart for the selected points

// Plot the fitted model and the original data at the ROI.
print(Chart.image.series(['fitted', 'EVI']), forest, ee.Reducer.mean(), 250)
      title: 'Harmonic model: original and fitted Forest',
      lineWidth: 1,
      pointSize: 3,

// Plot the fitted model and the original data at the ROI.
print(Chart.image.series(['fitted', 'EVI']), rice1, ee.Reducer.mean(), 250)
      title: 'Harmonic model: original and fitted rice 1',
      lineWidth: 1,
      pointSize: 3,

// Plot the fitted model and the original data at the ROI.
print(Chart.image.series(['fitted', 'EVI']), rice2, ee.Reducer.mean(), 250)
      title: 'Harmonic model: original and fitted rice 2',
      lineWidth: 1,
      pointSize: 3,

Assignment 1:

Do the analysis for different number of trends.

Assignment 2:

Do the analysis for different locations.



Find the script here.



  1. Hi, Thank you for sharing your excellent work.
    I have a question about smoothing since I know less about Java.
    I first exclude some images with high percent of cloud, so it produce gap for time series.
    I am thinking if you could provide some other smoothing methods to fill the gap, since I feel the harmonic way tend to more overfiting.

  2. Hi thank you so much for this really helped and reduced my work.
    can you please put my coding about how to calculate VHI TCI indices in GEE.

Leave a Reply