Develop a Harmonic Model: Original vs Fitted NDVI values and Export the Output to your Drive(Study Area: Thailand)

The Normalized Difference Vegetation Index (NDVI) is a simple graphical indicator that can be used to analyze remote sensing measurements, typically, but not necessarily, from a space platform, and assess whether the target being observed contains live green vegetation or not.

The NDVI is calculated from these individual measurements as follows:

NDVI= (NIR-Red) \ (NIR+Red)

There are a number of derivatives and alternatives to NDVI that have been proposed in the scientific literature to address the limitations of NDVI. These alternatives attempt to include intrinsic correction(s) for one or more perturbing factors. In this tutorial, we won’t be looking at these alternatives. Instead, we will see how NDVI may act differently at the same location in different time (year). For this we chose a location (101.548, 16.523) in Thailand. The steps are:

  1. Using the Landsat 8 images, compute the NDVI values
  2. Develop a harmonic model for original NDVI values
  3. Develop a harmonic model for fitted NDVI values
  4. Overlay these models to see the difference.

Code disclaimer: This code was originally published in Google Earth Engine under example section. The code was modified and edited so as to enable it to calculate NDVI at point (101.548, 16.523) in Thailand.

Here is the code:

var countries = ee.FeatureCollection("ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw")
var country_name = ['Thailand']
var country = countries.filter(ee.Filter.inList('Country', country_name));
//Map.addLayer(country,{},"Nepal");
Map.centerObject(country,5);  //Zoom to Study area

// Load a collection of Landsat TOA reflectance images..
var landsatCollection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterBounds(country);

// Set the region of interest to a point.
var roi = ee.Geometry.Point([101.548, 16.523]);

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

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

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

// Function to get a sequence of band names for harmonic terms.
var constructBandNames = 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 = constructBandNames('cos_', harmonicFrequencies);
var sinNames = constructBandNames('sin_', harmonicFrequencies);

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

// Function to mask clouds in Landsat 8 imagery.
var maskClouds = function(image) {
  var score = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var mask = score.lt(10);
  return image.updateMask(mask);
};

// Function to add an NDVI band, the dependent variable.
var addNDVI = function(image) {
  return image
    .addBands(image.normalizedDifference(['B5', 'B4'])
    .rename('NDVI'))
    .float();
};

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

// 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);
  };
};

// Filter to the area of interest, mask clouds, add variables.
var harmonicLandsat = landsatCollection
  //.filterBounds(roi)
  .filterBounds(country)
  .map(maskClouds)
  .map(addNDVI)
  .map(addDependents)
  .map(addHarmonics(harmonicFrequencies));

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

// Turn the array image into a multi-band image of coefficients.
var harmonicTrendCoefficients = harmonicTrend.select('coefficients')
  .arrayProject([0])
  .arrayFlatten([independents]);

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

// Plot the fitted model and the original data at the ROI.
print(ui.Chart.image.series(fittedHarmonic.select(['fitted','NDVI']), roi, ee.Reducer.mean(), 30)
    .setOptions({
      title: 'Harmonic model: original and fitted values',
      lineWidth: 1,
      pointSize: 3,
}));



// Pull out the three bands we're going to visualize.
var sin = harmonicTrendCoefficients.select('sin_1');
var cos = harmonicTrendCoefficients.select('cos_1');

// Do some math to turn the first-order Fourier model into
// hue, saturation, and value in the range[0,1].
var magnitude = cos.hypot(sin).multiply(5);
var phase = sin.atan2(cos).unitScale(-Math.PI, Math.PI);
var val = harmonicLandsat.select('NDVI').reduce('mean');


// Turn the HSV data into an RGB image and add it to the map.
var seasonality = ee.Image.cat(phase, magnitude, val).hsvToRgb();
Map.centerObject(country, 5);
Map.addLayer(seasonality.clip(country), {}, 'Seasonality');
Map.addLayer(roi, {}, 'ROI');

//Exporting the final output to your drive
var projection = "EPSG:5070"; //" "EPSG:4269"; 
var filename = "Harmonic Model (Seasonality) for "+country_name;
Export.image.toDrive({
  image: seasonality,
  description: filename,
  folder: 'Harmonic Model',
  scale: 30,
  region: country,
  crs: projection,
  maxPixels: 800000000
  
})
Seasonality – Thailand

22 thoughts on “Develop a Harmonic Model: Original vs Fitted NDVI values and Export the Output to your Drive(Study Area: Thailand)”

  1. Pingback: Using Hansen Global Forest Change Data to determine Forest Cover, Forest Gain, and Forest Loss (Study Area: Malaysia) – DINESH SHRESTHA

  2. Pingback: Make your map look Artistic – Add Grid Lines in your Map (Study Area: Canada) – DINESH SHRESTHA

      1. Its showing still downloading since last two days. Not sure how to speed up the process

        Also does this product mean ndvi during 1970 – 2018?

        var years = image.date().difference(‘1970-01-01’, ‘year’);

        How to prepare annual or quarter year product?

      2. I am sorry to hear that it took you so long. The processing time depends on different variables including the drive storage, file size, computing speed, and more. I have heard people waiting for 15 days in order to get their files exported. Hoping this is not the case here.
        As far as the years are concerned, the Landsat 8 dates back to 2014 only. No matter what dates you put there it will go back until 2014 only. You might need to use Landsat 7 and Landsat 5 images for your project. I believe Landsat 5 dates back until 1984. Be mindful that while using two or more sensors, you might have to run some algorithms to sender corrections and stuffs like that.
        In response to how to prepare annual product, I will be posting another tutorial that allows you to calculate annual 90th percentile NDVI from 1984 through 2019. Stay tuned.

  3. Pingback: Night Light Maps (Study Area: South Asia) – DINESH SHRESTHA

  4. Pingback: How to generate Histogram in Google Earth Engine? – DINESH SHRESTHA

  5. Thanks for everything. I would like to built NDVI composite layer for every four months of a year. any suggestion on the modification of the code?

  6. Pingback: Using 2001 MODIS Land Cover Type Yearly Global 500m Data to Determine the Land-Cover Analysis (Study Area: Mexico) – DINESH SHRESTHA

  7. Pingback: Learn How to Add two charts next to the map to interactively display a time-series of NDVI and reflectance for each click on the map? – DINESH SHRESTHA

  8. Pingback: Create a Chart to show Landsat 8 TOA Spectra at three points near Mumbai City – DINESH SHRESTHA

  9. Pingback: Learn how to calculate and export 90th Percentile Annual NDVI for years 1985-2019 using Landsat 8, 7, and 5 scenes (Study Area: Turkey) – DINESH SHRESTHA

  10. Pingback: Use 2018 Landsat 8 Scenes for Land Cover Classification in Google Earth Engine (Study Area: Dubai) – DINESH SHRESTHA

  11. Pingback: Land Cover Change Analysis between 1999 and 2018 in GEE (Study Area: Mumbai) – DINESH SHRESTHA

  12. Pingback: Time-lapse: Animate 30m Landsat images generated 90th percentile Annual NDVI from 1999 to 2018 (Study Area: South America) – DINESH SHRESTHA

Leave a Reply

%d bloggers like this: