2. Time-series analysis

Previous1. Introduction to Earth Engine Next3. Image classification

Last updated 26 days ago

Was this helpful?

2. Time-series analysis

Objective

The objective of today's session is to provide you with the critical skills needed to analyze image collections. In optical image processing, cloud masking is an important first step. Based on this pre-processing operation, we can create data statistics, visualizations and analysis.

Cloud masking

Reliable identification of clouds is necessary for many types of optical remote sensing image analysis.

Landsat

Landsat, a joint program of the USGS and NASA, has been observing the Earth continuously from 1972 through the present day. Today the Landsat satellites image the entire Earth's surface at a 30-meter resolution about once every two weeks, including multispectral and thermal data.

Simple cloud score

For scoring Landsat pixels by their relative cloudiness, Earth Engine provides a rudimentary cloud scoring algorithm in the ee.Algorithms.Landsat.simpleCloudScore() method. (For details on the implementation, see ). The following example uses the cloud scoring algorithm to mask clouds in a Landsat 8 image:

///////////////////////////////////////////////////
// 1. LOAD AND FILTER LANDSAT 8 COLLECTION

var L8_toa_ic = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA")
  .filterDate('2018-06-01', '2019-06-01')
  .filterBounds(ee.Geometry.Point(23.6, 37.9))  // Athens, Greece
  .filter('CLOUD_COVER < 70')
  .filter('CLOUD_COVER > 20');  // Keep moderately cloudy images

// Select one image from the collection for analysis
var image = ee.Image(L8_toa_ic.first());
print('Selected image:', image);


///////////////////////////////////////////////////
// 2. VISUALIZE ORIGINAL IMAGE (with clouds)

var vizParams = {
  bands: ['B7', 'B5', 'B3'],  // False-color (SWIR-NIR-Green)
  min: 0.05,
  max: 0.3
};

Map.setOptions('SATELLITE');
Map.centerObject(image, 9);
Map.addLayer(image, vizParams, 'Original (Cloudy) Image');


///////////////////////////////////////////////////
// 3. APPLY CLOUD MASK (Simple Cloud Score)

var scored = ee.Algorithms.Landsat.simpleCloudScore(image);
var cloudMask = scored.select('cloud').lte(20);  // Keep only low-cloud pixels

var imageMasked = image.updateMask(cloudMask);
Map.addLayer(imageMasked, vizParams, 'Masked Image (clouds removed)');

Map.addLayer(scored.select('cloud'), {min: 0, max: 100, palette: ['white', 'blue']}, 
'Cloud Score (0–100)', false);

We can use this algorithm to wrap it into a function and map it over a Landsat image collection.

///////////////////////////////////////////////////
// 1. LOAD LANDSAT 8 TOA IMAGE COLLECTION

var L8_toa_ic = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA")
  .filterDate('2018-06-01', '2019-06-01')
  .filterBounds(ee.Geometry.Point(23.6, 37.9));  // Athens, Greece
print('Initial image collection:', L8_toa_ic);


///////////////////////////////////////////////////
// 2. CLOUD MASKING FUNCTION

var maskClouds = function(image) {
  // Apply simple cloud score algorithm
  var scored = ee.Algorithms.Landsat.simpleCloudScore(image);
  // Create cloud mask (threshold: cloud score < 20)
  var mask = scored.select('cloud').lt(20);
  // Apply cloud mask and return masked image
  return image.updateMask(mask);
};


///////////////////////////////////////////////////
// 3. APPLY MASKING TO COLLECTION

var L8_toa_ic_cloudfree = L8_toa_ic.map(maskClouds);
print('Cloud-masked image collection:', L8_toa_ic_cloudfree);


///////////////////////////////////////////////////
// 4. CREATE A MEDIAN COMPOSITE

var L8_toa_median = L8_toa_ic_cloudfree.median();

var visParams = {
  bands: ['B7', 'B5', 'B3'],  // SWIR-NIR-Green (False color)
  min: 0.05,
  max: 0.4
};


///////////////////////////////////////////////////
// 5. VISUALIZATION

Map.setOptions('SATELLITE');
Map.centerObject(ee.Geometry.Point(23.6, 37.9), 9);  // Center over Athens
Map.addLayer(L8_toa_median, visParams, 'Cloud-free Median Composite');

Fmask

// Landsat 8 Collection 2, Level 2 – QA Masking & Scaling

///////////////////////////////////////////////////
// 1. FUNCTION TO MASK CLOUDS AND APPLY SCALING

function maskL8sr(image) {
  // Extract QA_PIXEL bitmask:
  // Bit 0: Fill
  // Bit 1: Dilated Cloud
  // Bit 2: Cirrus
  // Bit 3: Cloud
  // Bit 4: Cloud Shadow
  var qaMask = image.select('QA_PIXEL')
    .bitwiseAnd(parseInt('11111', 2))  // Mask bits 0–4
    .eq(0);  // Keep only pixels where all bits are zero (clear)

  // Exclude saturated pixels
  var satMask = image.select('QA_RADSAT').eq(0);

  // Apply scale factors (from USGS documentation)
  var optical = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermal = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Return masked, scaled image
  return image
    .addBands(optical, null, true)
    .addBands(thermal, null, true)
    .updateMask(qaMask)
    .updateMask(satMask);
}


///////////////////////////////////////////////////
// 2. LOAD AND PROCESS IMAGE COLLECTION

var L8_sr = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterDate('2020-01-01', '2021-01-01')  // Select year
  .map(maskL8sr);  // Apply cloud & saturation mask


///////////////////////////////////////////////////
// 3. COMPOSITE & VISUALIZATION

var composite = L8_sr.median();  // Cloud-free median composite

var visParams = {
  bands: ['SR_B7', 'SR_B5', 'SR_B4'],  // False-color infrared
  min: 0,
  max: 0.3
};

Map.setOptions('SATELLITE');
Map.setCenter(114.12466, -21.93527, 7);  // Exmouth, Australia
Map.addLayer(composite, visParams, 'Cloud-free L8 SR Composite');

MODIS

// MODIS Cloud Masking Example
// Purpose: Compute the frequency of clear-sky observations from MODIS QA flags.

///////////////////////////////////////////////////
// 1. DEFINE MASKING FUNCTIONS

// Mask pixels without any observations
function maskEmptyPixels(img) {
  return img.updateMask(img.select('num_observations_1km').gt(0));
}

// Mask cloudy pixels based on the internal cloud algorithm flag (bit 10)
function maskClouds(img) {
  var QA = img.select('state_1km');
  var cloudBitMask = 1 << 10;  // Bit 10 = internal cloud flag
  return img.updateMask(QA.bitwiseAnd(cloudBitMask).eq(0));
}

///////////////////////////////////////////////////
// 2. LOAD COLLECTION & APPLY MASKING

// Define time period
var startDate = '2010-01-01';
var endDate = '2011-01-01';

// Load MODIS surface reflectance data
var modis = ee.ImageCollection('MODIS/061/MOD09GA')
  .filterDate(startDate, endDate)
  .map(maskEmptyPixels);

// Count number of total valid observations
var totalObs = modis
  .select('num_observations_1km')
  .count();

// Apply cloud mask
var modisClear = modis.map(maskClouds);

// Count number of cloud-free observations
var clearObs = modisClear
  .select('num_observations_1km')
  .count()
  .unmask(0);  // Fill no-data with 0 for ratio calculation

// Compute clear observation ratio
var clearRatio = clearObs.toFloat().divide(totalObs)
  .rename('clear_ratio');

///////////////////////////////////////////////////
// 3. VISUALIZATION

// Center map over Europe (optional)
Map.setCenter(10, 50, 4);

// Load Gennadii Donchyts' palettes package from the Earth Engine community
var palettes = require('users/gena/packages:palettes');
var devon = palettes.crameri.devon[25].reverse();   

// RGB Composite of median surface reflectance
Map.addLayer(modisClear.median(),
  {bands: ['sur_refl_b01', 'sur_refl_b04', 'sur_refl_b03'], gain: 0.07, gamma: 1.4},
  'MODIS Median (Cloud Masked)', false);

// Total number of observations
Map.addLayer(totalObs, {min: 300, max: 365}, 'Total Observations', false);

// Number of clear observations
Map.addLayer(clearObs, {min: 0, max: 365}, 'Clear Observations', false);

// Ratio of clear to total observations
Map.addLayer(clearRatio, {min: 0, max: 1, palette: devon}, 'Clear-Sky Frequency');

Playtime

Task: Find the sunniest place on Earth [make us of ee.Image.pixelLonLat() and reduceRegion()].

Sentinel-2

Sentinel-2 is a wide-swath, high-resolution, multi-spectral imaging mission supporting Copernicus Land Monitoring studies, including the monitoring of vegetation, soil and water cover, as well as observation of inland waterways and coastal areas.

Cloud probability

////////////////////////////////////////////////////////////
// 1. LOAD INPUT DATA

// Sentinel-2 Surface Reflectance (Level-2A)
var s2SR = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');

// Sentinel-2 Cloud Probability dataset
var s2CloudProb = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY');

////////////////////////////////////////////////////////////
// 2. SETTINGS

// Date range and cloud threshold
var startDate = '2018-01-01';
var endDate = '2020-01-01';
var maxCloudProb = 40;

// Area of interest: a rectangle over part of Colombia
var region = ee.Geometry.Rectangle([-76.5, 2.0, -74, 4.0]);

// Visualization: Shortwave IR, NIR, Red (good for vegetation & moisture)
var rgbVis = {
  min: 0,
  max: [3000, 4000, 2500],
  bands: ['B12', 'B8', 'B4']
};

////////////////////////////////////////////////////////////
// 3. FUNCTIONS

// Cloud masking function using the 'probability' band from cloud probability dataset
function maskClouds(image) {
  var cloudProb = image.select('probability');
  var mask = cloudProb.lt(maxCloudProb);
  return image.updateMask(mask);
}

// Edge masking to avoid artifacts from invalid pixels (common in S2 data)
function maskSceneEdges(image) {
  var mask = image.select('B8A').mask().updateMask(image.select('B9').mask());
  return image.updateMask(mask);
}

////////////////////////////////////////////////////////////
// 4. FILTER AND PREPARE COLLECTIONS

// Common filter: by region and date
var filter = ee.Filter.and(
  ee.Filter.bounds(region),
  ee.Filter.date(startDate, endDate)
);

// Filter surface reflectance and cloud probability collections
var s2SR_filtered = s2SR.filter(filter).map(maskSceneEdges);
var s2CloudProb_filtered = s2CloudProb.filter(filter);

////////////////////////////////////////////////////////////
// 5. LINK CLOUD PROBABILITY + MASK CLOUDS

// Link corresponding cloud probability image to each S2 SR image
var s2_cloudMasked = s2SR_filtered
  .linkCollection(s2CloudProb_filtered, ['probability'])
  .map(maskClouds)
  .median();  // Reduce time series to median composite

print('Cloud-masked composite image', s2_cloudMasked);

////////////////////////////////////////////////////////////
// 6. VISUALIZATION

Map.setCenter(-75.4, 2.6, 12);  // Center on the area of interest
Map.addLayer(s2_cloudMasked, rgbVis, 'Sentinel-2 SR (Masked, ' + maxCloudProb + '% CloudThreshold)');

Cloud Score+

////////////////////////////////////////////////////////////
// 1. DATASETS
// Sentinel-2 Surface Reflectance
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');
// Cloud Score+  
var s2Csp = ee.ImageCollection('GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED');  

////////////////////////////////////////////////////////////
// 2. USER PARAMETERS

var startDate = ee.Date('2018-01-01');
var endDate = ee.Date('2020-01-01');
var minClearProbability = 0.6;

var region = ee.Geometry.Rectangle([-76.5, 2.0, -74, 4.0], null, false);  // Colombia

////////////////////////////////////////////////////////////
// 3. FUNCTIONS

// Cloud masking using Cloud Score+ cumulative density function (cs_cdf)
function maskClouds(image) {
  var clearSkyProb = image.select('cs_cdf');
  return image.updateMask(clearSkyProb.gte(minClearProbability));
}

// Edge masking: mask pixels based on invalid 20m and 60m band masks
function maskEdges(image) {
  return image.updateMask(
    image.select('B8A').mask().updateMask(image.select('B9').mask())
  );
}

////////////////////////////////////////////////////////////
// 4. PROCESSING

// Filter both collections by date and region
var filter = ee.Filter.and(
  ee.Filter.bounds(region),
  ee.Filter.date(startDate, endDate)
);

var s2SrFiltered = s2Sr.filter(filter).map(maskEdges);
var s2CspFiltered = s2Csp.filter(filter);

// Link S2 surface reflectance with Cloud Score+ via cs_cdf, mask, and reduce
var s2Masked = s2SrFiltered
  .linkCollection(s2CspFiltered, ['cs_cdf'])  // Join CSP data
  .map(maskClouds)                            // Apply cloud mask
  .median();                                  // Composite

print('Cloud-masked Sentinel-2 composite', s2Masked);

////////////////////////////////////////////////////////////
// 5. VISUALIZATION

var visParams = {
  bands: ['B12', 'B8', 'B4'],  // SWIR - NIR - Red (False Color Composite)
  min: 0,
  max: [3000, 4000, 2500]
};

Map.setCenter(-75.4, 2.6, 11);  // Zoom to area of interest
Map.addLayer(s2Masked, visParams, 
'S2 Cloud-Masked (' + (minClearProbability * 100) + '% CloudThreshold)');

Cloud Score+ not only masks clouds, it also helps detect cloud shadows. The script below additionally detects and masks cloud shadows using clouds and solar shading geometry.

Cloud displacement

The Cloud Displacement Index (CDI) makes use of the three highly correlated near infrared bands of Sentinel-2 that are observed with different view angles. Hence, elevated objects like clouds are observed under a parallax and can be reliably separated from bright ground objects.

Data availability

Before you start any study, you should verify what kind of data is available for your purposes. Here is one script for Landsat 4, 5, 7, 8 and 9 data that indicates for a given area and time span the average per-pixel revisit time and average number of observation per satellite:

The same principle can be applied to other satellites, such as Sentinel-1 or Sentinel-2.

Playtime

Task: Navigate to a location of interest and find out more about Landsat and Sentinel data availability.

Forest monitoring

Forests around the world are in crisis. Rapidly expanding global human footprint and climate change have caused extensive deforestation, and have left remnant forests fragmented and significantly altered in condition. Yet, these remaining forests harbour a dazzling array of unique ecosystems and an overwhelming majority of our planet’s current terrestrial biodiversity, thanks to in-situ conservation efforts including those by local and indigenous communities. Valuing these forests and ensuring their continued conservation depends on our understanding their long term dynamics and responses to dominant anthropogenic and climatic drivers. These scripts illustrates the use of Earth Engine to investigate forest vegetation condition over time.

Vegetation indices

Remotely sensed indices such as enhanced vegetation index (EVI), normalized difference vegetation index (NDVI), normalized burn ratio (NBR) and normalized difference water index (NDWI) are widely used to estimate vegetation status from satellite imagery. EVI and NDVI estimate vegetation chlorophyll content while the NBR estimates the severity of fires and the NDWI estimates vegetation moisture content. These indices can be derived from Landsat and other satellites with similar spectral settings. The following script illustrates the combined use of Landsat 5, 7 and 8 for time-series analysis of these indices in forest areas.

Task: In 2003 one of the largest forest fires in recent history (in Switzerland) occurred close to Leuk (7.651, 46.333) in Valais. Investigate based on available Landsat data how many years it took for the vegetation to recover to the pre-event level.

Trend analysis

To infer vegetation conditions, we infer a linear trend at each pixel by calculating its Sen's slope of maximum summer EVI with time. For this analysis we use the MODIS 250m/pixel 16-day composite vegetation indices dataset and build an image collection with an image for each year from 2000 to 2020. Each of these images is calculated to be the mean EVI in the summer months of its corresponding year. This value is our measure of the status of the vegetation for each year. Furthermore, we add the year as a band, in preparation for linear trend analysis. As a result, we can infer the pixel-wise vegetation greening or browning based on the sign of the slope value. Finally, we visualize results in forest areas and calculate summarizing stats for our area of interest.

Task: Check out the vegetation trends in you area of interest.

Climatologies

Information on climatic conditions is essential to many applications in environmental and ecological sciences. Working with remote sensing images is a global endeavor. Using climatological information helps to better understand local conditions in remote areas. The attached script provides information on precipitation, cloud cover, air temperatures, snow cover and vegetation cover.

Task: Derive climatological information of your favorite area of interest.

Thumbnails

The getVideoThumbURL() function generates an animation from all images in an ImageCollection where each image represents a frame. The general workflow for producing an animation is as follows:

Define a Geometry whose bounds determine the regional extent of the animation.
Define an ImageCollection.
Consider image visualization: either map an image visualization function over the collection or add image visualization arguments to the set of animation arguments.
Define animation arguments and call the getVideoThumbURL method.

The result of getVideoThumbURL is a URL. Print the URL to the console and click it to start Earth Engine servers generating the animation on-the-fly in a new browser tab. Alternatively, view the animation in the Code Editor console by calling the ui.Thumbnail function on the collection and its corresponding animation arguments. Upon rendering, the animation is available for downloading by right clicking on it and selecting appropriate options from its context menu.

The following example illustrates generating an animation depicting global temperatures over the course of 24 hours. Note that this example includes visualization arguments along with animation arguments, as opposed to first mapping a visualization function over the ImageCollection. Upon running this script, an animation similar to the Figure below should appear in the Code Editor console.

//////////////////////////////////////////////////////////////////////////////////////
// GFS Predicted Temperature Animation (2m Above Ground)
// NOAA Global Forecast System (GFS) hourly data visualization
//////////////////////////////////////////////////////////////////////////////////////

// 1. DEFINE AREA OF INTEREST (AOI)

// Global non-polar extent to avoid projection issues near poles
var aoi = ee.Geometry.Polygon(
  [[[-179.0, 78.0], [-179.0, -58.0], [179.0, -58.0], [179.0, 78.0]]],
  null, false
);

//////////////////////////////////////////////////////////////////////////////////////
// 2. LOAD GFS TEMPERATURE DATA

// Load hourly NOAA GFS 0.25° forecast data
var gfsTemp = ee.ImageCollection('NOAA/GFS0P25')
  .filterDate('2021-06-16', '2021-06-18') // Select forecast for a 2-day period
  .limit(48)                              // Limit to the first 48 hourly predictions
  .select('temperature_2m_above_ground'); // Select surface air temperature at 2 meters

//////////////////////////////////////////////////////////////////////////////////////
// 3. VISUALIZATION SETTINGS

// Define rendering parameters for temperature
var visParams = {
  min: -20,                // Lower temperature bound (°C)
  max: 35,                 // Upper temperature bound (°C)
  palette: [               // Gradient palette from cold to warm
    '#042333', '#2c3395', '#744992', '#b15f82',
    '#eb7958', '#fbb43d', '#e8fa5b'
  ]
};

// Define video thumbnail export parameters
var videoParams = {
  dimensions: 768,         // Output width in pixels
  region: aoi,             // Clipping region for the video
  framesPerSecond: 7,      // Playback speed
  crs: 'EPSG:3857',        // Projection for rendering
  min: visParams.min,
  max: visParams.max,
  palette: visParams.palette
};

//////////////////////////////////////////////////////////////////////////////////////
// 4. DISPLAY ANIMATION

// Show the temperature animation in the console
print('GFS Temperature Forecast Animation (°C)', ui.Thumbnail(gfsTemp, videoParams));

// Optionally, generate a direct link to the animation video
print('Downloadable Animation URL:', gfsTemp.getVideoThumbURL(videoParams));

Playtime

Task: Use the "ECMWF/ERA5_LAND/MONTHLY_AGGR" dataset to animate annual temperature variations.

The following sections describe how to use clipping and layer compositing to enhance visualizations by adding polygon borders and comparing images within a collection.

Multiple images can be overlaid using the blend Image method where overlapping pixels from two images are blended based on their masks (opacity). Multiple images can be overlaid using the blend Image method where overlapping pixels from two images are blended based on their masks (opacity). Vector data (Features) are drawn to images by applying the paint method. Features can be painted to an existing image, but the better practice is to paint them to a blank image, style it, and then blend the result with other styled image layers. You can also blend image data with a hillshade base layer to indicate terrain and give the visualization some depth.

... and here comes another blended animation of cumulative precipitation in the Indus catchment from 14.07.2022 to 30.08.2022.

... and here a "cool" example from the Panmah Glacier in the Karakorum featuring cloud-filtered annual Sentinel-2 composites with minimum snow cover between 2017 and 2023.

Assignment

Select one of the following topics dealing with time series analysis or define your own topic according to your interests. Your results should be summarized in a chart that features time on the x-axis, accompanied by a brief description. You will find more details in the assignment. This exercise is meant to apply and deepen the knowledge and skills you have acquired so far. You can work in teams or individually.

Previous1. Introduction to Earth Engine Next3. Image classification

Last updated 26 days ago

Was this helpful?