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 = EVI.map(multiply);

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
Map.addLayer(points);
Map.centerObject(points,8);

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 = EVI.map(addConstant)

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 = harminics.map(addTime)

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()
      .rename(cosNames);
    // Get the sin terms.
    var sines = time.multiply(frequencies).sin()
      .rename(sinNames);
    return image.addBands(cosines).addBands(sines);
  };
};

harmonics = harmonics.map(addHarmonics(harmonicFrequencies));

Step 11: Apply regression reduction.

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

// The output of the regression reduction is a 4x1 array image.
var harmonicTrend = harmonics
  .select(independents.add(dependent))
  .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 = harmonicTrend.select('coefficients')
  .arrayProject([0])
  .arrayFlatten([independents]);

Step 13: compute the fitted values

// Compute fitted values.
var fittedHarmonic = harmonics.map(function(image) {
  return image.addBands(
    image.select(independents)
      .multiply(harmonicTrendCoefficients)
      .reduce('sum')
      .rename('fitted'));
});

Step 14: compute the chart for the selected points

// Plot the fitted model and the original data at the ROI.
print(Chart.image.series(fittedHarmonic.select(['fitted', 'EVI']), forest, ee.Reducer.mean(), 250)
    .setOptions({
      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(fittedHarmonic.select(['fitted', 'EVI']), rice1, ee.Reducer.mean(), 250)
    .setOptions({
      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(fittedHarmonic.select(['fitted', 'EVI']), rice2, ee.Reducer.mean(), 250)
    .setOptions({
      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.

 

One thought on “Map number of agricultural seasons”

  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.

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s