# 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.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 score = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
};

// Function to add an NDVI band, the dependent variable.
return image
.rename('NDVI'))
.float();
};

// Function to add a time band.
// 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);
};

// Function to compute the specified number of harmonics
// and add them as bands.  Assumes the time band is present.
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);
};
};

var harmonicLandsat = landsatCollection
//.filterBounds(roi)
.filterBounds(country)

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

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

// Compute fitted values.
var fittedHarmonic = harmonicLandsat.map(function(image) {
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);

//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

})``````

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

1. Extremely good

1. 2. it isd very nices

1. 3. could you add code on how to export the data into local drive.

1. 2. The code has been updated. Please have a look. Please let me know if you have any questions.

1. 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.

4. 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?

1. Yes. Just change the time frame to 4 months like date(‘2019-01-01’, ‘2019-04-31’).