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 this Code Editor sample script). The following example uses the cloud scoring algorithm to mask clouds in a Landsat 8 image:

Open in Code Editor

// Load an ImageCollection.
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))
  .filter('CLOUD_COVER < 70')
  .filter('CLOUD_COVER > 20');
    
var image = ee.Image(L8_toa_ic.toList(5).get(0));    

// Define the visualization parameters.
var vizParams = {
  bands: ['B7', 'B5', 'B3'],
  min: 0.05,
  max: 0.3,
};

// Center the map and display the image.
Map.centerObject(image, 9); 
Map.addLayer(image, vizParams, 'cloudy L8 image', true);

// Add a cloud score band.  It is automatically called 'cloud'.
var scored = ee.Algorithms.Landsat.simpleCloudScore(image);

// Create a mask from the cloud score and combine it with the image mask.
var mask = scored.select(['cloud']).lte(20);

// Apply the mask to the image and display the result.
var image_masked = image.updateMask(mask);
Map.addLayer(image_masked, vizParams, 'image masked');

If you run this example in the Code Editor, try toggling the visibility of the TOA layers to compare the difference between the masked and unmasked imagery. Observe that the input to simpleCloudScore() is a single Landsat TOA scene. Also note that simpleCloudScore() adds a band called ‘cloud’ to the input image. The cloud band contains the cloud score from 0 (not cloudy) to 100 (most cloudy). The previous example uses an arbitrary threshold (20) on the cloud score to mask cloudy pixels.

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

Open in Code Editor

// Load an 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))

// This function masks clouds in Landsat 8 imagery.
var maskClouds = function(image) {
  var scored = ee.Algorithms.Landsat.simpleCloudScore(image);
  return image.updateMask(scored.select(['cloud']).lt(20));
};

// Map the cloud masking function over the collection.
var L8_toa_ic_cloudfree = L8_toa_ic.map(maskClouds);

// Compute a cloud-filtered median image and display.
var L8_toa_ic_cloudfree_median = L8_toa_ic_cloudfree.median();
Map.addLayer(L8_toa_ic_cloudfree_median, 
    {bands: ['B7', 'B5', 'B3'], max: 0.4}, 
    'L8_toa_ic_cloudfree_median');

Fmask

The ImageCollection Landsat 8 Surface Reflectance Tier 1 have been atmospherically corrected using LaSRC and includes a cloud, shadow, water and snow mask produced using CFMASK, as well as a per-pixel saturation mask.

Open in Code Editor

// This example demonstrates the use of the Landsat 8 Collection 2, Level 2
// QA_PIXEL band (CFMask) to mask unwanted pixels.

function maskL8sr(image) {
  // 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)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones and apply the masks.
  return image.addBands(opticalBands, null, true)
      .addBands(thermalBands, null, true)
      .updateMask(qaMask)
      .updateMask(saturationMask);
}

// Map the function over one year of data.
var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                     .filterDate('2020-01-01', '2021-01-01')
                     .map(maskL8sr);

var composite = collection.median();

// Display the results.
Map.setCenter(114.12466, -21.93527, 7);  // Exmouth, Australia
Map.addLayer(composite, {bands: ['SR_B7', 'SR_B5', 'SR_B4'], min: 0, max: 0.3});

MODIS

The MODIS (Moderate Resolution Imaging Spectroradiometer) program offers a large collection of Earth observation products. The foundation for many of these derivatives is the MODIS Surface Reflectance product MOD09GA. The MODIS Surface Reflectance products provide an estimate of the surface spectral reflectance as it would be measured at ground level in the absence of atmospheric scattering or absorption. Low-level data are corrected for atmospheric gases and aerosols. MOD09GA version 6 provides bands 1-7 in a daily gridded L2G product in the sinusoidal projection, including 500m reflectance values and 1km observation and geolocation statistics.

Open in Code Editor

// Modis Cloud Masking example.
// Calculate how frequently a location is labeled as clear (i.e. non-cloudy)
// according to the "internal cloud algorithm flag" of the MODIS "state 1km"
// QA band.

// A function to mask out pixels that did not have observations.
var maskEmptyPixels = function(image) {
  // Find pixels that had observations.
  var withObs = image.select('num_observations_1km').gt(0)
  return image.updateMask(withObs)
}

// A function to mask out cloudy pixels.
var maskClouds = function(image) {
  // Select the QA band.
  var QA = image.select('state_1km')
  // Make a mask to get bit 10, the internal_cloud_algorithm_flag bit.
  var bitMask = 1 << 10;
  // Return an image masking out cloudy areas.
  return image.updateMask(QA.bitwiseAnd(bitMask).eq(0))
}

// Start with an image collection for a 1 month period.
// and mask out areas that were not observed.
var collection = ee.ImageCollection('MODIS/061/MOD09GA')
        .filterDate('2010-01-01', '2011-01-01')
        .map(maskEmptyPixels)

// Get the total number of potential observations for the time interval.
var totalObsCount = collection
        .select('num_observations_1km')
        .count()

// Map the cloud masking function over the collection.
var collectionCloudMasked = collection.map(maskClouds)

// Get the total number of observations for non-cloudy pixels for the time
// interval.  The result is unmasked to set to unity so that all locations
// have counts, and the ratios later computed have values everywhere.
var clearObsCount = collectionCloudMasked
        .select('num_observations_1km')
        .count()
        .unmask(0)

Map.addLayer(collectionCloudMasked.median(),
    {bands: ['sur_refl_b01', 'sur_refl_b04', 'sur_refl_b03'], gain: 0.07,  gamma: 1.4}, 
    'median of masked collection')
    
Map.addLayer(totalObsCount,
    {min: 300, max: 365},
    'count of total observations', false)

Map.addLayer(clearObsCount,
    {min: 0, max: 365},
    'count of clear observations', false)
    
Map.addLayer(clearObsCount.toFloat().divide(totalObsCount),
    {min: 0, max: 1},
    'ratio of clear to total observations')

Playtime

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

Use the Code Editor

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

The S2 cloud probability is created with the sentinel2-cloud-detector library (using LightGBM). All bands are upsampled using bilinear interpolation to 10m resolution before the gradient boost base algorithm is applied. The resulting 0..1 floating point probability is scaled to 0..100 and stored as a UINT8. Areas missing any or all of the bands are masked out. Higher values are more likely to be clouds or highly reflective surfaces (e.g. roof tops or snow).

Open in Code Editor

////////////////////////////////////////////////////////////
// Data
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');
var s2Clouds = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY');

////////////////////////////////////////////////////////////
// Definitions
var start_date = ee.Date('2018-01-01');
var end_date = ee.Date('2020-01-01');
var max_cloud_probability = 40;
var region =
    ee.Geometry.Rectangle({coords: [-76.5, 2.0, -74, 4.0], geodesic: false});

////////////////////////////////////////////////////////////
// Functions
function maskClouds(img) {
  var clouds = img.select('probability');
  var isNotCloud = clouds.lt(max_cloud_probability);
  return img.updateMask(isNotCloud);
}

// The masks for the 10m bands sometimes do not exclude bad data at
// scene edges, so we apply masks from the 20m and 60m bands as well.
// Example asset that needs this operation:
// COPERNICUS/S2_CLOUD_PROBABILITY/20190301T000239_20190301T000238_T55GDP
function maskEdges(s2_img) {
  return s2_img.updateMask(
      s2_img.select('B8A').mask().updateMask(s2_img.select('B9').mask()));
}

////////////////////////////////////////////////////////////
// Processing

// Filter input collections by desired data range and region.
var criteria = ee.Filter.and(
    ee.Filter.bounds(region), 
    ee.Filter.date(start_date, end_date));

// apply the filter
s2Sr = s2Sr.filter(criteria).map(maskEdges);
s2Clouds = s2Clouds.filter(criteria);

// link collections, mask clouds and apply reducer
var s2CloudMasked = s2Sr.linkCollection(s2Clouds, ['probability'])
                        .map(maskClouds)
                        .median();
print(s2CloudMasked)
////////////////////////////////////////////////////////////
// Display
var rgbVis = {min: 0, max: [3000, 4000, 2500], bands: ['B12', 'B8', 'B4']};

Map.addLayer(
    s2CloudMasked, rgbVis, 'S2 SR masked at ' + max_cloud_probability + '%',
    true);
Map.setCenter(-75.4, 2.6, 12);

Cloud Score+

Cloud Score+ is a quality assessment (QA) processor for medium-to-high resolution optical satellite imagery. It includes two QA bands, "cs" and "cs_cdf", that both grade the usability of individual pixels with respect to surface visibility on a continuous scale between 0 and 1, where 0 represents "not clear" (occluded), while 1 represents "clear" (unoccluded) observations. For more information about the Cloud Score+ dataset and modelling approach, see this Medium post.

Open in Code Editor

////////////////////////////////////////////////////////////
// Data
var s2Sr = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED');
var S2_csp = ee.ImageCollection('GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED')

////////////////////////////////////////////////////////////
// Definitions
var start_date = ee.Date('2018-01-01');
var end_date = ee.Date('2020-01-01');
var min_clear_probability = 0.6;
var region =
    ee.Geometry.Rectangle({coords: [-76.5, 2.0, -74, 4.0], geodesic: false});

////////////////////////////////////////////////////////////
// Functions
function maskClouds(img) {
  var clouds_csp = img.select('cs_cdf');
  // the minimum clear sky probability threshold is set at 0.6
  var isNotCloud = clouds_csp.gte(min_clear_probability);
  return img.updateMask(isNotCloud);
}

// The masks for the 10m bands sometimes do not exclude bad data at
// scene edges, so we apply masks from the 20m and 60m bands as well.
// Example asset that needs this operation:
// COPERNICUS/S2_CLOUD_PROBABILITY/20190301T000239_20190301T000238_T55GDP
function maskEdges(s2_img) {
  return s2_img.updateMask(
      s2_img.select('B8A').mask().updateMask(s2_img.select('B9').mask()));
}

////////////////////////////////////////////////////////////
// Processing

// Filter input collections by desired data range and region.
var criteria = ee.Filter.and(
    ee.Filter.bounds(region), 
    ee.Filter.date(start_date, end_date));

// apply the filter
s2Sr = s2Sr.filter(criteria).map(maskEdges);
S2_csp = S2_csp.filter(criteria);


// link collections, mask clouds and apply reducer
var s2CloudMasked = s2Sr.linkCollection(S2_csp, ['cs_cdf'])
                        .map(maskClouds)
                        .median();
print(s2CloudMasked)
////////////////////////////////////////////////////////////
// Display
var rgbVis = {min: 0, max: [3000, 4000, 2500], bands: ['B12', 'B8', 'B4']};

Map.addLayer(
    s2CloudMasked, rgbVis, 'S2 SR masked at ' + min_clear_probability*100 + '%',
    true);
Map.setCenter(-75.4, 2.6, 12);

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.

Open in Code Editor (Sentinel-2 cloud and cloud shadow masking using Cloud Score+)

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.

Open in Code Editor (Sentinel-2 cloud displacement index)

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:

Open Code in Editor (Landsat data availability in Space)

Open Code in Editor (Landsat data availability in Time)

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

Open Code in Editor (Sentinel-1 data availability)

Open Code in Editor (Sentinel-2 data availability)

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.

Open in Code Editor (Landsat vegetation indices time series)

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.

Open in Code Editor (MODIS EVI trends using Sens Slope)

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.

Open in Code Editor (Script to derive climatologies)

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:

  1. Define a Geometry whose bounds determine the regional extent of the animation.

  2. Define an ImageCollection.

  3. Consider image visualization: either map an image visualization function over the collection or add image visualization arguments to the set of animation arguments.

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

Open in Code Editor

// Define an area of interest geometry with a global non-polar extent.
var aoi = ee.Geometry.Polygon(
  [[[-179.0, 78.0], [-179.0, -58.0], [179.0, -58.0], [179.0, 78.0]]], null,
  false);

// Import hourly predicted temperature image collection. 
// Note that predictions extend for 384 hours; 
// limit the collection to the first 48 hours. 
var tempCol = ee.ImageCollection('NOAA/GFS0P25')
  .filterDate('2021-06-16', '2021-06-18')
  .limit(48)
  .select('temperature_2m_above_ground');

// Define arguments for animation function parameters.
var videoArgs = {
  dimensions: 768,
  region: aoi,
  framesPerSecond: 7,
  crs: 'EPSG:3857',
  min: -20.0,
  max: 35.0,
  palette: ['042333', '2c3395', '744992', 'b15f82', 'eb7958', 'fbb43d', 'e8fa5b']
};

// Print the animation to the console as a ui.Thumbnail using the above defined
// arguments. Note that ui.Thumbnail produces an animation when the first input
// is an ee.ImageCollection instead of an ee.Image.
print(ui.Thumbnail(tempCol, videoArgs));

// Alternatively, print a URL that will produce the animation when accessed.
print(tempCol.getVideoThumbURL(videoArgs));

Playtime

Task: Use the "ECMWF/ERA5/MONTHLY" dataset to animate annual temperature variations.

Use the Code Editor

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.

Open in Code Editor (Thumbnails - Temperature Overlay)

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

Open in Code Editor (Thumbnails - Precipitation Overlay)

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

Open in Code Editor

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 one plot/chart that features time on the x-axis. This exercise is meant to apply and deepen the knowledge and skills you have acquired so far. You can work in groups of two or three, or individually.

Previous1. Introduction to Earth EngineNext3. Image classification

Last updated