7 MODIS Time-Series Analysis using Google Earth Engine (2008–2024)
7.1 Introduction
While Sentinel-2 provides high spatial resolution (10m), the MODIS sensor allows us to look much further back in time with high temporal frequency. This chapter demonstrates how to analyze vegetation dynamics across Nepal over a 16-year period.
7.2 Data Selection: MODIS MOD09Q1
We use the MODIS 250m 8-day composite. This product is ideal for monitoring regional changes without the heavy processing load of high-resolution data[cite: 7].
| Band | GEE Name | Description |
|---|---|---|
| B1 | sur_refl_b01 | Red (620-670 nm) |
| B2 | sur_refl_b02 | NIR (841-876 nm) |
7.3 Uploading Your Study Area (Shapefile)
To use your own boundaries (like a specific District or Watershed) in GEE, follow these steps to upload them as an Asset:
- Prepare your File:
- On your computer, select all files related to your shapefile.
- Right-click and Compress them into a single .zip file.
- Open GEE Assets:
- In the Google Earth Engine Code Editor, click the Assets tab (top-left panel).
- Click the red NEW button and select Table Upload Shape files.
- Upload:
- Select your
.zipfile from your computer. - Give it an Asset ID (e.g.,
nepal_study_area) and click Upload.
- Import to Script:
- Once the “Tasks” tab shows a blue checkmark (it may take 1–2 minutes), go back to Assets.
- Click the Arrow icon () next to your asset to import it into your script as
table.
Note: GEE will reject the upload if the (projection) file is missing. Always ensure your shapefile is in WGS 84 (EPSG:4326) for best compatibility.
7.3.1 Before starting to code: GEE Terms and their meanings
| Technical Term | Plain English Translation | What it does in your exercise |
|---|---|---|
ee.Image |
A single snapshot. | A single map of Nepal from one specific week. |
ee.ImageCollection |
A photo album. | All the snapshots from 2008 to 2024 stacked together. |
ee.FeatureCollection |
A digital stencil. | The boundary of Nepal or its provinces used to “cookie-cut” the data. |
.map() |
The “Do this to all” command. | Tells GEE: “Calculate NDVI for every single year in my album at once.” |
.clip() |
The Cookie Cutter. | Cuts out the rest of the world so you only see Nepal. |
7.4 Automated Time-Series Recipe
7.4.1 1. Defining the Period and Area
Instead of downloading dozens of files, we define a list of years and a function to “fetch” them automatically.
var districts = ee.FeatureCollection('projects/ee-00001-484109/assets/modis/nepal_district');
var years = [2008, 2012, 2016, 2020, 2024];
var nepalRoi = ee.Geometry.Rectangle([80.0, 26.0, 89.0, 31.0]);
Map.centerObject(districts, 7);
Map.addLayer(districts, {color: 'red'}, 'Nepal Districts');
// 3. Function: Get MODIS composite for a year
function getModisImage(year) {
// Find the 8-day composite that overlaps your desired period (e.g., early May)
// Example: Target composites around day-of-year ~129 (May 9 start, but adjust)
var start = ee.Date.fromYMD(year, 1, 1);
var end = ee.Date.fromYMD(year + 1, 1, 1);
var collection = ee.ImageCollection('MODIS/061/MOD09Q1')
.filterDate(start, end)
.filterBounds(nepalRoi);
// Select one composite per year — e.g., the one closest to mid-May
// Or manually pick by composite_date property if needed
var targetDay = 129;
var image = collection
.filter(ee.Filter.calendarRange(targetDay - 8, targetDay + 8, 'day_of_year'))
.sort('system:time_start', false) // Latest first
.first()
.select(['sur_refl_b01', 'sur_refl_b02']) // Same band names!
.clip(districts.geometry());
return image.set('year', year);
}
// 4. Process all years
var modisImages = years.map(function(year) {
return getModisImage(year);
});
// Convert to ImageCollection for easier handling (optional but useful)
var modisCollection = ee.ImageCollection.fromImages(modisImages);
7.4.2 Breaking Down the Code
var modisCollection = ee.ImageCollection.fromImages(years.map(modisImages));
The Input (
years): We start with a simple list of years we care about.The Action (
.map+getModis): We tell the computer to take that list and, for every year, run our “Function” (getModis). This funnction finds the satellite data, cuts it to the Nepal border, and labels it.The Result (
fromImages): We take all those individual processed maps and stack them into one organized “Image Collection” ready for analysis.
7.4.3 2. Temporal NDVI Calculation
We apply the NDVI formula to every image in our collection simultaneously using the .map() function.
// 5. Calculate NDVI
function calculateNdvi(image) {
var red = image.select('sur_refl_b01').multiply(0.0001);
var nir = image.select('sur_refl_b02').multiply(0.0001);
var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
return image.addBands(ndvi);
}
var ndviCollection = modisCollection.map(calculateNdvi);7.5 Counting Number of pixels (or areas)
In Google Earth Engine, .reduceRegions() is the primary tool for performing Zonal Statistics. It allows us to take a layer of data (like NDVI) and summarize it based on geographic boundaries (like Provinces or Districts).
We can consider it as an “Information Extractor” that turns pixels into a spreadsheet. Here are the different ways we can use it:
7.5.1 1. Central Tendency (The “Average” View)
This is the most common use. We can calculate the mean, median, or mode of a value within a region.
- Example: Finding the average NDVI (greenness) for each province in Nepal to see which is the “greenest” overall.
7.5.2 2. Extremes (The “Highest and Lowest” View)
We can identify the min or max values within an area.
- Example: Finding the highest elevation point or the maximum temperature recorded within a specific district.
7.5.3 3. Diversity and Spread (The “Variation” View)
We can calculate the standard deviation or variance to see how much the data fluctuates.
- Example: Checking if the forest density is uniform across a province or if there are huge gaps (high variation).
7.5.4 4. Pixel Counting (The “Area” View)
As seen in the previous slide, we can use it to count pixels that fall into specific categories after reclassification.
- Example: Counting how many “Dense Forest” pixels (Class 5) exist in a province versus “Water” pixels (Class 1).
7.5.5 5. Summation (The “Total” View)
We can sum values together.
- Example: Calculating the total amount of rainfall (in mm) that fell over a specific watershed by adding up the values of every pixel within that boundary.
7.5.6 Summary Table: Which Reducer to Choose?
| Goal | GEE Reducer Component |
|---|---|
| Find the Average | ee.Reducer.mean() |
| Find the Peak Value | ee.Reducer.max() |
| Area of Land Cover | ee.Reducer.count() |
| Typical Value | ee.Reducer.median() |
7.6 Zonal Statistics by Province
One of the most powerful features of GEE is the ability to extract statistics for administrative boundaries (Provinces) across time.
7.6.1 Reclassification for Tabular Analysis
We group NDVI into 5 classes (Water, Barren, Shrub, Forest, Dense Forest) to count the pixels in each category.
// 6. Reclassify NDVI
function reclassifyNdvi(image) {
var ndvi = image.select('NDVI');
var reclass = ndvi.lt(0).multiply(1)
.add(ndvi.gte(0).and(ndvi.lt(0.2)).multiply(2))
.add(ndvi.gte(0.2).and(ndvi.lt(0.5)).multiply(3))
.add(ndvi.gte(0.5).and(ndvi.lt(0.7)).multiply(4))
.add(ndvi.gte(0.7).multiply(5))
.int8(); // ← Key fix: Cast to signed 8-bit integer (or use .toInt() / .toUint8())
return image.addBands(reclass.rename('reclass'));
}
var reclassCollection = ndviCollection.map(reclassifyNdvi);7.7 Visualization
var ndviVis = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
var reclassVis = {min: 1, max: 5, palette: ['blue', 'yellow', 'orange', 'green', 'darkgreen']};
years.forEach(function(year) {
var img = reclassCollection.filter(ee.Filter.eq('year', year)).first();
Map.addLayer(img.select('NDVI'), ndviVis, 'NDVI ' + year, false);
Map.addLayer(img.select('reclass'), reclassVis, 'Reclass ' + year, year === 2024);
});7.8 Zonal Histogram:
When we have 5 land-cover classes and 7 provinces, we need to count how many pixels of each class fall into each province. We use a Frequency Histogram.
- The Logic: Imagine each province is a large bucket.
- The Process: GEE looks at every pixel inside a province’s “bucket,” identifies its class (1 to 5), and adds it to a tally.
- The Result: Instead of one single number, each province gets a Dictionary (a list of pairs) that looks like this:
{ "1": 500, "2": 1200, "3": 4500, ... }
Analogy: It’s like counting the different colored marbles (classes) inside several different jars (provinces).
7.9 Flattening: From “Folders” to a “Single List”
After running the analysis for multiple years, your data is structured like a “List of Lists” (each year is a folder containing a list of provinces).
var allHist = ee.FeatureCollection(histCollections).flatten();- The Problem: Most charting tools and CSV exports can’t read “folders.” They need one long, continuous list.
- The Solution:
.flatten()takes all the items out of their year-based folders and puts them into one single, flat FeatureCollection. - Why we do it: This allows the charting tool to see every province for every year in one go, so it can plot the 2008-2024 trends correctly.
7.10 Expanding the Dictionary: Naming the Columns
The raw histogram gives us a computer-friendly “Dictionary,” but we need a Spreadsheet-friendly format for our chart.
return f.set({
'count_1': roundCount(histDict.get('1', 0)),
'count_5': roundCount(histDict.get('5', 0))
});histDict.get('1', 0): This tells the computer: “Look for the count of Class 1. If you don’t find any (e.g., no water in that province), just put 0”.roundCount: Satellite pixels are sometimes calculated as fractions during resampling; we round them to the nearest whole number to represent real “units” of land.
7.11 Zonal Histogram (pixel counts per reclass value per province)
// 8. Zonal Histogram (pixel counts per reclass value per province)
var histCollections = years.map(function(year) {
// Get the single reclassified image for this year
var image = reclassCollection
.filter(ee.Filter.eq('year', year))
.first()
.select('reclass'); // Only the reclass band (integer 1-5)
// Compute frequency histogram per province
var hist = image.reduceRegions({
collection: provinces,
reducer: ee.Reducer.frequencyHistogram(),
scale: 250,
maxPixelsPerRegion: 1e14 // Correct name; increase further if needed (e.g., 1e14)
});
// Add year property
return hist.map(function(f) {
return f.set('year', year);
});
});
7.12 Flatten into one FeatureCollection
var allHist = ee.FeatureCollection(histCollections).flatten();
// Extract counts for all classes (1 to 5)
var allHistExpanded = allHist.map(function(f) {
var histDict = ee.Dictionary(f.get('histogram')); // ← Correct property name!
// Round to nearest integer (most common approach for zonal pixel counts)
var roundCount = function(val) {
return ee.Number(val).round();
};
return f.set({
'count_1': roundCount(histDict.get('1', 0)),
'count_2': roundCount(histDict.get('2', 0)),
'count_3': roundCount(histDict.get('3', 0)),
'count_4': roundCount(histDict.get('4', 0)),
'count_5': roundCount(histDict.get('5', 0)),
// 'total_pixels': roundCount(histDict.values().reduce(ee.Reducer.sum(), [0])) // Optional: sum all for validation
});
});7.13 Visualizing Change
By creating a stacked column chart, we can visualize if “Dense Forest” (Class 5) is increasing or decreasing in specific provinces over the last 16 years.
// Replace 'PROV_NAME' below with the actual property (e.g., 'NAME_1', 'province', etc.)
var chart = ui.Chart.feature.byFeature({
features: allHistExpanded,
xProperty: 'year',
yProperties: ['count_1', 'count_2', 'count_3', 'count_4', 'count_5']
})
.setChartType('ColumnChart')
.setOptions({
title: 'NDVI Reclass Pixel Counts by Province and Year',
hAxis: {title: 'Year', ticks: years.map(function(y) { return {v: y}; })},
vAxis: {title: 'Pixel Count'},
isStacked: true,
legend: {position: 'right'},
series: {
0: {color: 'blue', labelInLegend: 'Class 1 (<0)'},
1: {color: 'yellow', labelInLegend: 'Class 2 (0-0.2)'},
2: {color: 'orange', labelInLegend: 'Class 3 (0.2-0.5)'},
3: {color: 'green', labelInLegend: 'Class 4 (0.5-0.7)'},
4: {color: 'darkgreen', labelInLegend: 'Class 5 (>=0.7)'}
}
});7.14 Exporting Data to Google Drive
To analyze your results in Excel or R, you can export the processed table as a CSV file directly to your Google Drive.
// Create clean table WITHOUT geometry
var cleanForExport = allHistExpanded
.map(function(feature) {
// Keep only the properties we want, set geometry to null
return ee.Feature(null).copyProperties(feature, [
'year',
'count_1', 'count_2', 'count_3', 'count_4', 'count_5',
'FIRST_STAT' // Your province name property
]);
});
print(cleanForExport);
// Export the clean, geometry-free table
Export.table.toDrive({
collection: cleanForExport,
description: 'Nepal_NDVI_Histogram_SuperClean',
folder: 'GEE_Exports',
fileNamePrefix: 'NDVI_Hist_NoGeo',
fileFormat: 'CSV'
});7.14.1 Steps to Complete the Export:
- Run the Script: Click the Run button in the GEE Editor.
- Open Tasks: Go to the Tasks tab in the right-hand panel.
- Initiate Task: Click the Run button next to
Annual_LandCover_Area_Nepal. - Confirm: A dialog will appear; click Run again to start the cloud processing.
7.15 Practical Exercise
- Modify the
yearslist to include every year from 2008 to 2024. - Change the
startandenddates to October to compare the “Pre-Monsoon” vs “Post-Monsoon” vegetation health.
By comparing these sensors, you can learn how to choose the right tool for different environmental challenges in Nepal, such as glacial monitoring, urban expansion, or forest fire detection.
8 Choosing the Right Satellite for Nepal
Different environmental problems require different “eyes” in the sky. While high resolution is great for detail, high frequency is better for monitoring fast-changing events like floods.
8.1 Sensor Comparison Table
| Feature | MODIS | Landsat (8/9) | Sentinel-2 |
|---|---|---|---|
| Spatial Resolution | 250m to 1km (NASA LP DAAC 2021) | 30m | 10m to 20m |
| Temporal Resolution | Daily | 8 to 16 days | 5 days |
| Record Starts | 2000 | 1972 (Landsat 1) | 2015 (Drusch et al. 2012) |
| Best For | Regional / National | Historical Trends | Local / Detailed |
8.2 Application Guide for Nepal
| Environmental Problem | Recommended Sensor | Reason |
|---|---|---|
| Glacial Lake Outburst (GLOF) | Sentinel-2 | High resolution (10m) is needed to see small mountain lakes(Drusch et al. 2012; Ishtiaque, Shrestha, and Chhetri 2017). |
| Forest Fire Monitoring | MODIS | Daily revisits allow us to see fires almost in real-time.(Giglio, Schroeder, and Justice 2016; Davies et al. 2009). |
| Urban Growth (Kathmandu) | Landsat | We need the long history (since the 80s) to see the growth of the city. (Rimal et al. 2017; Ishtiaque, Shrestha, and Chhetri 2017). |
| Agricultural Monitoring | Sentinel-2 | 10m resolution is small enough for the terrace farming found in Nepal’s hills(Drusch et al. 2012).. |
| Snow Cover Mapping | MODIS | Clouds often block the view; daily revisits increase the chance of a clear day. (Hall et al. 2002) |
8.3 Quick GEE Reference
| Task | SNAP / QGIS Step | GEE Snippet | Citation Source |
|---|---|---|---|
| Resampling | r.resample |
Handled automatically by Scale. | (Gorelick et al. 2017) |
| NDVI | Band Maths | .normalizedDifference(['NIR', 'Red']) |
(Rouse et al. 1974) |
| Zonal Stats | Zonal Statistics |
.reduceRegions() |
(Google Earth Engine, n.d.) or (Gorelick et al. 2017) |
| Clipping | Clip by Mask |
.clip() |
(Google Earth Engine, n.d.) |