End-to-End Google Earth Engine (Full Course Material)
A hands-on introduction to applied remote sensing using Google Earth Engine.
Ujaval Gandhi
- Introduction
- Setting up the Environment
- Module 1: Earth Engine Basics
- Module 2: Earth Engine Intermediate
- Module 3: Supervised Classification
- Module 4: Change Detection
- Module 5: Earth Engine Apps
- Module 6: Google Earth Engine Python API
- Supplement
- Guided Projects
- Learning Resources
- Data Credits
- License
- Citing and Referencing
Introduction
Google Earth Engine is a cloud-based platform that enables large-scale processing of satellite imagery to detect changes, map trends, and quantify differences on the Earth’s surface. This course covers the full range of topics in Earth Engine to give the participants practical skills to master the platform and implement their remote sensing projects.
Setting up the Environment
Sign-up for Google Earth Engine
If you already have a Google Earth Engine account, you can skip this step.
Visit https://signup.earthengine.google.com/ and sign-up with your Google account. You can use your existing gmail account to sign-up. It usually takes 1-2 days for approval. Hence do this step as soon as possible.
Tips to ensure smooth sign-up process:
- Use Google Chrome browser.
- When signing up for Earth Engine, please log out of all Google accounts and ensure you are logged into only 1 account which you want associated with Earth Engine.
- Access to Google Earth Engine is granted via Google Groups. The default settings should be fine, but verify you have the correct setting enabled.
- Visit groups.google.com
- Click on Settings (gear icon) and select Global Settings.
- Make sure the option Allow group managers to add me to their groups is checked.
Complete the Class Pre-Work
This class needs about 2-hours of pre-work. Please watch the following videos to get a good understanding of remote sensing and how Earth Engine works. Videos are available online and can be streamed using the video links below.
Introduction to Remote Sensing
This video introduces the remote sensing concepts, terminology and techniques.
Introduction to Google Earth Engine
This video gives a broad overview of Google Earth Engine with selected case studies and application. The video also covers the Earth Engine architecture and how it is different than traditional remote sensing software.
Take the Quizes
After you watch the videos, please complete the following 2 Quizzes
- Quiz-1 Remote Sensing Fundamentals.
- Quiz-2 Google Earth Engine Fundamentals.
Get the Course Materials
The course material and exercises are in the form of Earth Engine scripts shared via a code repository.
- Click this link to open Google Earth Engine code editor and add the repository to your account.
- If successful, you will have a new repository named
users/ujavalgandhi/End-to-End-GEEin the Scripts tab in the Reader section. - Verify that your code editor looks like below
If you do not see the repository in the Reader section, refresh your browser tab and it will show up.
Code Editor After Adding the Class Repository
Module 1: Earth Engine Basics
Module 1 is designed to give you basic skills to be able to find datasets you need for your project, filter them to your region of interest, apply basic processing and export the results. Mastering this will allow you to start using Earth Engine for your project quickly and save a lot of time pre-processing the data.
01. Hello World
This script introduces the basic Javascript syntax and the video covers the programming concepts you need to learn when using Earth Engine. To learn more, visit Introduction to JavaScript for Earth Engine section of the Earth Engine User Guide.
The Code Editor is an Integrated Development Environment (IDE) for Earth Engine Javascript API.. It offers an easy way to type, debug, run and manage code. Type the code below and click Run to execute it and see the output in the Console tab.
Tip: You can use the keyboard shortcut Ctrl+Enter to run the code in the Code Editor
Hello World
print('Hello World');
// Variables
var city = 'Bengaluru';
var country = 'India';
print(city, country);
var population = 8400000;
print(population);
// List
var majorCities = ['Mumbai', 'Delhi', 'Chennai', 'Kolkata'];
print(majorCities);
// Dictionary
var cityData = {
'city': city,
'population': 8400000,
'elevation': 930
};
print(cityData);
// Function
var greet = function(name) {
return 'Hello ' + name;
};
print(greet('World'));
// This is a commentExercise
// These are the 5 largest cities in the world:
// Tokyo, Delhi, Shanghai, Mexico City, Sao Paulo
// Create a list named 'largeCities'
// The list should have names of all the above cities
// Print the list Saving Your Work
When you modify any script for the course repository, you may want to save a copy for yourself. If you try to click the Save button, you will get an error message like below
This is because the shared class repository is a Read-only repository. You can click Yes to save a copy in your repository. If this is the first time you are using Earth Engine, you will be prompted to choose the name of your home folder. Choose the name carefully, as it cannot be changed once created.
02. Working with Image Collections
Most datasets in Earth Engine come as a ImageCollection. An ImageCollection is a dataset that consists of images takes at different time and locations - usually from the same satellite or data provider. You can load a collection by searching the Earth Engine Data Catalog for the ImageCollection ID. Search for the Sentinel-2 Level 1C dataset and you will find its id COPERNICUS/S2_SR. Visit the Sentinel-2, Level 1C page and see Explore in Earth Engine section to find the code snippet to load and visualize the collection. This snippet is a great starting point for your work with this dataset. Click the Copy Code Sample button and paste the code into the code editor. Click Run and you will see the image tiles load in the map.
In the code snippet, You will see a function Map.setCenter() which sets the viewport to a specific location and zoom level. The function takes the X coordinate (longitude), Y coordinate (latitude) and Zoom Level parameters. Replace the X and Y coordinates with the coordinates of your city and click Run to see the images of your city.
Exercise
// Find the 'Sentinel-2 Level-1C' dataset page
// https://developers.google.com/earth-engine/datasets
// Copy/page the code snippet
// Change the code to display images for your home city03. Filtering Image Collections
The collection contains all imagery ever collected by the sensor. The entire collections are not very useful. Most applications require a subset of the images. We use filters to select the appropriate images. There are many types of filter functions, look at ee.Filter... module to see all available filters. Select a filter and then run the filter() function with the filter parameters.
We will learn about 3 main types of filtering techniques
- Filter by metadata: You can apply a filter on the image metadata using filters such as
ee.Filter.eq(),ee.Filter.lt()etc. You can filter by PATH/ROW values, Orbit number, Cloud cover etc. - Filter by date: You can select images in a particular date range using filters such as
ee.Filter.date(). - Filter by location: You can select the subset of images with a bounding box, location or geometry using the
ee.Filter.bounds(). You can also use the drawing tools to draw a geometry for filtering.
After applying the filters, you can use the size() function to check how many images match the filters.
var geometry = ee.Geometry.Point([77.60412933051538, 12.952912912328241])
Map.centerObject(geometry, 10)
var s2 = ee.ImageCollection("COPERNICUS/S2");
// Filter by metadata
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
// Filter by date
var filtered = s2.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
// Filter by location
var filtered = s2.filter(ee.Filter.bounds(geometry))
// Let's apply all the 3 filters together on the collection
// First apply metadata fileter
var filtered1 = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
// Apply date filter on the results
var filtered2 = filtered1.filter(
ee.Filter.date('2019-01-01', '2020-01-01'))
// Lastly apply the location filter
var filtered3 = filtered2.filter(ee.Filter.bounds(geometry))
// Instead of applying filters one after the other, we can 'chain' them
// Use the . notation to apply all the filters together
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
print(filtered.size())Exercise
// Add one more filter in the script below to select images
// from only one of the satellites - Sentinel-2A - from the
// Sentinel-2 constellation
// Hint1: Use the 'SPACECRAFT_NAME' property
// Hint2: Use the ee.Filter.eq() filter
var geometry = ee.Geometry.Point([77.60412933051538, 12.952912912328241])
Map.centerObject(geometry, 10)
var s2 = ee.ImageCollection("COPERNICUS/S2");
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
// Remove this comment and add a filter here
print(filtered.size())04. Creating Mosaics and Composites from ImageCollections
The default order of the collection is by date. So when you display the collection, it implicitly creates a mosaic with the latest pixels on top. You can call .mosaic() on a ImageCollection to create a mosaic image from the pixels at the top.
We can also create a composite image by applying selection criteria to each pixel from all pixels in the stack. Here we use the median() function to create a composite where each pixel value is the median of all pixels from the stack.
Tip: If you need to create a mosaic where the images are in a specific order, you can use the
.sort()function to sort your collection by a property first.
Mosaic vs. Composite
var geometry = ee.Geometry.Point([77.60412933051538, 12.952912912328241])
var s2 = ee.ImageCollection("COPERNICUS/S2");
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var mosaic = filtered.mosaic()
var medianComposite = filtered.median();
Map.addLayer(filtered, rgbVis, 'Filtered Collection');
Map.addLayer(mosaic, rgbVis, 'Mosaic');
Map.addLayer(medianComposite, rgbVis, 'Median Composite')Exercise
// Create a median composite for the year 2020 and load it to the map
// Compare both the composites to see the changes in your city05. Working with Feature Collections
Feature Collections are similar to Image Collections - but they contain Features, not images. They are equivalent to Vector Layers in a GIS. We can load, filter and display Feature Collections using similar techniques that we have learned so far.
Search for GAUL Second Level Administrative Boundaries and load the collection. This is a global collection that contains all Admin2 boundaries. We can apply a filter using the ADM1_NAME property to get all Admin2 boundaries (i.e. Districts) from a state.
var admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");
var karnataka = admin2.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var visParams = {'color': 'red'}
Map.addLayer(karnataka, visParams, 'Karnataka Districts')Exercise
// Apply a filter to select only the 'Bangalore Urban' district
// Display only the selected district
// Hint: The district names are in ADM2_NAME property
var admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");
var karnataka = admin2.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var visParams = {'color': 'red'}06. Importing Data
You can import vector or raster data into Earth Engine. We will now import a shapefile of Urban Areas for Natural Earth. Unzip the ne_10m_urban_areas.zip into a folder on your computer. In the Code Editor, go to Assets → New → Table Upload → Shape Files. Select the .shp, .shx, .dbf and .prj files. Enter ne_10m_urban_areas as the Asset Name and click Upload. Once the upload and ingest finishes, you will have a new asset in the Assets tab. The shapefile is imported as a Feature Collection in Earth Engine. Select the ne_10m_urban_areas asset and click Import. You can then visualize the imported data.
Importing a Shapefile
// Let's import some data to Earth Engine
// Upload the ne_10m_urban_areas.shp shapefile
// Import the collection
var urban = ee.FeatureCollection("users/ujavalgandhi/e2e/ne_10m_urban_areas");
// Visualize the collection
Map.addLayer(urban, {color: 'blue'}, 'Urban Areas')Exercise
// Apply a filter to select only large urban areas
// Use the property 'area_sqkm' and select all features
// with area > 400 sq.km
var urban = ee.FeatureCollection("users/ujavalgandhi/e2e/ne_10m_urban_areas");07. Clipping Images
It is often desirable to clip the images to your area of interest. You can use the clip() function to mask out an image using a geometry.
While in a Desktop software, clipping is desirable to remove unnecessary portion of a large image and save computation time, in Earth Engine clipping can actually increase the computation time. As described in the Earth Engine Coding Best Practices guide, avoid clipping the images or do it at the end of your script.
Original vs. Clipped Image
var s2 = ee.ImageCollection("COPERNICUS/S2")
var urban = ee.FeatureCollection("users/ujavalgandhi/e2e/ne_10m_urban_areas")
// Find the feature id by adding the layer to the map and using Inspector.
var filtered = urban.filter(ee.Filter.eq('system:index', '00000000000000002bf8'))
var geometry = filtered.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var image = filtered.median();
var clipped = image.clip(geometry)
Map.addLayer(clipped, rgbVis, 'Clipped')Exercise
// Add the imported table to the Map
// Use the Inspector to find the id of your home city or any urban area of your choice
// Change the filter to use the id of the selected feature08. Exporting Data
Earth Engine allows for exporting both vector and raster data to be used in an external program. Vector data can be exported as a CSV or a Shapefile, while Rasters can be exported as GeoTIFF files. We will now export the Sentinel-2 Composite as a GeoTIFF file.
Tip: Code Editor supports autocompletion of API functions using the combination Ctrl+Space. Type a few characters of a function and press Ctrl+Space to see autocomplete suggestions. You can also use the same key combination to fill all parameters of the function automatically.
Once you run this script, the Tasks tab will be highlighted. Switch to the tab and you will see the tasks waiting. Click Run next to each task to start the process.
On clicking the Run button, you will be prompted for a confirmation dialog. Verify the settings and click Run to start the export.
Once the Export finishes, a GeoTiff file for each export task will be added to your Google Drive in the specified folder. You can download them and use it in a GIS software.
var s2 = ee.ImageCollection("COPERNICUS/S2")
var urban = ee.FeatureCollection("users/ujavalgandhi/e2e/ne_10m_urban_areas")
var filtered = urban.filter(ee.Filter.eq('system:index', '00000000000000002bf8'))
var geometry = filtered.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var image = filtered.median();
var clipped = image.clip(geometry)
Map.addLayer(clipped, rgbVis, 'Clipped')
var exportImage = clipped.select('B.*')
Export.image.toDrive({
image: exportImage,
description: 'Bangalore_Composite_Raw',
folder: 'earthengine',
fileNamePrefix: 'bangalore_composite_raw',
region: geometry,
scale: 20,
maxPixels: 1e9
})
// Rather than exporting raw bands, we can apply a rendered image
// visualize() function allows you to apply the same parameters
// that are used in earth engine which exports a 3-band RGB image
print(clipped)
var visualized = clipped.visualize(rgbVis)
print(visualized)
// Now the 'visualized' image is RGB image, no need to give visParams
Map.addLayer(visualized, {}, 'Visualized Image')
Export.image.toDrive({
image: visualized,
description: 'Bangalore_Composite_Visualized',
folder: 'earthengine',
fileNamePrefix: 'bangalore_composite_visualized',
region: geometry,
scale: 20,
maxPixels: 1e9
})Exercise
// Write the export function to export the results for your chosen urban areaAssignment 1
Load the Night Lights Data for 2019 and 2020. Compare the imagery for your region and find the changes in the city.
Assignment1 Expected Output
// Assignment
// Export the Night Lights images for May,2015 and May,2020
// Workflow:
// Load the VIIRS Nighttime Day/Night Band Composites collection
// Filter the collection to the date range
// Extract the 'avg_rad' band which represents the nighttime lights
// Clip the image to the geometry of your city
// Export the resulting image as a GeoTIFF file.
// Hint1:
// There are 2 VIIRS Nighttime Day/Night collections
// Use the one that corrects for stray light
// Hint2:
// The collection contains 1 global image per month
// After filtering for the month, there will be only 1 image in the collection
// You can use the following technique to extract that image
// var image = ee.Image(filtered.first())Module 2: Earth Engine Intermediate
Module 2 builds on the basic Earth Engine skills you have gained. This model introduces the parallel programming concepts using Map/Reduce - which is key in effectively using Earth Engine for analyzing large volumes of data. You will learn how to use the Earth Engine API for calculating various spectral indices, do cloud masking and then use map/reduce to do apply these computations to collections of imagery. You will also learn how to take long time-series of data and create charts.
01. Earth Engine Objects
This script introduces the basics of the Earth Engine API. When programming in Earth Engine, you must use the Earth Engine API so that your computations can use the Google Earth Engine servers. To learn more, visit Earth Engine Objects and Methods section of the Earth Engine User Guide.
// Let's see how to take a list of numbers and add 1 to each element
var myList = ee.List.sequence(1, 10);
// Define a function that takes a number and adds 1 to it
var myFunction = function(number) {
return number + 1;
}
print(myFunction(1));
//Re-Define a function using Earth Engine API
var myFunction = function(number) {
return ee.Number(number).add(1);
}
// Map the function of the list
var newList = myList.map(myFunction);
print(newList);
// Extracting value from a list
var value = newList.get(0)
print(value)
// Casting
// Let's try to do some computation on the extracted value
//var newValue = value.add(1)
//print(newValue)
// You get an error because Earth Engine doesn't know what is the type of 'value'
// We need to cast it to appropriate type first
var value = ee.Number(value)
var newValue = value.add(1)
print(newValue)
// Dictionary
// Convert javascript objects to EE Objects
var data = {'city': 'Bengaluru', 'population': 8400000, 'elevation': 930};
var eeData = ee.Dictionary(data);
// Once converted, you can use the methods from the
// ee.Dictionary module
print(eeData.get('city'))
// Dates
// For any date computation, you should use ee.Date module
var date = ee.Date('2019-01-01')
var futureDate = date.advance(1, 'year')
print(futureDate)As a general rule, you should always use Earth Engine API methods in your code, there is one exception where you will need to use client-side Javascript method. If you want to get the current time, the server doesn’t know your time. You need to use javascript method and cast it to an Earth Engine object.
var now = Date.now() print(now) var now = ee.Date(now) print(now)
Exercise
var s2 = ee.ImageCollection("COPERNICUS/S2");
var geometry = ee.Geometry.Point([77.60412933051538, 12.952912912328241]);
var now = Date.now()
var now = ee.Date(now)
// Apply another filter to the collection below to filter images
// collected in the last 1-month
// Do not hard-code the dates, it should always show images
// from the past 1-month whenever you run the script
// Hint: Use ee.Date.advance() function
// to compute the date 1 month before now
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))02. Calculating Indices
Spectral Indices are central to many aspects of remote sensing. Whether you are studying vegetation or tracking fires - you will need to compute a pixel-wise ratio of 2 or more bands. The most commonly used formula for calculating an index is the Normalized Difference between 2 bands. Earth Engine provides a helper function normalizedDifference() to help calculate normalized indices, such as Normalized Difference Vegetation Index (NDVI). For more complex formulae, you can also use the expression() function to describe the calculation.
MNDWI, SAVI and NDVI images
var s2 = ee.ImageCollection("COPERNICUS/S2");
var admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");
var bangalore = admin2.filter(ee.Filter.eq('ADM2_NAME', 'Bangalore Urban'))
var geometry = bangalore.geometry()
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var image = filtered.median();
// Calculate Normalized Difference Vegetation Index (NDVI)
// 'NIR' (B8) and 'RED' (B4)
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
// Calculate Modified Normalized Difference Water Index (MNDWI)
// 'GREEN' (B3) and 'SWIR1' (B11)
var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']);
// Calculate Soil-adjusted Vegetation Index (SAVI)
// 1.5 * ((NIR - RED) / (NIR + RED + 0.5))
// For more complex indices, you can use the expression() function
// Note:
// For the SAVI formula, the pixel values need to converted to reflectances
// Multiplyng the pixel values by 'scale' gives us the reflectance value
// The scale value is 0.0001 for Sentinel-2 dataset
var savi = image.expression(
'1.5 * ((NIR - RED) / (NIR + RED + 0.5))', {
'NIR': image.select('B8').multiply(0.0001),
'RED': image.select('B4').multiply(0.0001),
}).rename('savi');
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
var ndviVis = {min:0, max:1, palette: ['white', 'green']}
var ndwiVis = {min:0, max:0.5, palette: ['white', 'blue']}
Map.addLayer(image.clip(geometry), rgbVis, 'Image');
Map.addLayer(mndwi.clip(geometry), ndwiVis, 'mndwi')
Map.addLayer(savi.clip(geometry), ndviVis, 'savi')
Map.addLayer(ndvi.clip(geometry), ndviVis, 'ndvi')Exercise
var s2 = ee.ImageCollection("COPERNICUS/S2");
var admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");
var bangalore = admin2.filter(ee.Filter.eq('ADM2_NAME', 'Bangalore Urban'))
var geometry = bangalore.geometry()
Map.centerObject(geometry)
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var image = filtered.median();
// Exercise
// Calculate the Normalized Difference Built-Up Index (NDBI) for the image
// Hint: NDBI = (SWIR1 – NIR) / (SWIR1 + NIR)
// Visualize the built-up area using a 'red' palette03. Computation on ImageCollections
So far we have learnt how to run computation on single images. If you want to apply some computation - such as calculating an index - to many images, you need to use map(). You first define a function that takes 1 image and returns the result of the computation on that image. Then you can map() that function over the ImageCollection which results in a new ImageCollection with the results of the computation. This is similar to a for-loop that you maybe familiar with - but using map() allows the computation to run in parallel. Learn more at Mapping over an ImageCollection
NDVI computation on an ImageCollection
var s2 = ee.ImageCollection("COPERNICUS/S2");
var admin1 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level1");
var karnataka = admin1.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var geometry = karnataka.geometry()
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var composite = filtered.median().clip(geometry)
Map.addLayer(composite, rgbVis, 'Karnataka Composite')
// Write a function that computes NDVI for an image and adds it as a band
function addNDVI(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
return image.addBands(ndvi);
}
// Map the function over the collection
var withNdvi = filtered.map(addNDVI);
var composite = withNdvi.median()
var ndviComposite = composite.select('ndvi').clip(karnataka)
var palette = [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
'74A901', '66A000', '529400', '3E8601', '207401', '056201',
'004C00', '023B01', '012E01', '011D01', '011301'];
var ndviVis = {min:0, max:0.5, palette: palette }
Map.addLayer(ndviComposite, ndviVis, 'ndvi')Exercise
var s2 = ee.ImageCollection("COPERNICUS/S2");
var admin1 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level1");
var karnataka = admin1.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var geometry = karnataka.geometry()
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
var composite = filtered.median().clip(geometry)
Map.addLayer(composite, rgbVis, 'Karnataka Composite')
// This function calculates both NDVI an d NDWI indices
// and returns an image with 2 new bands added to the original image.
function addIndices(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
var ndwi = image.normalizedDifference(['B3', 'B8']).rename('ndwi');
return image.addBands(ndvi).addBands(ndwi);
}
// Map the function over the collection
var withIndices = filtered.map(addIndices);
// Composite
var composite = withIndices.median()
print(composite)
// Extract the 'ndwi' band and display a NDWI map
// use the palette ['white', 'blue']
// Hint: Use .select() function to select a band04. Cloud Masking
Masking pixels in an image makes those pixels transparent and excludes them from analysis and visualization. To mask an image, we can use the updateMask() function and pass it an image with 0 and 1 values. All pixels where the mask image is 0 will be masked.
Most remote sensing datasets come with a QA or Cloud Mask band that contains the information on whether pixels is cloudy or not. Your Code Editor contains pre-defined functions for masking clouds for popular datasets under Scripts Tab → Examples → Cloud Masking.
The script below takes the Sentinel-2 masking function and shows how to apply it on an image.
Applying pixel-wise QA bitmask
var image = ee.Image('COPERNICUS/S2/20190703T050701_20190703T052312_T43PGP')
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
Map.centerObject(image)
Map.addLayer(image, rgbVis, 'Full Image', false)
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var maskedImage = ee.Image(maskS2clouds(image))
Map.addLayer(maskedImage, rgbVis, 'Masked Image')Exercise
// Get the Level-2A Surface Reflectance image
var imageSR = ee.Image('COPERNICUS/S2_SR/20190703T050701_20190703T052312_T43PGP')
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
Map.centerObject(imageSR)
Map.addLayer(imageSR, rgbVis, 'SR Image')
// Function to remove cloud and snow pixels from Sentinel-2 SR image
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(5);
var snow = snowProb.lt(5);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 5% or cloud shadow classification
var mask = (cloud.and(snow)).and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask);
}
// Exercise
// Apply the above cloud masking function to SR image
// Add the masked image to the map
// Hint: After adding the masked image to the map, turn-off
// the original image layer to see the result of the masking functionIf you are using Sentinel-2 data, do check out the an alternative cloud masking techninque using the S2 Cloudless dataset. Learn more
var imageSR = ee.Image('COPERNICUS/S2_SR/20190703T050701_20190703T052312_T43PGP') var s2Cloudless = ee.Image('COPERNICUS/S2_CLOUD_PROBABILITY/20190703T050701_20190703T052312_T43PGP') var clouds = s2Cloudless.lt(50) var cloudlessMasked = imageSR.mask(clouds) var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']}; Map.addLayer(cloudlessMasked, rgbVis, 'S2 Cloudless Masked Image')
05. Reducers
When writing parallel computing code, a Reduce operation allows you to compute statistics on a large amount of inputs. In Earth Engine, you need to run reduction operation when creating composites, calculating statistics, doing regression analysis etc. The Earth Engine API comes with a large number of built-in reducer functions (such as ee.Reducer.sum(), ee.Reducer.histogram(), ee.Reducer.linearFit() etc.) that can perform a variety of statistical operations on input data. You can run reducers using the reduce() function. Earth Engine supports running reducers on all data structures that can hold multiple values, such as Images (reducers run on different bands), ImageCollection, FeatureCollection, List, Dictionary etc. The script below introduces basic concepts related to reducers.
// Computing stats on a list
var myList = ee.List.sequence(1, 10);
print(myList)
// Use a reducer to compute min and max in the list
var mean = myList.reduce(ee.Reducer.mean());
print(mean);
var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
var s2 = ee.ImageCollection("COPERNICUS/S2");
Map.centerObject(geometry)
// Apply a reducer on a image collection
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry))
.select('B.*')
print(filtered.size())
var collMean = filtered.reduce(ee.Reducer.mean());
print('Reducer on Collection', collMean);
var image = ee.Image('COPERNICUS/S2/20190223T050811_20190223T051829_T44RPR')
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
Map.addLayer(image, rgbVis, 'Image')
Map.addLayer(geometry, {color: 'red'}, 'Farm')
// If we want to compute the average value in each band,
// we can use reduceRegion instead
var stats = image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: geometry,
scale: 100,
maxPixels: 1e10
})
print(stats);
// Result of reduceRegion is a dictionary.
// We can extract the values using .get() function
print('Average value in B4', stats.get('B4'))Exercise
var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
var image = ee.Image('COPERNICUS/S2/20190223T050811_20190223T051829_T44RPR')
Map.addLayer(image, rgbVis, 'Image')
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
// Exercise
// Compute the average NDVI for the farm from the given image
// Hint: Use the reduceRegion() function06. Time-Series Charts
Now we can put together all the skills we have learnt so far - filter, map, reduce, and cloud-masking to create a chart of average NDVI values for a given farm over 1 year. Earth Engine API comes with support for charting functions based on the Google Chart API. Here we use the ui.Chart.image.series() function to create a time-series chart.
Computing NDVI Time-series for a Farm
NDVI Time-series showing Dual-Cropping Cycle
var s2 = ee.ImageCollection("COPERNICUS/S2");
var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
var filtered = s2
.filter(ee.Filter.date('2017-01-01', '2018-01-01'))
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)//.divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = filtered.map(maskS2clouds)
// Write a function that computes NDVI for an image and adds it as a band
function addNDVI(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
return image.addBands(ndvi);
}
// Map the function over the collection
var withNdvi = filtered.map(addNDVI);
// Display a time-series chart
var chart = ui.Chart.image.series({
imageCollection: withNdvi.select('ndvi'),
region: geometry,
reducer: ee.Reducer.mean(),
scale: 20
}).setOptions({
lineWidth: 1,
title: 'NDVI Time Series',
interpolateNulls: true,
vAxis: {title: 'NDVI'},
hAxis: {title: '', format: 'YYYY-MMM'}
})
print(chart);Exercise
// Delete the farm boundary from the previous script
// and add another farm at a location of your choice
// Print the chart.Assignment 2
Assignment2 Expected Output
var terraclimate = ee.ImageCollection("IDAHO_EPSCOR/TERRACLIMATE");
var geometry = ee.Geometry.Point([77.54849920033682, 12.91215102400037]);
// Assignment
// Use Gridded Climate dataset to chart a 40+ year time series
// if temparature at any location
// Workflow
// Load the TerraClimate collection
// Select the 'tmmx' band
// Scale the band values
// Filter the scaled collection
// Use ui.Chart.image.series() function to create the chart
// Hint1
// Data needed to be scaled by 0.1
// map() a function and multiply each image
// Multiplying creates a new image that doesn't have the same properties
// Use copyProperties() function to copy timestamp to new image
var tmax = terraclimate.select('tmmx')
var tmaxScaled = tmax.map(function(image) {
return image.multiply(0.1)
.copyProperties(image,['system:time_start']);
})
// Hint2
// You will need to specify a scale in meters for charting
// Use projection().nominalScale() to find the
// image resolution in meters
var image = ee.Image(terraclimate.first())
print(image.projection().nominalScale())Module 3: Supervised Classification
Introduction to Machine Learning and Supervised Classification
Supervised classification is arguably the most important classical machine learning techniques in remote sensing. Applications range from generating Land Use/Land Cover maps to change detection. Google Earth Engine is unique suited to do supervised classification at scale. The interactive nature of Earth Engine development allows for iterative development of supervised classification workflows by combining many different datasets into the model. This module covers basic supervised classification workflow, accuracy assessment, hyperparameter tuning and change detection.
01. Basic Supervised Classification
We will learn how to do a basic land cover classification using training samples collected from the Code Editor using the High Resolution basemap imagery provided by Google Maps. This method requires no prior training data and is quite effective to generate high quality classification samples anywhere in the world. The goal is to classify each source pixel into one of the following classes - urban, bare, water or vegetation. Using the drawing tools in the code editor, you create 4 new feature collection with points representing pixels of that class. Each feature collection has a property called landcover with values of 0, 1, 2 or 3 indicating whether the feature collection represents urban, bare, water or vegetation respectively. We then train a Random Forest classifier using these training set to build a model and apply it to all the pixels of the image to create a 4 class image.
Fun fact: The classifiers in Earth Engine API have names starting with smile - such as
ee.Classifier.smileRandomForest(). The smile part refers to the Statistical Machine Intelligence and Learning Engine (SMILE) JAVA library which is used by Google Earth Engine to implement these algorithms.
Supervised Classification Output
var bangalore = ee.FeatureCollection("users/ujavalgandhi/public/bangalore_boundary")
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
// The following collections were created using the
// Drawing Tools in the code editor
var urban = ee.FeatureCollection("users/ujavalgandhi/e2e/urban_gcps")
var bare = ee.FeatureCollection("users/ujavalgandhi/e2e/bare_gcps")
var water = ee.FeatureCollection("users/ujavalgandhi/e2e/water_gcps")
var vegetation = ee.FeatureCollection("users/ujavalgandhi/e2e/vegetation_gcps")
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(bangalore))
.select('B.*')
var composite = filtered.median().clip(bangalore)
// Display the input composite.
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
Map.addLayer(composite, rgbVis, 'image');
var gcps = urban.merge(bare).merge(water).merge(vegetation)
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: gcps,
properties: ['landcover'],
scale: 10
});
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50).train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// // Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
// Display the GCPs
// We use the style() function to style the GCPs
var palette = ee.List(['gray','brown','blue','green'])
var landcover = ee.List([0, 1, 2, 3])
var gcpsStyled = ee.FeatureCollection(
landcover.map(function(lc){
var color = palette.get(landcover.indexOf(lc));
var markerStyle = { color: 'white', pointShape: 'diamond',
pointSize: 4, width: 1, fillColor: color}
return gcps.filter(ee.Filter.eq('landcover', lc))
.map(function(point){
return point.set('style', markerStyle)
})
})).flatten();
Map.addLayer(gcpsStyled.style({styleProperty:"style"}), {}, 'GCPs')
Map.centerObject(gcpsStyled)Exercise
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
var urbanAreas = ee.FeatureCollection("users/ujavalgandhi/e2e/ne_10m_urban_areas")
// Perform supervised classification for your city
// Find the feature id by adding the layer to the map and using Inspector.
var city = urbanAreas.filter(ee.Filter.eq('system:index', '00000000000000002bf8'))
var geometry = city.geometry()
Map.centerObject(geometry)
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(geometry))
.select('B.*')
var composite = filtered.median().clip(geometry)
// Display the input composite.
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
Map.addLayer(composite, rgbVis, 'image');
// Exercise
// Add training points for 4 classes
// Assign the 'landcover' property as follows
// urban: 0
// bare: 1
// water: 2
// vegetation: 3
// After adding points, uncomments lines below
// var gcps = urban.merge(bare).merge(water).merge(vegetation)
// // Overlay the point on the image to get training data.
// var training = composite.sampleRegions({
// collection: gcps,
// properties: ['landcover'],
// scale: 10,
// tileScale: 16
// });
// print(training)
// // Train a classifier.
// var classifier = ee.Classifier.smileRandomForest(50).train({
// features: training,
// classProperty: 'landcover',
// inputProperties: composite.bandNames()
// });
// // // Classify the image.
// var classified = composite.classify(classifier);
// Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019'); 02. Accuracy Assessment
It is important to get a quantitative estimate of the accuracy of the classification. To do this, a common strategy is to divide your training samples into 2 random fractions - one used for training the model and the other for validation of the predictions. Once a classifier is trained, it can be used to classify the entire image. We can then compare the classified values with the ones in the validation fraction. We can use the ee.Classifier.confusionMatrix() method to calculate a Confusion Matrix representing expected accuracy.
Classification results are evaluated based on the following metrics
- Overall Accuracy: How many samples were classified correctly.
- Producer’s Accuracy: How well did the classification predict each class.
- Consumer’s Accuracy (Reliability): How reliable is the prediction in each class.
- Kappa Coefficient: How well the classification performed as compared to random assignment.
Accuracy Assessment
Don’t get carried away tweaking your model to give you the highest validation accuracy. You must use both qualitative measures (such as visual inspection of results) along with quantitative measures to assess the results.
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7");
var gcp = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_gcps");
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.select('B.*')
var composite = filtered.median().clip(boundary)
// Display the input composite.
Map.addLayer(composite, rgbVis, 'image');
// Add a random column and split the GCPs into training and validation set
var gcp = gcp.randomColumn()
// This being a simpler classification, we take 60% points
// for validation. Normal recommended ratio is
// 70% training, 30% validation
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: trainingGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
//**************************************************************************
// Accuracy Assessment
//**************************************************************************
// Use classification map to assess accuracy using the validation fraction
// of the overall training set created above.
var test = classified.sampleRegions({
collection: validationGcp,
properties: ['landcover'],
tileScale: 16,
scale: 10,
});
var testConfusionMatrix = test.errorMatrix('landcover', 'classification')
// Printing of confusion matrix may time out. Alternatively, you can export it as CSV
print('Confusion Matrix', testConfusionMatrix);
print('Test Accuracy', testConfusionMatrix.accuracy());Exercise
var classified = ee.Image('users/ujavalgandhi/e2e/arkavathy_base_classification');
var gcp = ee.FeatureCollection('users/ujavalgandhi/e2e/arkavathy_gcps');
var gcp = gcp.randomColumn()
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
//**************************************************************************
// Accuracy Assessment
//**************************************************************************
// Use classification map to assess accuracy using the validation fraction
// of the overall training set created above.
var test = classified.sampleRegions({
collection: validationGcp,
properties: ['landcover'],
tileScale: 16,
scale: 10,
});
var testConfusionMatrix = test.errorMatrix('landcover', 'classification')
print('Confusion Matrix', testConfusionMatrix);
print('Test Accuracy', testConfusionMatrix.accuracy());
// Exercise
// Calculate and print the following assessment metrics
// 1. Producer's accuracy
// 2. Consumer's accuracy
// 3. Kappa coefficient
// Hint: Look at the ee.ConfusionMatrix module for appropriate methods03. Improving the Classification
The Earth Engine data-model is especially well suited for machine learning tasks because of its ability to easily incorporate data sources of different spatial resolutions, projections and data types together By giving additional information to the classifier, it is able to separate different classes easily. Here we take the same example and augment it with the following techniques
- Apply Cloud Masking
- Add Spectral Indices: We add bands for different spectral indices such as - NDVI, NDBI, MNDWI and BSI.
- Add Elevation and Slope: We also add slope and elevation bands from the ALOS DEM.
- Normalize the Inputs: Machine learning models work best when all the inputs have the same scale. We will divide each band with the maximum value. This method ensures that all input values are between 0-1. A more complete and robust technique for image normalization is provided in the course Supplement.
Our training features have more parameters and contain values of the same scale. The result is a much improved classification.
Improved Classification Accuracy with use of Spectral Indices and Elevation Data
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7");
var gcp = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_gcps");
var alos = ee.Image("JAXA/ALOS/AW3D30/V2_2");
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Function to remove cloud and snow pixels from Sentinel-2 SR image
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(10);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 10% or cloud shadow classification
var mask = cloud.and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask);
}
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.map(maskCloudAndShadowsSR)
.select('B.*')
var composite = filtered.median().clip(boundary)
var addIndices = function(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']);
var bsi = image.expression(
'(( X + Y ) - (A + B)) /(( X + Y ) + (A + B)) ', {
'X': image.select('B11'), //swir1
'Y': image.select('B4'), //red
'A': image.select('B8'), // nir
'B': image.select('B2'), // blue
}).rename('bsi');
return image.addBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi)
}
var composite = addIndices(composite);
// Calculate Slope and Elevation
var elev = alos.select('AVE_DSM').rename('elev');
var slope = ee.Terrain.slope(alos.select('AVE_DSM')).rename('slope');
var composite = composite.addBands(elev).addBands(slope);
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.2};
Map.addLayer(composite, visParams, 'RGB');
// Normalize the image
// Machine learning algorithms work best on images when all features have
// the same range
// Function to Normalize Image
// Pixel Values should be between 0 and 1
// Formula is (x - xmin) / (xmax - xmin)
//**************************************************************************
function normalize(image){
var bandNames = image.bandNames();
// Compute min and max of the image
var minDict = image.reduceRegion({
reducer: ee.Reducer.min(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var maxDict = image.reduceRegion({
reducer: ee.Reducer.max(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var mins = ee.Image.constant(minDict.values(bandNames));
var maxs = ee.Image.constant(maxDict.values(bandNames));
var normalized = image.subtract(mins).divide(maxs.subtract(mins))
return normalized
}
var composite = normalize(composite);
// Add a random column and split the GCPs into training and validation set
var gcp = gcp.randomColumn()
// This being a simpler classification, we take 60% points
// for validation. Normal recommended ratio is
// 70% training, 30% validation
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: trainingGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
print(training)
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
//**************************************************************************
// Accuracy Assessment
//**************************************************************************
// Use classification map to assess accuracy using the validation fraction
// of the overall training set created above.
var test = classified.sampleRegions({
collection: validationGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
var testConfusionMatrix = test.errorMatrix('landcover', 'classification')
// Printing of confusion matrix may time out. Alternatively, you can export it as CSV
print('Confusion Matrix', testConfusionMatrix);
print('Test Accuracy', testConfusionMatrix.accuracy());
Exercise
// Exercise
// Improve your classification from Exercise 01c
// Add different spectral indicies to your composite
// by using the function below
var addIndices = function(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']);
var bsi = image.expression(
'(( X + Y ) - (A + B)) /(( X + Y ) + (A + B)) ', {
'X': image.select('B11'), //swir1
'Y': image.select('B4'), //red
'A': image.select('B8'), // nir
'B': image.select('B2'), // blue
}).rename('bsi');
return image.addBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi)
}
04. Exporting Classification Results
When working with complex classifiers over large regions, you may get a User memory limit exceeded or Computation timed out error in the Code Editor. The reason for this is that there is a fixed time limit and smaller memory allocated for code that is run with the On-Demand Computation mode. For larger computations, you can use the Batch mode with the Export functions. Exports run in the background and can run longer than 5-minutes time allocated to the computation code run from the Code Editor. This allows you to process very large and complex datasets. Here’s an example showing how to export your classification results to Google Drive.
We can only export Images or FeatureCollections. What if you wanted to export a number that is the result of a long computation? A useful hack is to create a FeatureCollection with just 1 feature containing
nullgeometry and a property containing the number you want to export.
Exported Classification Outputs
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7");
var gcp = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_gcps");
var alos = ee.Image("JAXA/ALOS/AW3D30/V2_2");
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Function to remove cloud and snow pixels from Sentinel-2 SR image
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(10);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 10% or cloud shadow classification
var mask = cloud.and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask);
}
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.map(maskCloudAndShadowsSR)
.select('B.*')
var composite = filtered.median().clip(boundary)
var addIndices = function(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']);
var bsi = image.expression(
'(( X + Y ) - (A + B)) /(( X + Y ) + (A + B)) ', {
'X': image.select('B11'), //swir1
'Y': image.select('B4'), //red
'A': image.select('B8'), // nir
'B': image.select('B2'), // blue
}).rename('bsi');
return image.addBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi)
}
var composite = addIndices(composite);
// Calculate Slope and Elevation
var elev = alos.select('AVE_DSM').rename('elev');
var slope = ee.Terrain.slope(alos.select('AVE_DSM')).rename('slope');
var composite = composite.addBands(elev).addBands(slope);
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.2};
Map.addLayer(composite, visParams, 'RGB');
// Normalize the image
// Machine learning algorithms work best on images when all features have
// the same range
// Function to Normalize Image
// Pixel Values should be between 0 and 1
// Formula is (x - xmin) / (xmax - xmin)
//**************************************************************************
function normalize(image){
var bandNames = image.bandNames();
// Compute min and max of the image
var minDict = image.reduceRegion({
reducer: ee.Reducer.min(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var maxDict = image.reduceRegion({
reducer: ee.Reducer.max(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var mins = ee.Image.constant(minDict.values(bandNames));
var maxs = ee.Image.constant(maxDict.values(bandNames));
var normalized = image.subtract(mins).divide(maxs.subtract(mins))
return normalized
}
var composite = normalize(composite);
// Add a random column and split the GCPs into training and validation set
var gcp = gcp.randomColumn()
// This being a simpler classification, we take 60% points
// for validation. Normal recommended ratio is
// 70% training, 30% validation
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
Map.addLayer(validationGcp)
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: trainingGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
print(training)
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
//**************************************************************************
// Accuracy Assessment
//**************************************************************************
// Use classification map to assess accuracy using the validation fraction
// of the overall training set created above.
var test = classified.sampleRegions({
collection: validationGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
var testConfusionMatrix = test.errorMatrix('landcover', 'classification')
print('Confusion Matrix', testConfusionMatrix);
print('Test Accuracy', testConfusionMatrix.accuracy());
//**************************************************************************
// Exporting Results
//**************************************************************************
// Create a Feature with null geometry and the value we want to export.
// Use .array() to convert Confusion Matrix to an Array so it can be
// exported in a CSV file
var fc = ee.FeatureCollection([
ee.Feature(null, {
'accuracy': testConfusionMatrix.accuracy(),
'matrix': testConfusionMatrix.array()
})
])
print(fc)
Export.table.toDrive({
collection: fc,
description: 'Accuracy_Export',
folder: 'earthengine',
fileNamePrefix: 'accuracy',
fileFormat: 'CSV'
})Exercise
It is also a good idea to export the classified image as an Asset. This will allows you to import the classified image in another script without running the whole classification workflow. Use the Export.image.toAsset() function to export the classified image as an asset.
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7");
var gcp = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_gcps");
var alos = ee.Image("JAXA/ALOS/AW3D30/V2_2");
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Function to remove cloud and snow pixels from Sentinel-2 SR image
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(10);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 10% or cloud shadow classification
var mask = cloud.and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask);
}
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.map(maskCloudAndShadowsSR)
.select('B.*')
var composite = filtered.median().clip(boundary)
var addIndices = function(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']);
var bsi = image.expression(
'(( X + Y ) - (A + B)) /(( X + Y ) + (A + B)) ', {
'X': image.select('B11'), //swir1
'Y': image.select('B4'), //red
'A': image.select('B8'), // nir
'B': image.select('B2'), // blue
}).rename('bsi');
return image.addBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi)
}
var composite = addIndices(composite);
// Calculate Slope and Elevation
var elev = alos.select('AVE_DSM').rename('elev');
var slope = ee.Terrain.slope(alos.select('AVE_DSM')).rename('slope');
var composite = composite.addBands(elev).addBands(slope);
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.2};
Map.addLayer(composite, visParams, 'RGB');
// Normalize the image
// Machine learning algorithms work best on images when all features have
// the same range
// Function to Normalize Image
// Pixel Values should be between 0 and 1
// Formula is (x - xmin) / (xmax - xmin)
//**************************************************************************
function normalize(image){
var bandNames = image.bandNames();
// Compute min and max of the image
var minDict = image.reduceRegion({
reducer: ee.Reducer.min(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var maxDict = image.reduceRegion({
reducer: ee.Reducer.max(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var mins = ee.Image.constant(minDict.values(bandNames));
var maxs = ee.Image.constant(maxDict.values(bandNames));
var normalized = image.subtract(mins).divide(maxs.subtract(mins))
return normalized
}
var composite = normalize(composite);
// Add a random column and split the GCPs into training and validation set
var gcp = gcp.randomColumn()
// This being a simpler classification, we take 60% points
// for validation. Normal recommended ratio is
// 70% training, 30% validation
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
Map.addLayer(validationGcp)
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: trainingGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
// Exercise
// Use the Export.image.toAsset() function to export the
// classified image as a Earth Engine Asset.
// This will allows you to import the classified image in another script
// without running the whole classification workflow.
// Hint: For images with discrete pixel values, we must set the
// pyramidingPolicy to 'mode'.
// The pyramidingPolicy parameter should a dictionary specifying
// the policy for each band. A simpler way to specify it for all
// bands is to use {'.default': 'mode'}05. Calculating Area
Now that we have the results of our classification, we will learn how to calculate the area for pixels in each class. Calculating area for features is done using the area() function and for images using the ee.Image.pixelArea() function. The ee.Image.pixelArea() function creates an image where each pixel’s value is the area of the pixel. We multiply this pixel area image with our image and sum up the area using the reduceRegion() function.
The
ee.Image.pixelArea()function uses a custom equal-area projection for area calculation. The result is area in square meters regardless of the projection of the input image. Learn more
Calculating Green Cover from Classified Image
var classified = ee.Image("users/ujavalgandhi/e2e/bangalore_classified");
var bangalore = ee.FeatureCollection("users/ujavalgandhi/public/bangalore_boundary");
Map.addLayer(bangalore, {color: 'blue'}, 'Bangalore City')
Map.addLayer(classified,
{min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']},
'Classified Image 2019');
// Calling .geometry() on a feature collection gives the
// dissolved geometry of all features in the collection
// .area() function calculates the area in square meters
var cityArea = bangalore.geometry().area()
// We can cast the result to a ee.Number() and calculate the
// area in square kilometers
var cityAreaSqKm = ee.Number(cityArea).divide(1e6).round()
print(cityAreaSqKm)
// Area Calculation for Images
var vegetation = classified.eq(3)
// If the image contains values 0 or 1, we can calculate the
// total area using reduceRegion() function
// The result of .eq() operation is a binary image with pixels
// values of 1 where the condition matched and 0 where it didn't
Map.addLayer(vegetation, {min:0, max:1, palette: ['white', 'green']}, 'Green Cover')
// Since our image has only 0 and 1 pixel values, the vegetation
// pixels will have values equal to their area
var areaImage = vegetation.multiply(ee.Image.pixelArea())
// Now that each pixel for vegetation class in the image has the value
// equal to its area, we can sum up all the values in the region
// to get the total green cover.
var area = areaImage.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: bangalore.geometry(),
scale: 10,
maxPixels: 1e10
})
// The result of the reduceRegion() function is a dictionary with the key
// being the band name. We can extract the area number and convert it to
// square kilometers
var vegetationAreaSqKm = ee.Number(area.get('classification')).divide(1e6).round()
print(vegetationAreaSqKm)If you want to compute area covered by each class, you can use a Grouped Reducer. See the Supplement to see a code snippet.
Exercise
// Exercise
// Compute and print the percentage green cover of the cityModule 4: Change Detection
Introduction to Change Detection
Many earth observation datasets are available at regular intervals over long periods of time. This enables us to detect changes on the Earth’s surface. Change detection technique in remote sensing fall in the following categories
- Single Band Change: Measuring change in a single band image or a spectral index using a threshold
- Multi Band Change: Measuring spectral distance and spectral angle between two multiband images
- Classification of Change: One-pass classification using stacked image containing bands from before and after an event
- Post Classification Comparison: Comparing two classified images and computing class transitions
In the following sections, we will apply the supervised classification techniques for change detection using multi-temporal images.
01. Spectral Distance Change
When you want to detect changes from multi-band images, a useful technique is to compute the Spectral Distance and Spectral Angle between the two images. Pixels that exhibit a large change will have a larger distance compared to those that did not change. This technique is particularly useful when there are no suitable index to detect the change. It can be applied to detect change after natural disasters or human conflicts.
Here we use this technique to detect landslides using before/after composites. You may learn more about this technique at Craig D’Souza’s Change Detection presentation.
Spectral Distance Change Detection
var geometry = ee.Geometry.Polygon([[
[75.70357667713435, 12.49723970868507],
[75.70357667713435, 12.470171844429931],
[75.7528434923199, 12.470171844429931],
[75.7528434923199, 12.49723970868507]
]]);
Map.centerObject(geometry)
var s2 = ee.ImageCollection("COPERNICUS/S2")
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)//.divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = s2
.filter(ee.Filter.bounds(geometry))
.map(maskS2clouds)
.select('B.*')
var dateOfIncident = ee.Date('2018-12-15')
var before = filtered
.filter(ee.Filter.date(dateOfIncident.advance(-2, 'year'), dateOfIncident))
.filter(ee.Filter.calendarRange(12, 12, 'month'))
.median()
var after = filtered.filter(ee.Filter.date(
dateOfIncident, dateOfIncident.advance(1, 'month'))).median()
Map.addLayer(before, rgbVis, 'Before')
Map.addLayer(after, rgbVis, 'After')
// Calculate Distance
// Formula at https://www.varsitytutors.com/calculus_3-help/distance-between-vectors
var magnitude = function(image) {
return image.pow(2).reduce(ee.Reducer.sum()).sqrt();
};
var distance = magnitude(after.subtract(before))
// Calculate Angle
// https://byjus.com/angle-between-two-vectors-formula/
var dot = before.multiply(after).reduce(ee.Reducer.sum());
var angle = dot.divide(magnitude(after))
.divide(magnitude(before))
.acos()
Map.addLayer(distance, {min: 0, max: 1500, palette: ['white', 'red']}, 'spectral distance');
Map.addLayer(angle, {min: 0, max: 1, palette: ['white', 'purple']}, 'angle');Exercise
var geometry = ee.Geometry.Polygon([[
[75.70357667713435, 12.49723970868507],
[75.70357667713435, 12.470171844429931],
[75.7528434923199, 12.470171844429931],
[75.7528434923199, 12.49723970868507]
]]);
Map.centerObject(geometry)
var s2 = ee.ImageCollection("COPERNICUS/S2")
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)//.divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = s2
.filter(ee.Filter.bounds(geometry))
.map(maskS2clouds)
.select('B.*')
var dateOfIncident = ee.Date('2018-12-15')
var before = filtered
.filter(ee.Filter.date(dateOfIncident.advance(-2, 'year'), dateOfIncident))
.filter(ee.Filter.calendarRange(12, 12, 'month'))
.median()
var after = filtered.filter(ee.Filter.date(
dateOfIncident, dateOfIncident.advance(1, 'month'))).median()
Map.addLayer(before, rgbVis, 'Before')
Map.addLayer(after, rgbVis, 'After')
var magnitude = function(image) {
return image.pow(2).reduce(ee.Reducer.sum()).sqrt();
};
var distance = magnitude(after.subtract(before))
Map.addLayer(distance, {min: 0, max: 1500, palette: ['white', 'red']}, 'spectral distance');
// Exercise
// Inspect the distance image and find a suitable threshold
// that signifies damage after the landslides
// Apply the threshold and create a new image showing landslides
// Display the results
// Hint: Use the .gt() method to apply the threshold02. Direct Classification of Change
This technique of change detection is also known as One-pass Classification or Direct Multi-date Classification. Here we create a single stacked image containing bands from before and after images. We train a classifier with training data sampled from the stacked image and apply the classifier on the stacked image to find all change pixels.
All pixels that changed from bare ground to built-up
var bangalore = ee.FeatureCollection('users/ujavalgandhi/public/bangalore_boundary');
var change = ee.FeatureCollection('users/ujavalgandhi/e2e/bangalore_change_gcps');
var nochange = ee.FeatureCollection('users/ujavalgandhi/e2e/bangalore_nochange_gcps')
var s2 = ee.ImageCollection("COPERNICUS/S2")
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)//.divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = s2
.filter(ee.Filter.date('2019-01-01', '2019-02-01'))
.filter(ee.Filter.bounds(bangalore))
.map(maskS2clouds)
var image2019 = filtered.median().clip(bangalore)
// Display the input composite.
Map.addLayer(image2019, rgbVis, '2019');
var filtered = s2
.filter(ee.Filter.date('2020-01-01', '2020-02-01'))
.filter(ee.Filter.bounds(bangalore))
.map(maskS2clouds)
var image2020 = filtered.median().clip(bangalore)
Map.addLayer(image2020, rgbVis, '2020');
var stackedImage = image2019.addBands(image2020)
// Overlay the point on the image to get training data.
var training = stackedImage.sampleRegions({
collection: change.merge(nochange),
properties: ['class'],
scale: 10
});
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50).train({
features: training,
classProperty: 'class',
inputProperties: stackedImage.bandNames()
});
// Classify the image.
var classified = stackedImage.classify(classifier);
Map.addLayer(classified, {min: 0, max: 1, palette: ['white', 'red']}, 'change'); Exercise
// Add an NDBI band to improve the detection of changes.
var addNDBI = function(image) {
var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
return image.addBands(ndbi)
}
// use addNDBI() function to add the NDBI band to both 2019 and 2020 composite images
// Hint1: You can save the resulting image in the same variable to avoid changing
// a lot of code.
// var image = addNDBI(image)03. Post-classification Comparison
We dealing with multi-class images, a useful metric for change detection is to know how many pixels from class X changed to class Y. This can be accomplished using the ee.Reducer.frequencyHistogram() reducer as shown below.
var bangalore = ee.FeatureCollection("users/ujavalgandhi/public/bangalore_boundary")
var urban = ee.FeatureCollection("users/ujavalgandhi/e2e/urban_gcps")
var bare = ee.FeatureCollection("users/ujavalgandhi/e2e/bare_gcps")
var water = ee.FeatureCollection("users/ujavalgandhi/e2e/water_gcps")
var vegetation = ee.FeatureCollection("users/ujavalgandhi/e2e/vegetation_gcps")
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
Map.centerObject(bangalore)
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
var filtered = s2
.filter(ee.Filter.date('2019-01-01', '2019-01-31'))
.filter(ee.Filter.bounds(bangalore))
.select('B.*')
var before = filtered.median().clip(bangalore)
// Display the input composite.
Map.addLayer(before, rgbVis, 'before');
var training = urban.merge(bare).merge(water).merge(vegetation)
// Overlay the point on the image to get training data.
var training = before.sampleRegions({
collection: training,
properties: ['landcover'],
scale: 10
});
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50).train({
features: training,
classProperty: 'landcover',
inputProperties: before.bandNames()
});
// // Classify the image.
var beforeClassified = before.classify(classifier);
Map.addLayer(beforeClassified,
{min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, 'before_classified');
// 2020 Jan
var after = s2
.filter(ee.Filter.date('2020-01-01', '2020-01-31'))
.filter(ee.Filter.bounds(bangalore))
.select('B.*')
.median()
.clip(bangalore)
Map.addLayer(after, rgbVis, 'after');
// Classify the image.
var afterClassified= after.classify(classifier);
Map.addLayer(afterClassified,
{min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, 'after_classified');
// Reclassify from 0-3 to 1-4
var beforeClasses = beforeClassified.remap([0, 1, 2, 3], [1, 2, 3, 4])
var afterClasses = afterClassified.remap([0, 1, 2, 3], [1, 2, 3, 4])
// Show all changed areas
var changed = afterClasses.subtract(beforeClasses).neq(0)
Map.addLayer(changed, {min:0, max:1, palette: ['white', 'red']}, 'Change')
// We multiply the before image with 100 and add the after image
// The resulting pixel values will be unique and will represent each unique transition
// i.e. 102 is urban to bare, 103 urban to water etc.
var merged = beforeClasses.multiply(100).add(afterClasses).rename('transitions')
// Use a frequencyHistogram to get a pixel count per class
var transitionMatrix = merged.reduceRegion({
reducer: ee.Reducer.frequencyHistogram(),
geometry: bangalore,
maxPixels: 1e10,
scale:10,
tileScale: 16
})
// This prints number of pixels for each class transition
print(transitionMatrix.get('transitions')) Exercise
Lost water pixels between 2019 and 2020
// Exercise
// Show all areas where water became other classes and display the result
// Hint1: Select class 3 pixels from before image and NOT class 3 pixels from after image
// Hint2: use the .and() operation to select pixels matching both conditionsModule 5: Earth Engine Apps
This module is focused the concepts related to client vs. server that will help you in creating web apps. We will be building an app using the Earth Engine User Interface API and publishing it to Google Cloud.
01. Client vs. Server
The User Interface elements in your Code Editor - Map View, Drawing Tools etc. are ‘client-side’ elements. They run in YOUR browser. Image Collections, feature collections, calculations on Earth Engine objects etc. are ‘server-side’ elements. They run in Google’s data center. You cannot mix both these objects. To learn more, visit the Client vs. Server section of the Earth Engine User Guide.
- To convert client-side objects to server-side objects, you can use the appropriate API function. Server-side functions start with
ee., suchee.Date(),ee.Image()etc. - To convert server-side objects to client-side objects, you can call
.getInfo()on am Earth Engine object. For the Python API, this is the only way to extract information from a server-side object, but the Javascript API provides a better (and preferred) - method for bring server-side objects to client-side using theevaluate()method. This method asynchronously retrieves the value of the object, without blocking the user interface - meaning it will let your code continue to execute while it fetches the value.
Tip: You can use
ee.Algorithms.ObjectType()to get the type of a server-side object
var date = '2020-01-01' // This is client-side
print(typeof(date))
var eedate = ee.Date('2020-01-01').format() // This is server-side
print(typeof(eedate))
// To bring server-side objects to client-side, you can call .getInfo()
// var clientdate = eedate.getInfo()
// print(clientdate)
// print(typeof(clientdate))
// getInfo() blocks the execution of your code till the value is fetched
// If the value takes time to compute, your code editor will freeze
// This is not a good user experience
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
var filtered = s2.filter(ee.Filter.date('2020-01-01', '2020-01-31'))
//var numImages = filtered.size().getInfo()
//print(numImages)
// A better approach is to use evaluate() function
// You need to define a 'callback' function which will be called once the
// value has been computed and ready to be used.
var myCallback = function(object) {
print(object)
}
print('Computing the size of the collection')
var numImages = filtered.size().evaluate(myCallback)Exercise
var date = ee.Date.fromYMD(2019, 1, 1)
print(date)
// We can use the format() function to create
// a string from a date object
var dateString = date.format('dd MMM, YYYY')
print(dateString)
// Exercise
// The print statement below combines a client-side string
// with a server-side string - resulting in an error.
// Fix the code so that the following message is printed
// 'The date is 01 Jan, 2019'
var message = 'The date is ' + dateString
print(message)
// Hint:
// Convert the client-side string to a server-side string
// Use ee.String() to create a server-side string
// Use the .cat() function instead of + to combine 2 strings02. Building and Publishing an App
Earth Engine comes with a User Interface API that allows you to build an interactive web application powered by Earth Engine. Building a web mapping application typically requires the skills of a full stack developer and are out of reach for most analysts and scientists. But the Earth Engine User Interface API makes this process much more accessible by providing ready-to-use widgets, such as Buttons, Drop-down Menus, Sliders etc. - and free cloud hosting to allow anyone to publish an app with just a few lines of code.
The apps run in your browser, so they need to use client-side functions. All the user interface functions are contained in the ui. package - such as ui.Select(), ui.Button(). You can create those elements by calling these functions with appropriate parameters. The main container object is the ui.Panel() which can contain different types of widgets. Learn more in the Earth Engine User Interface API section of the Earth Engine User Guide.
The code below shows how to build an app called Night Lights Explorer that allows anyone to pick a year/month and load the VIIRS Nighttime Day/Night Band Composite for the selected month. Copy/paste the code below to your Code Editor and click Run.
You will see a panel on the right-hand side with 2 drop-down boxes and a button. These are User Interface (UI) widgets provided by the Earth Engine API that allows the user to interactively select the values. You can select the values for year and month and click Load button to see the image for the selected month. We will now publish this app. Click on the Apps button.
App with UI Elements
In the Manage Apps window, click New App.
Enter the name of your app. The app will be hosted on Google Cloud, so you will need to create and link a Google Cloud project with the app. If you don’t have a Google Cloud account, you can click the Not Yet button to create a new project.
When prompted to Choose a Cloud Project for your apps, you can select Create a new Cloud Project and enter an unique id and click Select.
You may get an error asking you to accept the terms of service. Click the Cloud Terms of Service link and follow the instructions. Once done, click OK.
Back in the Publish New App dialog, leave all other settings to default and click Publish.
The app will be hosted on Google Cloud and you can access it by clicking on the App Name of your app in the Manage Apps dialog.
You will see your Earth Engine powered app running in the browser. Anyone can access and interact with the app by just visiting the App URL.
The app publishing process takes a few minutes. So if you get an error that your app is not yet ready, check back in a few minutes.
// Panels are the main container widgets
var mainPanel = ui.Panel({
style: {width: '300px'}
});
var title = ui.Label({
value: 'Night Lights Explorer',
style: {'fontSize': '24px'}
});
// You can add widgets to the panel
mainPanel.add(title)
// You can even add panels to other panels
var dropdownPanel = ui.Panel({
layout: ui.Panel.Layout.flow('horizontal'),
});
mainPanel.add(dropdownPanel);
var yearSelector = ui.Select({
placeholder: 'please wait..',
})
var monthSelector = ui.Select({
placeholder: 'please wait..',
})
var button = ui.Button('Load')
dropdownPanel.add(yearSelector)
dropdownPanel.add(monthSelector)
dropdownPanel.add(button)
// Let's add a dropdown with the years
var years = ee.List.sequence(2014, 2020)
var months = ee.List.sequence(1, 12)
// Dropdown items need to be strings
var yearStrings = years.map(function(year){
return ee.Number(year).format('%04d')
})
var monthStrings = months.map(function(month){
return ee.Number(month).format('%02d')
})
// Evaluate the results and populate the dropdown
yearStrings.evaluate(function(yearList) {
yearSelector.items().reset(yearList)
yearSelector.setPlaceholder('select a year')
})
monthStrings.evaluate(function(monthList) {
monthSelector.items().reset(monthList)
monthSelector.setPlaceholder('select a month')
})
// Define a function that triggers when any value is changed
var loadComposite = function() {
var col = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG");
var year = yearSelector.getValue()
var month = monthSelector.getValue()
var startDate = ee.Date.fromYMD(
ee.Number.parse(year), ee.Number.parse(month), 1)
var endDate = startDate.advance(1, 'month')
var filtered = col.filter(ee.Filter.date(startDate, endDate))
var image = ee.Image(filtered.first()).select('avg_rad')
var nighttimeVis = {min: 0.0, max: 60.0}
var layerName = 'Night Lights ' + year + '-' + month
Map.addLayer(image, nighttimeVis, layerName)
}
button.onClick(loadComposite)
Map.setCenter(76.43, 12.41, 8)
ui.root.add(mainPanel);Exercise
// Exercise
// Add a button called 'Reset'
// Clicking the button should remove all loaded layers
// Hint: Use Map.clear() for removing the layersModule 6: Google Earth Engine Python API
Introduction to the Python API
Till this point in the course, we have used the Earth Engine Javascript API for all our analysis. Earth Engine also provides a Python API. If you are a Python programmer, you may prefer to use the Python API to integrate Earth Engine in your spatial analysis workflow. There are many options for running Python code that uses the Google Earth Engine API. We will use the following two methods in this course.
Google Colab
An easy way to start using the Google Earth Engine Python API is via Google Colab. Google Colaboratory provides a hosted environment to run Python notebooks without having to install Python locally. It also comes pre-installed with many useful packages - including the Google Earth Engine Python API. You can simply visit https://colab.research.google.com/ and start a new notebook.
Local Development Environment
An advantage of Python API is that you can use it in your own development environment - so you get a lot more flexibility to automate as well as combine other analysis and visualization libraries with Earth Engine. This requires installing Python and the Earth Engine Python API on your machine or server. You also need to do a one-time authentication and save the token on the machine. The preferred method for installing the Earth Engine Python API is via Anaconda. Please follow our Google Earth Engine Python API Installation Guide for step-by-step instructions.
01. Python API Syntax
Coming from the programming in Earth Engine through the Code Editor, you will need to slightly adapt your scripts to be able to run in Python. For the bulk of your code, you will be using Earth Engine API’s server-side objects and functions - which will be exactly the same in Python. You only need to make a few syntactical changes.
Here’s the full list of differences.
Initialization
First of all, you need to run the following cells to initialize the API and authorize your account. You will be prompted to sign-in to the account and allow access to View and Manage your Earth Engine data. Once you approve, you will get an a verification code that needs to be entered at the prompt. This step needs to be done just once per session.
import eeee.Authenticate()ee.Initialize()Variables
Python code doesn’t use the ‘var’ keyword
var city = 'San Fransico'
var state = 'California'
print(city, state)
var population = 881549
print(population)city = 'San Fransico'
state = 'California'
print(city, state)
population = 881549
print(population)Line Continuation
Python doesn’t use a semi-colon for line ending. To indicate line-continuation you need to use the \ character
var s2 = ee.ImageCollection('COPERNICUS/S2');
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'));s2 = ee.ImageCollection('COPERNICUS/S2')
filtered = s2 \
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \
.filter(ee.Filter.date('2019-01-01', '2019-12-31')) Functions
Instead of the function keyword, Python uses the def keyword. Also the in-line functions are defined using lambda anonymous functions.
In the example below, also now the and operator - which is capitalized as And in Python version to avoid conflict with the built-in and operator. The same applies to Or and Not operators. true, false, null in Python are also spelled as True, False and None.
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)//.divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}def maskS2clouds(image):
qa = image.select('QA60')
cloudBitMask = 1 << 10
cirrusBitMask = 1 << 11
mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask) \
.select("B.*") \
.copyProperties(image, ["system:time_start"])Function Arguments
Named arguments to Earth Engine functions need to be in quotes. Also when passing the named arguments as a dictionary, it needs to be passed using the ** keyword.
var gaul = ee.FeatureCollection("FAO/GAUL/2015/level1");
var gfsad = ee.Image("USGS/GFSAD1000_V0");
var wheatrice = gfsad.select('landcover').eq(1)
var uttarpradesh = gaul.filter(ee.Filter.eq('ADM1_NAME', 'Uttar Pradesh'))
var points = wheatrice.selfMask().stratifiedSample(
{numPoints:100, region:uttarpradesh, geometries: true})gaul = ee.FeatureCollection("FAO/GAUL/2015/level1")
gfsad = ee.Image("USGS/GFSAD1000_V0")
wheatrice = gfsad.select('landcover').eq(1)
uttarpradesh = gaul.filter(ee.Filter.eq('ADM1_NAME', 'Uttar Pradesh'))
points = wheatrice.selfMask().stratifiedSample(**
{'numPoints':100, 'region':uttarpradesh, 'geometries': True})
print(points.getInfo())In-line functions
The syntax for defining in-line functions is also slightly different. You need to use the lambda keyword
var points = points.map(function(feature) {
return ee.Feature(feature.geometry(), {'id': feature.id()})
})points = points.map(lambda feature: ee.Feature(feature.geometry(), {'id': feature.id()} ))Printing Values
The print() function syntax is the same. But you must remember that in the Code Editor when you cann print, the value of the server object is fetched and then printed. You must do that explicitely by calling getInfo() on any server-side object.
print(points.first()print(points.first().getInfo())Exercise
Take the Javascript code snippet below and write the equiavalent Python code in the cell below.
- Hint1: Chaining of filters require the use of line continuation character
\ - Hint2: Printing of server-side objects requires calling
.getInfo()on the object
The correct code should print the value 30.
var geometry = ee.Geometry.Point([77.60412933051538, 12.952912912328241]);
var s2 = ee.ImageCollection('COPERNICUS/S2');
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.bounds(geometry));
print(filtered.size());import ee
ee.Authenticate()
ee.Initialize()02. Automatic Conversion of Javascript Code to Python
Interactive leaflet map created by geemap
geemap is an open-source Python package that comes with many helpful features that help you use Earth Engine effectively in Python.
It comes with a function that can help you translate your javascript earth engine code to Python automatically.
Google Colab doesn’t come pre-installed with the package, so we install it via pip.
try:
import geemap
except ModuleNotFoundError:
if 'google.colab' in str(get_ipython()):
print('geemap not found, installing via pip in Google Colab...')
!pip install geemap --quiet
import geemap
else:
print('geemap not found, please install via conda in your environment')The automatic conversion of code is done by calling the geemap.js_snippet_to_py() function. We first create a string with the javascript code.
javascript_code = """
var geometry = ee.Geometry.Point([107.61303468448624, 12.130969369851766]);
Map.centerObject(geometry, 12)
var s2 = ee.ImageCollection("COPERNICUS/S2")
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = s2
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(geometry))
.map(maskS2clouds)
// Write a function that computes NDVI for an image and adds it as a band
function addNDVI(image) {
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('ndvi');
return image.addBands(ndvi);
}
var withNdvi = filtered.map(addNDVI);
var composite = withNdvi.median()
palette = [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
'74A901', '66A000', '529400', '3E8601', '207401', '056201',
'004C00', '023B01', '012E01', '011D01', '011301'];
ndviVis = {min:0, max:0.5, palette: palette }
Map.addLayer(withNdvi.select('ndvi'), ndviVis, 'NDVI Composite')
"""lines = geemap.js_snippet_to_py(
javascript_code, add_new_cell=False,
import_ee=True, import_geemap=True, show_map=True)
for line in lines:
print(line.rstrip())The automatic conversion works great. Review it and paste it to the cell below.
import ee
import geemap
Map = geemap.Map()
geometry = ee.Geometry.Point([107.61303468448624, 12.130969369851766])
Map.centerObject(geometry, 12)
s2 = ee.ImageCollection("COPERNICUS/S2")
rgbVis = {
'min': 0.0,
'max': 3000,
'bands': ['B4', 'B3', 'B2'],
}
# Write a function for Cloud masking
def maskS2clouds(image):
qa = image.select('QA60')
cloudBitMask = 1 << 10
cirrusBitMask = 1 << 11
mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask) \
.select("B.*") \
.copyProperties(image, ["system:time_start"])
filtered = s2 \
.filter(ee.Filter.date('2019-01-01', '2019-12-31')) \
.filter(ee.Filter.bounds(geometry)) \
.map(maskS2clouds)
# Write a function that computes NDVI for an image and adds it as a band
def addNDVI(image):
ndvi = image.normalizedDifference(['B5', 'B4']).rename('ndvi')
return image.addBands(ndvi)
withNdvi = filtered.map(addNDVI)
composite = withNdvi.median()
palette = [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
'74A901', '66A000', '529400', '3E8601', '207401', '056201',
'004C00', '023B01', '012E01', '011D01', '011301']
ndviVis = {'min':0, 'max':0.5, 'palette': palette }
Map.addLayer(withNdvi.select('ndvi'), ndviVis, 'NDVI Composite')
MapExercise
Take the Javascript code snippet below and use geemap to automatically convert it to Python.
var admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");
var karnataka = admin2.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
var visParams = {color: 'red'}
Map.centerObject(karnataka)
Map.addLayer(karnataka, visParams, 'Karnataka Districts')03. Batch Exports
One of the most commonly asked questions by Earth Engine users is - How do I download all images in a collection? The Earth Engine Python API comes with a ee.batch module that allows you to launch batch exports and manage tasks. The recommended way to do batch exports like this is to use the Python API’s ee.batch.Export functions and use a Python for-loop to iterate and export each image. The ee.batch module also gives you ability to control Tasks - allowing you to automate exports.
import eeee.Authenticate()ee.Initialize()Create a Collection
geometry = ee.Geometry.Point([107.61303468448624, 12.130969369851766])
s2 = ee.ImageCollection("COPERNICUS/S2")
rgbVis = {
'min': 0.0,
'max': 3000,
'bands': ['B4', 'B3', 'B2'],
}
# Write a function for Cloud masking
def maskS2clouds(image):
qa = image.select('QA60')
cloudBitMask = 1 << 10
cirrusBitMask = 1 << 11
mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask) \
.select("B.*") \
.copyProperties(image, ["system:time_start"])
filtered = s2 \
.filter(ee.Filter.date('2019-01-01', '2020-01-01')) \
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) \
.filter(ee.Filter.bounds(geometry)) \
.map(maskS2clouds)
# Write a function that computes NDVI for an image and adds it as a band
def addNDVI(image):
ndvi = image.normalizedDifference(['B5', 'B4']).rename('ndvi')
return image.addBands(ndvi)
withNdvi = filtered.map(addNDVI)Export All Images
Exports are done via the ee.batch module. A key difference between javascript and Python version is that the region parameter needs to be supplied with actual geometry coordinates.
image_ids = withNdvi.aggregate_array('system:index').getInfo()
print('Total images: ', len(image_ids))# Export with 100m resolution for this demo
for i, image_id in enumerate(image_ids):
image = ee.Image(withNdvi.filter(ee.Filter.eq('system:index', image_id)).first())
task = ee.batch.Export.image.toDrive(**{
'image': image.select('ndvi'),
'description': 'Image Export {}'.format(i+1),
'fileNamePrefix': image.id().getInfo(),
'folder':'earthengine',
'scale': 100,
'region': image.geometry().bounds().getInfo()['coordinates'],
'maxPixels': 1e10
})
task.start()
print('Started Task: ', i+1)Manage Running/Waiting Tasks
You can manage tasks as well. Get a list of tasks and get state information on them
tasks = ee.batch.Task.list()
for task in tasks:
task_id = task.status()['id']
task_state = task.status()['state']
print(task_id, task_state)You can cancel tasks as well
tasks = ee.batch.Task.list()
for task in tasks:
task_id = task.status()['id']
task_state = task.status()['state']
if task_state == 'RUNNING' or task_state == 'READY':
task.cancel()
print('Task {} canceled'.format(task_id))
else:
print('Task {} state is {}'.format(task_id, task_state))Exercise
The code below uses the TerraClimate data and creates an ImageCollection with 12 monthly images of maximum temperature. It also extract the geometry for Australia from the LSIB collection. Add the code to start an export task for each image in the collection for australia.
- Hint1: TerraClimate images have a scale of 4638.3m
import ee
lsib = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
australia = lsib.filter(ee.Filter.eq('country_na', 'Australia'))
geometry = australia.geometry()
terraclimate = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')
tmax = terraclimate.select('tmmx')
def scale(image):
return image.multiply(0.1) \
.copyProperties(image,['system:time_start'])
tmaxScaled = tmax.map(scale)
filtered = tmaxScaled \
.filter(ee.Filter.date('2020-01-01', '2021-01-01')) \
.filter(ee.Filter.bounds(geometry))
image_ids = filtered.aggregate_array('system:index').getInfo()
print('Total images: ', len(image_ids))Replace the comments with your code.
for i, image_id in enumerate(image_ids):
exportImage = ee.Image(filtered.filter(ee.Filter.eq('system:index', image_id)).first())
# Clip the image to the region geometry
clippedImage = exportImage.clip(geometry)
## Create the export task using ee.batch.Export.image.toDrive()
## Start the task
Launching multiple tasks using the Python API
04. Creating Charts in Python
The Google Earth Engine Python API does not come with a charting module. But you can use third-party modules to create interactive charts. You may also convert the Earth Engine objects to a Pandas dataframe and plot it using Python libraries like Matplotlib
This notebook shows how to use the geemap package to create a Time-Series Chart from a ImageCollection.
References:
- geemap Chart module
- geemap Example notebook
import eetry:
import geemap
except ModuleNotFoundError:
if 'google.colab' in str(get_ipython()):
print('geemap not found, installing via pip in Google Colab...')
!pip install geemap --quiet
import geemap
else:
print('geemap not found, please install via conda in your environment')ee.Authenticate()ee.Initialize()Load the TerraClimate collection and select the ‘tmmx’ band.
terraclimate = ee.ImageCollection("IDAHO_EPSCOR/TERRACLIMATE")
tmax = terraclimate.select('tmmx')Define a point location for the chart.
geometry = ee.Geometry.Point([77.57738128916243, 12.964758918835752])Scale the band values so they are in Degree Celcius.
def scale_image(image):
return ee.Image(image).multiply(0.1)\
.copyProperties(image, ['system:time_start'])
tmaxScaled = tmax.map(scale_image)Filter the collection.
filtered = tmaxScaled.filter(ee.Filter.date('2019-01-01', '2020-01-01')) \
.filter(ee.Filter.intersects('.geo', geometry))To chart an image series in Python, we must first extract the values from each image and create a FeatureCollection.
def extract_data(image):
stats = image.reduceRegion(**{
'reducer':ee.Reducer.mean(),
'geometry':geometry,
'scale':5000
})
properties = {
'month': image.get('system:index'),
'tmmx': stats.get('tmmx')
}
return ee.Feature(None, properties)
data = ee.FeatureCollection(filtered.map(extract_data))print(data.first().getInfo())Create an Interactive Chart using geemap
from geemap import chartoptions = {
'title': 'Max Monthly Temperature at Bangalore',
'legend_location': 'top-right',
'height': '500px',
'ylabel': 'Temperature (C)',
'xlabel': 'Date'
}chart.feature_byFeature(data, 'month', ['tmmx'], **options)Create a chart using Matplotlib
We can convert a FeatureCollection to a DataFrame using geemap helper function ee_to_pandas.
import geemap
df = geemap.ee_to_pandas(data)Now we have a regular Pandas dataframe that can be plotted with matplotlib.
%matplotlib inline
import matplotlib.pyplot as pltfig, ax = plt.subplots()
fig.set_size_inches(20,10)
df.plot(ax=ax,
title='Max Monthly Temperature at Bangalore',
x='month',
ylabel='Temperature (C)',
kind='bar')
plt.tight_layout()Exercise
Customize the above chart by plotting it as a Line Chart in red color.
- Hint1: Use
kind='line'along with acoloroption.
05. Automating Downloads
Another common use of the GEE Python API is to automate data processing and export. You can create a Python script that can be called from a server or launched on a schedule using tools such as Windows Scheduler or crontab.
This script below provides a complete example of automating a download using Google Earth Engine API. It uses the Google Earth Engine API to compute the average soil moisture for the past 2 weeks over all districts in a state. The result is then downloaded as a JSON file and saved locally.
Make sure you have completed the one-time authentication flow before running the script.
Follow the steps below to create a script to download data from GEE.
- Create a new file named
download_data.pywith the content shown below.
import datetime
import ee
import csv
import os
ee.Initialize()
# Get current date and convert to milliseconds
end_date = ee.Date(datetime.datetime.now().timestamp()*1000)
start_date = end_date.advance(-2, 'week')
date_string = end_date.format('YYYY_MM_dd')
filename = 'ssm_{}.csv'.format(date_string.getInfo())
# Saving to current directory. You can change the path to appropriate location
output_path = os.path.join(filename)
# Datasets
soilmoisture = ee.ImageCollection("NASA_USDA/HSL/SMAP10KM_soil_moisture")
admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2")
# Filter to a state
karnataka = admin2.filter(ee.Filter.eq('ADM1_NAME', 'Karnataka'))
# Select the ssm band
ssm = soilmoisture.select('ssm')
filtered = ssm .filter(ee.Filter.date(start_date, end_date))
mean = filtered.mean()
stats = mean.reduceRegions(**{
'collection': karnataka,
'reducer': ee.Reducer.mean().setOutputs(['meanssm']),
'scale': 10000,
'crs': 'EPSG:32643'
})
# Select columns to keep and remove geometry to make the result lightweight
# Change column names to match your uploaded shapefile
columns = ['ADM2_NAME', 'meanssm']
exportCollection = stats.select(**{
'propertySelectors': columns,
'retainGeometry': False})
features = exportCollection.getInfo()['features']
data = []
for f in features:
data.append(f['properties'])
field_names = ['ADM2_NAME', 'meanssm']
with open(output_path, 'w') as csvfile:
writer = csv.DictWriter(csvfile, fieldnames = field_names)
writer.writeheader()
writer.writerows(data)
print('Success: File written at', output_path)- From the terminal, navigate to the directory where you have created the file and type the command below to run the script.
python download_data.py
- The script will download the data from GEE and save a file to your current directory.
Supplement
This section contains useful scripts and code snippets that can be adapted for your projects.
Advanced Supervised Classification Techniques
Hyperparameter Tuning
A recommended best practice for improving the accuracy of your machine learning model is to tune different parameters. For example, when using the ee.Classifier.smileRandomForest() classifier, we must specify the Number of Trees. We know that higher number of trees result in more computation requirement, but it doesn’t necessarily result in better results. Instead of guessing, we programmatically try a range of values and choose the smallest value possible that results in the highest accuracy.
Supervised Classification Output
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7");
var gcp = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_gcps");
var alos = ee.Image("JAXA/ALOS/AW3D30/V2_2");
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
// Function to remove cloud and snow pixels from Sentinel-2 SR image
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(10);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 10% or cloud shadow classification
var mask = cloud.and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask);
}
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.map(maskCloudAndShadowsSR)
.select('B.*')
var composite = filtered.median().clip(boundary)
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000, gamma: 1.2};
Map.centerObject(boundary)
Map.addLayer(composite, visParams, 'RGB');
var addIndices = function(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
var ndbi = image.normalizedDifference(['B11', 'B8']).rename(['ndbi']);
var mndwi = image.normalizedDifference(['B3', 'B11']).rename(['mndwi']);
var bsi = image.expression(
'(( X + Y ) - (A + B)) /(( X + Y ) + (A + B)) ', {
'X': image.select('B11'), //swir1
'Y': image.select('B4'), //red
'A': image.select('B8'), // nir
'B': image.select('B2'), // blue
}).rename('bsi');
return image.addBands(ndvi).addBands(ndbi).addBands(mndwi).addBands(bsi)
}
var composite = addIndices(composite);
// Calculate Slope and Elevation
var elev = alos.select('AVE_DSM').rename('elev');
var slope = ee.Terrain.slope(alos.select('AVE_DSM')).rename('slope');
var composite = composite.addBands(elev).addBands(slope);
// Normalize the image
// Machine learning algorithms work best on images when all features have
// the same range
// Function to Normalize Image
// Pixel Values should be between 0 and 1
// Formula is (x - xmin) / (xmax - xmin)
//**************************************************************************
function normalize(image){
var bandNames = image.bandNames();
// Compute min and max of the image
var minDict = image.reduceRegion({
reducer: ee.Reducer.min(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var maxDict = image.reduceRegion({
reducer: ee.Reducer.max(),
geometry: boundary,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var mins = ee.Image.constant(minDict.values(bandNames));
var maxs = ee.Image.constant(maxDict.values(bandNames));
var normalized = image.subtract(mins).divide(maxs.subtract(mins))
return normalized
}
var composite = normalize(composite);
// Add a random column and split the GCPs into training and validation set
var gcp = gcp.randomColumn()
// This being a simpler classification, we take 60% points
// for validation. Normal recommended ratio is
// 70% training, 30% validation
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
Map.addLayer(validationGcp)
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: trainingGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
print(training)
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
//**************************************************************************
// Feature Importance
//**************************************************************************
// Run .explain() to see what the classifer looks like
print(classifier.explain())
// Calculate variable importance
var importance = ee.Dictionary(classifier.explain().get('importance'))
// Calculate relative importance
var sum = importance.values().reduce(ee.Reducer.sum())
var relativeImportance = importance.map(function(key, val) {
return (ee.Number(val).multiply(100)).divide(sum)
})
print(relativeImportance)
// Create a FeatureCollection so we can chart it
var importanceFc = ee.FeatureCollection([
ee.Feature(null, relativeImportance)
])
var chart = ui.Chart.feature.byProperty({
features: importanceFc
}).setOptions({
title: 'Feature Importance',
vAxis: {title: 'Importance'},
hAxis: {title: 'Feature'}
})
print(chart)
//**************************************************************************
// Hyperparameter Tuning
//**************************************************************************
var test = composite.sampleRegions({
collection: validationGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16
});
// Tune the numberOfTrees parameter.
var numTreesList = ee.List.sequence(10, 150, 10);
var accuracies = numTreesList.map(function(numTrees) {
var classifier = ee.Classifier.smileRandomForest(numTrees)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Here we are classifying a table instead of an image
// Classifiers work on both images and tables
return test
.classify(classifier)
.errorMatrix('landcover', 'classification')
.accuracy();
});
var chart = ui.Chart.array.values({
array: ee.Array(accuracies),
axis: 0,
xLabels: numTreesList
}).setOptions({
title: 'Hyperparameter Tuning for the numberOfTrees Parameters',
vAxis: {title: 'Validation Accuracy'},
hAxis: {title: 'Number of Tress', gridlines: {count: 15}}
});
print(chart)
// Tuning Multiple Parameters
// We can tune many parameters together using
// nested map() functions
// Let's tune 2 parameters
// numTrees and bagFraction
var numTreesList = ee.List.sequence(10, 150, 10);
var bagFractionList = ee.List.sequence(0.1, 0.9, 0.1);
var accuracies = numTreesList.map(function(numTrees) {
return bagFractionList.map(function(bagFraction) {
var classifier = ee.Classifier.smileRandomForest({
numberOfTrees: numTrees,
bagFraction: bagFraction
})
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Here we are classifying a table instead of an image
// Classifiers work on both images and tables
var accuracy = test
.classify(classifier)
.errorMatrix('landcover', 'classification')
.accuracy();
return ee.Feature(null, {'accuracy': accuracy,
'numberOfTrees': numTrees,
'bagFraction': bagFraction})
})
}).flatten()
var resultFc = ee.FeatureCollection(accuracies)
// Export the result as CSV
Export.table.toDrive({
collection: resultFc,
description: 'Multiple_Parameter_Tuning_Results',
folder: 'earthengine',
fileNamePrefix: 'numtrees_bagfraction',
fileFormat: 'CSV'})
Post-Processing Classification Results
Supervised classification results often contain salt-and-pepper noise caused by mis-classified pixels. It is usually preferable to apply some post-processing techniques to remove such noise. The following script contains the code for two popular techniques for post-processing classification results.
- Using un-supervised clustering to replacing classified value by majority value in each cluster.
- Replacing isolated pixels with surrounding value with a majority filter.
Remember that the neighborhood methods are scale-dependent so the results will change as you zoom in/out. Export the results at the desired scale to see the effect of post-processing.
// Sentinel-2 Median Composite
var composite = ee.Image("users/ujavalgandhi/e2e/arkavathy_2019_composite");
Map.addLayer(composite, {min: 0, max: 0.3, bands: ['B4', 'B3', 'B2']}, 'RGB Composite');
// Raw Supervised Classification Results
var classified = ee.Image("users/ujavalgandhi/e2e/arkavathy_final_classification");
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, 'Original');
Map.centerObject(classified, 10)
//**************************************************************************
// Post process by clustering
//**************************************************************************
// Cluster using Unsupervised Clustering methods
var seeds = ee.Algorithms.Image.Segmentation.seedGrid(5);
var snic = ee.Algorithms.Image.Segmentation.SNIC({
image: composite.select('B.*'),
compactness: 0,
connectivity: 4,
neighborhoodSize: 10,
size: 2,
seeds: seeds
})
var clusters = snic.select('clusters')
// Assign class to each cluster based on 'majority' voting (using ee.Reducer.mode()
var smoothed = classified.addBands(clusters);
var clusterMajority = smoothed.reduceConnectedComponents({
reducer: ee.Reducer.mode(),
labelBand: 'clusters'
});
Map.addLayer(clusterMajority, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']},
'Processed using Clusters');
//**************************************************************************
// Post process by replacing isolated pixels with surrounding value
//**************************************************************************
// count patch sizes
var patchsize = classified.connectedPixelCount(40, false);
// run a majority filter
var filtered = classified.focal_mode({
radius: 10,
kernelType: 'square',
units: 'meters',
});
// updated image with majority filter where patch size is small
var connectedClassified = classified.where(patchsize.lt(25),filtered);
Map.addLayer(connectedClassified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']},
'Processed using Connected Pixels');Principal Component Analysis (PCA)
PCA is a very useful technique in improving your supervised classification results. This is a statistical technique that compresses data from a large number of bands into fewer uncorrelated bands. You can run PCA on your image and add the first few (typically 3) principal component bands to the original composite before sampling training points. In the example below, you will notice that 97% of the variance from the 13-band original image is captured in the 3-band PCA image. This sends a stronger signal to the classifier and improves accuracy by allowing it to distinguish different classes better.
var composite = ee.Image("users/ujavalgandhi/e2e/arkavathy_2019_composite");
var boundary = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_boundary")
Map.addLayer(composite, {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3, gamma: 1.2}, 'RGB');
//**************************************************************************
// Function to calculate Principal Components
// Code adapted from https://developers.google.com/earth-engine/guides/arrays_eigen_analysis
//**************************************************************************
function PCA(maskedImage){
var image = maskedImage.unmask()
var scale = 20;
var region = boundary;
var bandNames = image.bandNames();
// Mean center the data to enable a faster covariance reducer
// and an SD stretch of the principal components.
var meanDict = image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: region,
scale: scale,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);
// This helper function returns a list of new band names.
var getNewBandNames = function(prefix) {
var seq = ee.List.sequence(1, bandNames.length());
return seq.map(function(b) {
return ee.String(prefix).cat(ee.Number(b).int());
});
};
// This function accepts mean centered imagery, a scale and
// a region in which to perform the analysis. It returns the
// Principal Components (PC) in the region as a new image.
var getPrincipalComponents = function(centered, scale, region) {
// Collapse the bands of the image into a 1D array per pixel.
var arrays = centered.toArray();
// Compute the covariance of the bands within the region.
var covar = arrays.reduceRegion({
reducer: ee.Reducer.centeredCovariance(),
geometry: region,
scale: scale,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
// Get the 'array' covariance result and cast to an array.
// This represents the band-to-band covariance within the region.
var covarArray = ee.Array(covar.get('array'));
// Perform an eigen analysis and slice apart the values and vectors.
var eigens = covarArray.eigen();
// This is a P-length vector of Eigenvalues.
var eigenValues = eigens.slice(1, 0, 1);
// Compute Percentage Variance of each component
var eigenValuesList = eigenValues.toList().flatten()
var total = eigenValuesList.reduce(ee.Reducer.sum())
var percentageVariance = eigenValuesList.map(function(item) {
return (ee.Number(item).divide(total)).multiply(100).format('%.2f')
})
// This will allow us to decide how many components capture
// most of the variance in the input
print('Percentage Variance of Each Component', percentageVariance)
// This is a PxP matrix with eigenvectors in rows.
var eigenVectors = eigens.slice(1, 1);
// Convert the array image to 2D arrays for matrix computations.
var arrayImage = arrays.toArray(1);
// Left multiply the image array by the matrix of eigenvectors.
var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
// Turn the square roots of the Eigenvalues into a P-band image.
var sdImage = ee.Image(eigenValues.sqrt())
.arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
// Turn the PCs into a P-band image, normalized by SD.
return principalComponents
// Throw out an an unneeded dimension, [[]] -> [].
.arrayProject([0])
// Make the one band array image a multi-band image, [] -> image.
.arrayFlatten([getNewBandNames('pc')])
// Normalize the PCs by their SDs.
.divide(sdImage);
};
var pcImage = getPrincipalComponents(centered, scale, region);
return pcImage.mask(maskedImage.mask());
}
// As you see from the printed results, ~97% of the variance
// from the original image is captured in the first 3 principal components
// We select those and discard others
var pca = PCA(composite).select(['pc1', 'pc2', 'pc3'])
// Add the 3 PCA bands to the original image
// Use this image for sampling training points for
// supervised classification
var composite = composite.addBands(pca)
// PCA computation is expensive and can time out when displaying on the map
// Export the results and import them back
Export.image.toAsset({
image: pca,
description: 'Principal_Components_Image',
assetId: 'users/ujavalgandhi/e2e/arkavathy_pca',
region: boundary,
scale: 20,
maxPixels: 1e10})
// Once the export finishes, import the asset and display
var pcaImported = ee.Image('users/ujavalgandhi/e2e/arkavathy_pca')
var pcaVisParams = {bands: ['pc1', 'pc2', 'pc3'], min: -2, max: 2};
Map.addLayer(pcaImported, pcaVisParams, 'Principal Components');Multi-temporal Composites for Crop Classification
Crop classification is a difficult problem. A useful technique that aids in clear distinction of crops is to account for crop phenology. This technique can be applied to detect a specific type of crop or distinguish crops from other forms of vegetation. You can create composite images for different periods of the cropping cycle and create a stacked image to be used for classification. This allows the classifier to learn the temporal pattern and detect pixels that exhibit similar patterns.
Capturing Crop Phenology through Seasonal Composites
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7")
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
Map.centerObject(boundary, 11)
// Function to remove cloud pixels from Sentinel-2 SR image
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(10);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 10% or cloud shadow classification
var mask = cloud.and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask).divide(10000)
.copyProperties(image, ['system:time_start']);
}
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.map(maskCloudAndShadowsSR)
// There are 3 distinct crop seasons in the area of interest
// Jan-March = Winter (Rabi) Crops
// April-June = Summer Crops / Harvest
// July-December = Monsoon (Kharif) Crops
var cropCalendar = ee.List([[1,3], [4,6], [7,12]])
// We create different composites for each season
var createSeasonComposites = function(months) {
var startMonth = ee.List(months).get(0)
var endMonth = ee.List(months).get(1)
var monthFilter = ee.Filter.calendarRange(startMonth, endMonth, 'month')
var seasonFiltered = filtered.filter(monthFilter)
var composite = seasonFiltered.median()
return composite.select('B.*').clip(boundary)
}
var compositeList = cropCalendar.map(createSeasonComposites)
var rabi = ee.Image(compositeList.get(0))
var harvest = ee.Image(compositeList.get(1))
var kharif = ee.Image(compositeList.get(2))
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3, gamma: 1.2};
Map.addLayer(rabi, visParams, 'Rabi')
Map.addLayer(harvest, visParams, 'Harvest')
Map.addLayer(kharif, visParams, 'Kharif')
// Create a stacked image with composites from all seasons
// This multi-temporal image is able capture the crop phenology
// Classifier will be able to detect crop-pixels from non-crop pixels
var composite = rabi.addBands(harvest).addBands(kharif)
// This is a 36-band image
// Use this image for sampling training points for
// to train a crop classifier
print(composite)Calculating Area by Class
This code snippet shows how to use a Grouped Reducer to calculate area covered by each class in a classified image. It also shows how to use the ui.Chart.image.byClass() function to create a chart showing the area for each class.
var classified = ee.Image("users/ujavalgandhi/e2e/bangalore_classified");
var bangalore = ee.FeatureCollection("users/ujavalgandhi/public/bangalore_boundary");
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
// Create a 2 band image with the area image and the classified image
// Divide the area image by 1e6 so area results are in Sq Km
var areaImage = ee.Image.pixelArea().divide(1e6).addBands(classified);
// Calculate Area by Class
// Using a Grouped Reducer
var areas = areaImage.reduceRegion({
reducer: ee.Reducer.sum().group({
groupField: 1,
groupName: 'classification',
}),
geometry: bangalore,
scale: 100,
tileScale: 4,
maxPixels: 1e10
});
var classAreas = ee.List(areas.get('groups'))
print(classAreas)
var areaChart = ui.Chart.image.byClass({
image: areaImage,
classBand: 'classification',
region: bangalore,
scale: 100,
reducer: ee.Reducer.sum(),
classLabels: ['urban', 'bare', 'water', 'vegetation'],
}).setOptions({
hAxis: {title: 'Classes'},
vAxis: {title: 'Area Km^2'},
title: 'Area by class',
series: {
0: { color: 'gray' },
1: { color: 'brown' },
2: { color: 'blue' },
3: { color: 'green' }
}
});
print(areaChart); Spectral Signature Plots
For supervised classification, it is useful to visualize average spectral responses for each band for each class. Such charts are called Spectral Response Curves or Spectral Signatures. Such charts helps determine separability of classes. If classes have very different signatures, a classifier will be able to separate them well.
We can also plot spectral signatures of all training samples for a class and check the quality of the training dataset. If all training samples show similar signatures - it indicates that you have done a good job of collecting appropriate samples. You can also catch potential outliers from these plots.
These charts provide a qualitative and visual methods for checking separability of classes. For quantitative methods, one can apply measures such as Spectral Distance, Mahalanobis distance, Bhattacharyya distance , Jeffreys-Matusita (JM) distance etc. You can find the code for these in this Stack Exchange answer.
Mean Signatures for All Classes
Spectral Signatures for All Training Points by Class
var gcps = ee.FeatureCollection("users/ujavalgandhi/e2e/bangalore_gcps");
var composite = ee.Image('users/ujavalgandhi/e2e/bangalore_composite');
// Overlay the point on the image to get bands data.
var training = composite.sampleRegions({
collection: gcps,
properties: ['landcover'],
scale: 10
});
// We will create a chart of spectral signature for all classes
// We have multiple GCPs for each class
// Use a grouped reducer to calculate the average reflectance
// for each band for each class
// We have 12 bands so need to repeat the reducer 12 times
// We also need to group the results by class
// So we find the index of the landcover property and use it
// to group the results
var bands = composite.bandNames()
var numBands = bands.length()
var bandsWithClass = bands.add('landcover')
var classIndex = bandsWithClass.indexOf('landcover')
// Use .combine() to get a reducer capable of
// computing multiple stats on the input
var combinedReducer = ee.Reducer.mean().combine({
reducer2: ee.Reducer.stdDev(),
sharedInputs: true})
// Use .repeat() to get a reducer for each band
// We then use .group() to get stats by class
var repeatedReducer = combinedReducer.repeat(numBands).group(classIndex)
var gcpStats = training.reduceColumns({
selectors: bands.add('landcover'),
reducer: repeatedReducer,
})
// Result is a dictionary, we do some post-processing to
// extract the results
var groups = ee.List(gcpStats.get('groups'))
var classNames = ee.List(['urban', 'bare', 'water', 'vegetation'])
var fc = ee.FeatureCollection(groups.map(function(item) {
// Extract the means
var values = ee.Dictionary(item).get('mean')
var groupNumber = ee.Dictionary(item).get('group')
var properties = ee.Dictionary.fromLists(bands, values)
var withClass = properties.set('class', classNames.get(groupNumber))
return ee.Feature(null, withClass)
}))
// Chart spectral signatures of training data
var options = {
title: 'Average Spectral Signatures',
hAxis: {title: 'Bands'},
vAxis: {title: 'Reflectance',
viewWindowMode:'explicit',
viewWindow: {
max:0.6,
min:0
}},
lineWidth: 1,
pointSize: 4,
series: {
0: {color: 'grey'},
1: {color: 'brown'},
2: {color: 'blue'},
3: {color: 'green'},
}};
// Default band names don't sort propertly
// Instead, we can give a dictionary with
// labels for each band in the X-Axis
var bandDescriptions = {
'B2': 'B02/Blue',
'B3': 'B03/Green',
'B4': 'B04/Red',
'B8': 'B08/NIR',
'B11': 'B11/SWIR-1',
'B12': 'B12/SWIR-2'
}
// Create the chart and set options.
var chart = ui.Chart.feature.byProperty({
features: fc,
xProperties: bandDescriptions,
seriesProperty: 'class'
})
.setChartType('ScatterChart')
.setOptions(options);
print(chart)
var classChart = function(landcover, label, color) {
var options = {
title: 'Spectral Signatures for ' + label + ' Class',
hAxis: {title: 'Bands'},
vAxis: {title: 'Reflectance',
viewWindowMode:'explicit',
viewWindow: {
max:0.6,
min:0
}},
lineWidth: 1,
pointSize: 4,
};
var fc = training.filter(ee.Filter.eq('landcover', landcover))
var chart = ui.Chart.feature.byProperty({
features: fc,
xProperties: bandDescriptions,
})
.setChartType('ScatterChart')
.setOptions(options);
print(chart)
}
classChart(0, 'Urban')
classChart(1, 'Bare')
classChart(2, 'Water')
classChart(3, 'Vegetation')Identify Misclassified GCPs
While doing accuracy assessment, you will see the validation features that were not classified correctly. It is useful to visually see the points that were misclassified. We can use ee.Filter.eq() and ee.Filter.neq() filters to filter the features where the actual and predicted classes were different. The code below shows how to implement this and also use the style() function visualize them effectively.
// Script that shows how to apply filters to identify
// validation points that were misclassified
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var basin = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_7");
var gcp = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_gcps");
Map.centerObject(gcp)
var arkavathy = basin.filter(ee.Filter.eq('HYBAS_ID', 4071139640))
var boundary = arkavathy.geometry()
var rgbVis = {
min: 0.0,
max: 3000,
bands: ['B4', 'B3', 'B2'],
};
var filtered = s2
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(boundary))
.select('B.*')
var composite = filtered.median().clip(boundary)
// Display the input composite.
Map.addLayer(composite, rgbVis, 'image');
// Add a random column and split the GCPs into training and validation set
var gcp = gcp.randomColumn()
var trainingGcp = gcp.filter(ee.Filter.lt('random', 0.6));
var validationGcp = gcp.filter(ee.Filter.gte('random', 0.6));
// Overlay the point on the image to get training data.
var training = composite.sampleRegions({
collection: trainingGcp,
properties: ['landcover'],
scale: 10,
tileScale: 16,
geometries:true
});
// Train a classifier.
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'landcover',
inputProperties: composite.bandNames()
});
// Classify the image.
var classified = composite.classify(classifier);
Map.addLayer(classified, {min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
var test = classified.sampleRegions({
collection: validationGcp,
properties: ['landcover'],
tileScale: 16,
scale: 10,
geometries:true
});
var testConfusionMatrix = test.errorMatrix('landcover', 'classification')
print('Confusion Matrix', testConfusionMatrix);
// Let's apply filters to find misclassified points
// We can find all points which are labeled landcover=0 (urban)
// but were not classified correctly
// We use ee.Filter.and() function to create a combined filter
var combinedFilter = ee.Filter.and(
ee.Filter.eq('landcover', 0), ee.Filter.neq('classification', 0))
var urbanMisclassified = test.filter(combinedFilter)
print('Urban Misclassified Points', urbanMisclassified)
// We can also apply a filter to select all misclassified points
// Since we are comparing 2 properties agaist each-other,
// we need to use a binary filter
var misClassified = test.filter(ee.Filter.notEquals({
leftField:'classification', rightField:'landcover'}))
print('All Misclassified Points', misClassified)
// Display the misclassified points by styling them
var landcover = ee.List([0, 1, 2, 3])
var palette = ee.List(['gray','brown','blue','green'])
var misclassStyled = ee.FeatureCollection(
landcover.map(function(lc){
var feature = misClassified.filter(ee.Filter.eq('landcover', lc))
var color = palette.get(landcover.indexOf(lc));
var markerStyle = {color:color}
return feature.map(function(point){
return point.set('style', markerStyle)
})
})).flatten();
Map.addLayer(misclassStyled.style({styleProperty:"style"}), {}, 'Misclassified Points')Image Normalization and Standardization
For machine learning, it is a recommended practice to either normalize or standardize your features. The code below shows how to implement these feature scaling techniques.
var image = ee.Image("users/ujavalgandhi/e2e/arkavathy_2019_composite");
var boundary = ee.FeatureCollection("users/ujavalgandhi/e2e/arkavathy_boundary")
var geometry = boundary.geometry()
//**************************************************************************
// Function to Normalize Image
// Pixel Values should be between 0 and 1
// Formula is (x - xmin) / (xmax - xmin)
//**************************************************************************
function normalize(image){
var bandNames = image.bandNames();
// Compute min and max of the image
var minDict = image.reduceRegion({
reducer: ee.Reducer.min(),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var maxDict = image.reduceRegion({
reducer: ee.Reducer.max(),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var mins = ee.Image.constant(minDict.values(bandNames));
var maxs = ee.Image.constant(maxDict.values(bandNames));
var normalized = image.subtract(mins).divide(maxs.subtract(mins))
return normalized
}
//**************************************************************************
// Function to Standardize Image
// (Mean Centered Imagery with Unit Standard Deviation)
// https://365datascience.com/tutorials/statistics-tutorials/standardization/
//**************************************************************************
function standardize(image){
var bandNames = image.bandNames();
// Mean center the data to enable a faster covariance reducer
// and an SD stretch of the principal components.
var meanDict = image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means)
var stdDevDict = image.reduceRegion({
reducer: ee.Reducer.stdDev(),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var stddevs = ee.Image.constant(stdDevDict.values(bandNames));
var standardized = centered.divide(stddevs);
return standardized
}
var standardizedImage = standardize(image)
var normalizedImage = normalize(image)
Map.addLayer(image,
{bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3, gamma: 1.2}, 'Original Image');
Map.addLayer(normalizedImage,
{bands: ['B4', 'B3', 'B2'], min: 0, max: 1, gamma: 1.2}, 'Normalized Image');
Map.addLayer(standardizedImage,
{bands: ['B4', 'B3', 'B2'], min: -1, max: 2, gamma: 1.2}, 'Standarized Image');
Map.centerObject(geometry)
// Verify Normalization
var beforeDict = image.reduceRegion({
reducer: ee.Reducer.minMax(),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var afterDict = normalizedImage.reduceRegion({
reducer: ee.Reducer.minMax(),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
print('Original Image Min/Max', beforeDict)
print('Normalized Image Min/Max', afterDict)
// Verify Standadization
// Verify that the means are 0 and standard deviations are 1
var beforeDict = image.reduceRegion({
reducer: ee.Reducer.mean().combine({
reducer2: ee.Reducer.stdDev(), sharedInputs: true}),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
var resultDict = standardizedImage.reduceRegion({
reducer: ee.Reducer.mean().combine({
reducer2: ee.Reducer.stdDev(), sharedInputs: true}),
geometry: geometry,
scale: 20,
maxPixels: 1e9,
bestEffort: true,
tileScale: 16
});
// Means are very small franctions close to 0
// Round them off to 2 decimals
var afterDict = resultDict.map(function(key, value) {
return ee.Number(value).format('%.2f')
})
print('Original Image Mean/StdDev', beforeDict)
print('Standadized Image Mean/StdDev', afterDict)Time-Series Smoothing and Gap-filling
Moving Window Smoothing
A technique applied to a time series for removal of the fine-grained variation between time steps is known as Smoothing. This example shows how a moving-window smoothing algorithm can be applied in Earth Engine. Using a Save-all Join, the collection is joined with itself and all images that fall within the temporal-window are added as a property of each image. Next, a mean reducer is applied on all the images, resulting in the average value of the pixel within the time-frame. The resulting time-series reduces the sharp peaks and valleys - and is more robust against outliers (such as cloudy pixels)
Moving Window Average Smoothing
// Moving-Window Temporal Smoothing
var s2 = ee.ImageCollection("COPERNICUS/S2");
var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var rgbVis = {min: 0.0, max: 3000, bands: ['B4', 'B3', 'B2']};
var filtered = s2
.filter(ee.Filter.date('2017-01-01', '2018-01-01'))
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask).divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = filtered.map(maskS2clouds)
// Write a function that computes NDVI for an image and adds it as a band
function addNDVI(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
return image.addBands(ndvi);
}
// Map the function over the collection
var withNdvi = filtered.map(addNDVI);
// Moving-Window Smoothing
// Specify the time-window
var days = 30
// We use a 'join' to find all images that are within
// the time-window
// The join will add all matching images into a
// new property called 'images'
var join = ee.Join.saveAll({
matchesKey: 'images'
});
// This filter will match all images that are captured
// within the specified day of the source image
var diffFilter = ee.Filter.maxDifference({
difference: 1000 * 60 * 60 * 24 * days,
leftField: 'system:time_start',
rightField: 'system:time_start'
});
// Select the 'ndvi' band
var ndviCol = withNdvi.select('ndvi')
var joinedCollection = join.apply({
primary: ndviCol,
secondary: ndviCol,
condition: diffFilter
});
// Each image in the joined collection will contain
// matching images in the 'images' property
// Extract and return the mean of matched images
var smoothedCollection = ee.ImageCollection(joinedCollection.map(function(image) {
var collection = ee.ImageCollection.fromImages(image.get('images'));
return ee.Image(image).addBands(collection.mean().rename('moving_average'));
}))
// Display a time-series chart
var chart = ui.Chart.image.series({
imageCollection: smoothedCollection.select(['ndvi', 'moving_average']),
region: geometry,
reducer: ee.Reducer.mean(),
scale: 20
}).setOptions({
lineWidth: 1,
title: 'NDVI Time Series',
interpolateNulls: true,
vAxis: {title: 'NDVI'},
hAxis: {title: '', format: 'YYYY-MMM'},
series: {
1: {color: 'gray', lineDashStyle: [1, 1]}, // Original NDVI
0: {color: 'red', lineWidth: 2 }, // Smoothed NDVI
},
})
print(chart);Temporal Interpolation
The code below shows how to do temporal gap-filling of time-series data. A detailed explanation of the code and other examples are described in our blog post Temporal Gap-Filling with Linear Interpolation in GEE.
// Temporal Interpolation (Gap-Filling Masked Pixels)
var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var filtered = s2
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
function maskCloudAndShadowsSR(image) {
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(5);
var snow = snowProb.lt(5);
var scl = image.select('SCL');
var shadow = scl.eq(3); // 3 = cloud shadow
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 5% or cloud shadow classification
var mask = (cloud.and(snow)).and(cirrus.neq(1)).and(shadow.neq(1));
return image.updateMask(mask).divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"]);
}
var filtered = filtered.map(maskCloudAndShadowsSR)
// Add a band containing timestamp to each image
// This will be used to do pixel-wise interpolation later
var filtered = filtered.map(function(image) {
var timeImage = image.metadata('system:time_start').rename('timestamp')
// The time image doesn't have a mask.
// We set the mask of the time band to be the same as the first band of the image
var timeImageMasked = timeImage.updateMask(image.mask().select(0))
return image.addBands(timeImageMasked)
})
// Specify the time-window
// This will determine how much backward and forward are we willing to
// look for an unmasked pixel in the time-series
var days = 30
// For each image in the collection, we need to find all images
// before and after the specified time-window
// This is accomplished using Joins
// We need to do 2 joins
// Join 1: Join the collection with itself to find all images before each image
// Join 2: Join the collection with itself to find all images after each image
// We first define the filters needed for the join
// Define a maxDifference filter to find all images within the specified days
// The filter needs the time difference in milliseconds
// Convert days to milliseconds
var millis = ee.Number(days).multiply(1000*60*60*24)
var maxDiffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
})
// We need a lessThanOrEquals filter to find all images after a given image
// This will compare the given image's timestamp against other images' timestamps
var lessEqFilter = ee.Filter.lessThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
// We need a greaterThanOrEquals filter to find all images before a given image
// This will compare the given image's timestamp against other images' timestamps
var greaterEqFilter = ee.Filter.greaterThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
// Apply the joins
// For the first join, we need to match all images that are after the given image.
// To do this we need to match 2 conditions
// 1. The resulting images must be within the specified time-window of target image
// 2. The target image's timestamp must be lesser than the timestamp of resulting images
// Combine two filters to match both these conditions
var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)
// This join will find all images after, sorted in descending order
// This will gives us images so that closest is last
var join1 = ee.Join.saveAll({
matchesKey: 'after',
ordering: 'system:time_start',
ascending: false})
var join1Result = join1.apply({
primary: filtered,
secondary: filtered,
condition: filter1
})
// Each image now as a property called 'after' containing
// all images that come after it within the time-window
print(join1Result.first())
// Do the second join now to match all images within the time-window
// that come before each image
var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)
// This join will find all images before, sorted in ascending order
// This will gives us images so that closest is last
var join2 = ee.Join.saveAll({
matchesKey: 'before',
ordering: 'system:time_start',
ascending: true})
var join2Result = join2.apply({
primary: join1Result,
secondary: join1Result,
condition: filter2
})
// Each image now as a property called 'before' containing
// all images that come after it within the time-window
print(join2Result.first())
// Do the interpolation
// We now write a function that will be used to interpolate all images
// This function takes an image and replaces the masked pixels
// with the interpolated value from before and after images.
var interpolateImages = function(image) {
var image = ee.Image(image)
// We get the list of before and after images from the image property
// Mosaic the images so we a before and after image with the closest unmasked pixel
var beforeImages = ee.List(image.get('before'))
var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
var afterImages = ee.List(image.get('after'))
var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
// Interpolation formula
// y = y1 + (y2-y1)*((t – t1) / (t2 – t1))
// y = interpolated image
// y1 = before image
// y2 = after image
// t = interpolation timestamp
// t1 = before image timestamp
// t2 = after image timestamp
// We first compute the ratio (t – t1) / (t2 – t1)
// Get image with before and after times
var t1 = beforeMosaic.select('timestamp').rename('t1')
var t2 = afterMosaic.select('timestamp').rename('t2')
var t = image.metadata('system:time_start').rename('t')
var timeImage = ee.Image.cat([t1, t2, t])
var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
't': timeImage.select('t'),
't1': timeImage.select('t1'),
't2': timeImage.select('t2'),
})
// You can replace timeRatio with a constant value 0.5
// if you wanted a simple average
// Compute an image with the interpolated image y
var interpolated = beforeMosaic
.add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
// Replace the masked pixels in the current image with the average value
var result = image.unmask(interpolated)
return result.copyProperties(image, ['system:time_start'])
}
// map() the function to interpolate all images in the collection
var interpolatedCol = ee.ImageCollection(join2Result.map(interpolateImages))
// interpolatedCol is the gap-filled collection that can be used
// for further analysis.
// Let's visualize the NDVI time-series
function addNDVI(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
return image.addBands(ndvi);
}
var ndviVis = {min:0, max: 0.9, palette: [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
'66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
'012E01', '011D01', '011301'
]}
var visualizeImage = function(image) {
return image.select('ndvi').visualize(ndviVis).clip(geometry)
}
var visCollectionOriginal = filtered
.map(addNDVI)
.map(visualizeImage)
var visCollectionInterpolated = interpolatedCol
.map(addNDVI)
.map(visualizeImage)
var videoArgs = {
dimensions: 400,
region: geometry,
framesPerSecond: 1,
crs: 'EPSG:3857',
};
print('Original NDVI Images', visCollectionOriginal.getVideoThumbURL(videoArgs))
print('Gap Filled NDVI Images', visCollectionInterpolated.getVideoThumbURL(videoArgs))Savitzky-Golay Smoothing
The Savitzky–Golay filter fits a polynomial to a set of data points in a time-series. The Open Earth Engine Library (OEEL) provides an efficient implementation of this filter that can be applied on an ImageCollection. However, the time-series must be pre-processed so there are images at a regular interval. We use the interpolation technique described in the previous section and prepare a continous time-series without any masked pixels. The result is a new ImageCollection containing images at a regular interval (5-day) and with pixel values smoothed using the Savitzky–Golay filter.
Savitzky-Golay Smoothing
// Aplying Savitzky-Golay Filter on a NDVI Time-Series
// This script uses the OEEL library to apply a
// Savitzky-Golay filter on a imagecollection
// We require a regularly-spaced time-series without
// any masked pixels. So this script applies
// linear interpolation to created regularly spaced images
// from the original time-series
// Step-1: Prepare a NDVI Time-Series
// Step-2: Create an empty Time-Series with images at n days
// Step-3: Use Joins to find before/after images
// Step-4: Apply linear interpolation to fill each image
// Step-5: Apply Savitzky-Golay filter
// Step-6: Visualize the results
//##############################################################
// Step-1: Prepare a NDVI time-series
//##############################################################
var s2 = ee.ImageCollection("COPERNICUS/S2");
var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
Map.addLayer(geometry, {color: 'red'}, 'Farm')
Map.centerObject(geometry)
var startDate = ee.Date('2017-01-01')
var endDate = startDate.advance(1, 'year')
var filtered = s2
.filter(ee.Filter.date(startDate, endDate))
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
function maskS2clouds(image) {
var qa = image.select('QA60')
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask)//.divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
var filtered = filtered.map(maskS2clouds)
// Write a function that computes NDVI for an image and adds it as a band
function addNDVI(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).toFloat().rename('ndvi');
return image.addBands(ndvi);
}
// Map the function over the collection
var withNdvi = filtered.map(addNDVI);
// Select 'ndvi' band
var ndviCol = withNdvi.select('ndvi')
print('Original Collection', ndviCol)
//##############################################################
// Step-2: Create an empty Time-Series with images at n days
//##############################################################
// Select the interval. We will have 1 image every n days
var n = 5;
var totalDays = endDate.difference(startDate, 'day');
var daysToInterpolate = ee.List.sequence(1, totalDays, n)
var initImages = daysToInterpolate.map(function(day) {
var image = ee.Image().rename('ndvi').toFloat().set({
'system:index': ee.Number(day).format('%d'),
'system:time_start': startDate.advance(day, 'day').millis(),
// Set a property so we can identify interpolated images
'type': 'interpolated'
})
return image
})
var initCol = ee.ImageCollection.fromImages(initImages)
print('Empty Collection', initCol)
//##############################################################
// Step-3: Use Joins to find before/after images
//##############################################################
// Merge empty collection with the original collection so we can
// find images to interpolate from
var mergedCol = ndviCol.merge(initCol)
var mergedCol = mergedCol.map(function(image) {
var timeImage = image.metadata('system:time_start').rename('timestamp')
var timeImageMasked = timeImage.updateMask(image.mask().select(0))
return image.addBands(timeImageMasked)
})
// Specify the time-window
// Set it so that we have at least 1 non-cloudy image in the period
var days = 60
var millis = ee.Number(days).multiply(1000*60*60*24)
var maxDiffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var lessEqFilter = ee.Filter.lessThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var greaterEqFilter = ee.Filter.greaterThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)
var join1 = ee.Join.saveAll({
matchesKey: 'after',
ordering: 'system:time_start',
ascending: false})
var join1Result = join1.apply({
primary: mergedCol,
secondary: mergedCol,
condition: filter1
})
var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)
var join2 = ee.Join.saveAll({
matchesKey: 'before',
ordering: 'system:time_start',
ascending: true})
var join2Result = join2.apply({
primary: join1Result,
secondary: join1Result,
condition: filter2
})
//##############################################################
// Step-4: Apply linear interpolation to fill each image
//##############################################################
// Once the joins are done, we don't need original NDVI images
// We keep only the blank images which now have matching NDVI images
// as properties
var filtered = join2Result.filter(ee.Filter.eq('type', 'interpolated'))
// Interpolatinon function
function interpolateImages(image) {
var image = ee.Image(image)
var beforeImages = ee.List(image.get('before'))
var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
var afterImages = ee.List(image.get('after'))
var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
var t1 = beforeMosaic.select('timestamp').rename('t1')
var t2 = afterMosaic.select('timestamp').rename('t2')
var t = image.metadata('system:time_start').rename('t')
var timeImage = ee.Image.cat([t1, t2, t])
var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
't': timeImage.select('t'),
't1': timeImage.select('t1'),
't2': timeImage.select('t2'),
})
var interpolated = beforeMosaic
.add((afterMosaic.subtract(beforeMosaic).multiply(timeRatio)))
var result = image.unmask(interpolated)
return result.copyProperties(image, ['system:time_start'])
}
var interpolatedCol = ee.ImageCollection(
filtered.map(interpolateImages)).select('ndvi')
print('Interpolated Collection', interpolatedCol)
//##############################################################
// Step-5: Apply Savitzky-Golay filter
//##############################################################
var oeel=require('users/OEEL/lib:loadAll');
// https://www.open-geocomputing.org/OpenEarthEngineLibrary/#.ImageCollection.SavatskyGolayFilter
// Use the same maxDiffFilter we used earlier
var maxDiffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
})
// Use the default distanceFunction
var distanceFunction = function(infromedImage, estimationImage) {
return ee.Image.constant(
ee.Number(infromedImage.get('system:time_start'))
.subtract(
ee.Number(estimationImage.get('system:time_start')))
);
}
// Apply smoothing of order=3
var order = 3;
var smoothed = oeel.ImageCollection.SavatskyGolayFilter(
interpolatedCol,
maxDiffFilter,
distanceFunction,
order)
// Select the d_0_ndvi band and rename it
var smoothed = smoothed.select(['d_0_ndvi'], ['smoothed'])
//##############################################################
// Step-6: Visualize the results
//##############################################################
// Chart the time-series at a single location
var title = 'Savitsky-Golay smoothing' +
'(order = '+ order + ', window_size = ' + days + ')'
// Plot the original and fitted NDVI time-series
var chart = ui.Chart.image.series({
imageCollection: ndviCol.merge(smoothed),
region: geometry,
reducer: ee.Reducer.mean(),
scale: 20
}).setOptions({
lineWidth: 1,
title: title,
interpolateNulls: true,
vAxis: {title: 'NDVI'},
hAxis: {title: '', format: 'YYYY-MMM'},
series: {
0: {color: 'blue', lineWidth: 1,
lineDashStyle: [1, 1], pointSize: 1,
}, // Original NDVI
1: {color: 'red', lineWidth: 2 }, // Smoothed NDVI
},
})
print(chart)
// Create an animation to visualize the smoothed collection
var ndviVis = {min:0, max: 0.7, palette: [
'#fdae61','#fee08b','#d9ef8b',
'#a6d96a','#66bd63','#1a9850','#006837']
}
var visualizeImage = function(image) {
return image.visualize(ndviVis).clip(geometry)
}
var visCollectionOriginal = ndviCol
.map(visualizeImage)
var visCollectionInterpolated = smoothed
.map(visualizeImage)
var videoArgs = {
dimensions: 400,
region: geometry,
framesPerSecond: 2,
crs: 'EPSG:3857',
};
print('Original NDVI Images', visCollectionOriginal.getVideoThumbURL(videoArgs))
var videoArgs = {
dimensions: 400,
region: geometry,
framesPerSecond: 5,
crs: 'EPSG:3857',
};
print('SG Filtered NDVI Images', visCollectionInterpolated.getVideoThumbURL(videoArgs))User Interface Templates
Adding a Discrete Legend
You may want to add a legend for a classified image to your map visualization in your App. Here’s a code snippet that shows how you can build it using the UI Widgets.
Creating a Discrete Map Legend
var classified = ee.Image("users/ujavalgandhi/e2e/bangalore_classified")
Map.centerObject(classified)
Map.addLayer(classified,
{min: 0, max: 3, palette: ['gray', 'brown', 'blue', 'green']}, '2019');
var legend = ui.Panel({style: {position: 'middle-right', padding: '8px 15px'}});
var makeRow = function(color, name) {
var colorBox = ui.Label({
style: {color: '#ffffff',
backgroundColor: color,
padding: '10px',
margin: '0 0 4px 0',
}
});
var description = ui.Label({
value: name,
style: {
margin: '0px 0 4px 6px',
}
});
return ui.Panel({
widgets: [colorBox, description],
layout: ui.Panel.Layout.Flow('horizontal')}
)};
var title = ui.Label({
value: 'Legend',
style: {fontWeight: 'bold',
fontSize: '16px',
margin: '0px 0 4px 0px'}});
legend.add(title);
legend.add(makeRow('gray','Built-up'))
legend.add(makeRow('brown','Bare Earth'))
legend.add(makeRow('blue','Water'))
legend.add(makeRow('green','Vegetation'))
Map.add(legend);Adding a Continous Legend
If you are displaying a raster layer in your app with a color palette, you can use the following technique to add a colorbar using the snipet below.
Creating a Continuous Raster Legend
var s2 = ee.ImageCollection("COPERNICUS/S2");
var admin2 = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level2");
var bangalore = admin2.filter(ee.Filter.eq('ADM2_NAME', 'Bangalore Urban'))
var geometry = bangalore.geometry()
var filtered = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.date('2019-01-01', '2019-12-31'))
.filter(ee.Filter.bounds(geometry))
var image = filtered.median();
// Calculate Normalized Difference Vegetation Index (NDVI)
// 'NIR' (B8) and 'RED' (B4)
var ndvi = image.normalizedDifference(['B8', 'B4']).rename(['ndvi']);
var palette = ['#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641']
var ndviVis = {min:0, max:0.5, palette: palette}
Map.centerObject(geometry, 12)
Map.addLayer(ndvi.clip(geometry), ndviVis, 'ndvi')
function createColorBar(titleText, palette, min, max) {
// Legend Title
var title = ui.Label({
value: titleText,
style: {fontWeight: 'bold', textAlign: 'center', stretch: 'horizontal'}});
// Colorbar
var legend = ui.Thumbnail({
image: ee.Image.pixelLonLat().select(0),
params: {
bbox: [0, 0, 1, 0.1],
dimensions: '200x20',
format: 'png',
min: 0, max: 1,
palette: palette},
style: {stretch: 'horizontal', margin: '8px 8px', maxHeight: '40px'},
});
// Legend Labels
var labels = ui.Panel({
widgets: [
ui.Label(min, {margin: '4px 10px',textAlign: 'left', stretch: 'horizontal'}),
ui.Label((min+max)/2, {margin: '4px 20px', textAlign: 'center', stretch: 'horizontal'}),
ui.Label(max, {margin: '4px 10px',textAlign: 'right', stretch: 'horizontal'})],
layout: ui.Panel.Layout.flow('horizontal')});
// Create a panel with all 3 widgets
var legendPanel = ui.Panel({
widgets: [title, legend, labels],
style: {position: 'bottom-center', padding: '8px 15px'}
})
return legendPanel
}
// Call the function to create a colorbar legend
var colorBar = createColorBar('NDVI Values', palette, 0, 0.5)
Map.add(colorBar)Split Panel App Template
A common use-case for Earth Engine Apps is to display 2 images in a split panel app. Below script contains a simple template that you can use to create an interactive split panel app. Here we have 2 map objects - leftMap and rightMap. You can add different images to each of the maps and users will be able to explore them side-by-side. [View Animated GIF ↗]
A Split Panel App that displays Pre- and Post-Storm Images
// On June 9, 2018 - A massive dust storm hit North India
// This example shows before and after imagery from Sentinel-2
// Display two visualizations of a map.
// Set a center and zoom level.
var center = {lon: 77.47, lat: 28.41, zoom: 12};
// Create two maps.
var leftMap = ui.Map(center);
var rightMap = ui.Map(center);
// Remove UI controls from both maps, but leave zoom control on the left map.
leftMap.setControlVisibility(false);
rightMap.setControlVisibility(false);
leftMap.setControlVisibility({zoomControl: true});
// Link them together.
var linker = new ui.Map.Linker([leftMap, rightMap]);
// Create a split panel with the two maps.
var splitPanel = ui.SplitPanel({
firstPanel: leftMap,
secondPanel: rightMap,
orientation: 'horizontal',
wipe: true
});
// Remove the default map from the root panel.
ui.root.clear();
// Add our split panel to the root panel.
ui.root.add(splitPanel);
var rgb_vis = {min: 0, max: 3200, bands: ['B4', 'B3', 'B2']};
var preStorm = ee.Image('COPERNICUS/S2/20180604T052651_20180604T053435_T43RGM')
var postStorm = ee.Image('COPERNICUS/S2/20180614T052651_20180614T053611_T43RGM')
// Add a RGB Landsat 8 layer to the left map.
leftMap.addLayer(preStorm, rgb_vis);
rightMap.addLayer(postStorm, rgb_vis);Code Sharing and Script Modules
As your Earth Engine project grows, you need a way to organize and share your code to collaborate with others. We will learn some best practices on how best to set-up your project in Earth Engine.
Sharing a Single Script
To share your code from a single script, you need to use the Get Link button in the code editor. As you click the button, the contents of your code editor is captured and encoded into a URL. When you share this URL with someone, they will be able see same content as your code editor. This is a great way to send a snapshot of your code so others can reproduce your output. Remember that the script links are just snapshots, if you change your code after sending the link to someone, they will not see the updates.
When trying to send someone a link, do NOT click the Copy Script Path button. Sending this path to someone will NOT give them access to your code. The script path only works only on public or shared repositories.
Code Sharing using Get Link button
While sharing the script using Get Link, you should also share any private Assets that you may have uploaded and are using in the script. You can share the asset with a specific email address, or check the Anyone can read box if you want anyone with the script link to be able to access it. Failing to do so will prevent others from running your script.
Sharing Uploaded Assets
Learn more in the Script links section of the Google Earth Engine User Guide.
Sharing Multiple Scripts
If you want to share a collection of scripts with other users or your collaborators, the best way is to create a new Repository.
Creating New Repository
You can put multiple scripts within the repository and share the repository with other users. You can grant them Reader or Writer access so they can view/add/modify/delete scripts in that repository. If you want to make it readable by Public, you can check the Anyone can read option. You will see a URL in the form of https://code.earthengine.google.co.in/?accept_repo=.... When you share this URL with other users and they visit that link, your repository will be added to their Code Editor under the Reader or Writer folder depending on their access.
Creating New Repository
Learn more in the Script Manager section of the Google Earth Engine User Guide.
Sharing Code between Scripts
For a large project, it is preferable to share commonly used functions between scripts. That way, each script doesn’t have to re-implement the same code. Earth Engine enables this using Script Modules. Using a special object called exports, you can expose a function to other scripts. Learn more in the Script modules section of the Google Earth Engine User Guide.
There are many Earth Engine users who have shared their repositories publicly and written script modules to perform a variety of tasks. Here’s an example of using the grid module from the users/gena/packages repository to create regularly spaced grids in Earth Engine.
Using a function from a script module
var karnataka = ee.FeatureCollection("users/ujavalgandhi/public/karnataka");
Map.addLayer(karnataka, {color: 'gray'}, 'State Boundary')
var bounds = karnataka.geometry().bounds()
var coords = ee.List(bounds.coordinates().get(0))
var xmin = ee.List(coords.get(0)).get(0)
var ymin = ee.List(coords.get(0)).get(1)
var xmax = ee.List(coords.get(2)).get(0)
var ymax = ee.List(coords.get(2)).get(1)
// source code for the grid package:
// https://code.earthengine.google.com/?accept_repo=users/gena/packages
// Import the module to our script using 'require'
var gridModule = require('users/gena/packages:grid')
// Now we can run any function from the module
// We try running the generateGrid function to create regularly spaced vector grid
// generateGrid(xmin, ymin, xmax, ymax, dx, dy, marginx, marginy, opt_proj)
var spacing = 0.5
var gridVector = gridModule.generateGrid(xmin, ymin, xmax, ymax, spacing, spacing, 0, 0)
Map.centerObject(gridVector)
Map.addLayer(gridVector, {color: 'blue'}, 'Grids')Useful Public Repositories
Please visit Awesome Earth Engine to see a curated list of Google Earth Engine resources.
We also have a few recommendations of a few selected packages, which have very useful functions to make you productive in Earth Engine.
General Purpose Packages
- eepackages: A set of Google Earth Engine utilities by maintained by Gennadii Donchyts for both Javascript and Python API.
- geetools: Tools for cloud masking, batch processing, and more
- ee-palettes: Module for generating color palettes
- spectral: A javascript module that provides ready-to-use curated list of spectral indices for GEE.
- eemont: Python package that provides utility methods to create a more fluid code by being friendly with the Python method chaining.
Application Specific Packages
- LEAF-Toolbox: Google Earth Engine application that produces various Vegetation Biophysical Products, including Leaf Area Index (LAI).
- RivWidthCloud: Package to automate extracting river centerline and width for both Javascript and Python API.
Guided Projects
Below are step-by-step video-based walkthrough of implementing real-world projects using Earth Engine. You can continue their learning journey by implementing these projects for their region of interest after the class.
Get the Code
- Click this link to open Google Earth Engine code editor and add the repository to your account.
- If successful, you will have a new repository named
users/ujavalgandhi/End-to-End-Projectsin the Scripts tab in the Reader section.
If you do not see the repository in the Reader section, refresh your browser tab and it will show up.
Code Editor After Adding the Projects Repository
Project 1: Drought Monitoring
Calculating Rainfall Deviation from the 30-year mean using CHIRPS Gridded Rainfall Data
Project 3: Extracting Time-Series
Extracting a 10-year NDVI time-series over multiple polygons using MODIS data.
Learning Resources
- Awesome Earth Engine: A curated list of Google Earth Engine resources.
Data Credits
- Sentinel-2 Level-1C, Level-2A and Sentinel-1 SAR GRD: Contains Copernicus Sentinel data.
- TerraClimate: Monthly Climate and Climatic Water Balance for Global Terrestrial Surfaces, University of Idaho: Abatzoglou, J.T., S.Z. Dobrowski, S.A. Parks, K.C. Hegewisch, 2018, Terraclimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015, Scientific Data 5:170191, doi:10.1038/sdata.2017.191
- VIIRS Stray Light Corrected Nighttime Day/Night Band Composites Version 1: C. D. Elvidge, K. E. Baugh, M. Zhizhin, and F.-C. Hsu, “Why VIIRS data are superior to DMSP for mapping nighttime lights,” Asia-Pacific Advanced Network 35, vol. 35, p. 62, 2013.
- FAO GAUL 500m: Global Administrative Unit Layers 2015, Second-Level Administrative Units: Source of Administrative boundaries: The Global Administrative Unit Layers (GAUL) dataset, implemented by FAO within the CountrySTAT and Agricultural Market Information System (AMIS) projects.
- CHIRPS Pentad: Climate Hazards Group InfraRed Precipitation with Station Data (version 2.0 final): Funk, Chris, Pete Peterson, Martin Landsfeld, Diego Pedreros, James Verdin, Shraddhanand Shukla, Gregory Husak, James Rowland, Laura Harrison, Andrew Hoell & Joel Michaelsen. “The climate hazards infrared precipitation with stations—a new environmental record for monitoring extremes”. Scientific Data 2, 150066. doi:10.1038/sdata.2015.66 2015.
- MOD13Q1.006 Terra Vegetation Indices 16-Day Global 250m: Didan, K. (2015). MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 [Data set]. NASA EOSDIS Land Processes DAAC. Accessed 2021-05-06 from https://doi.org/10.5067/MODIS/MOD13Q1.006
- Natural Earth Urban Areas: Urban Areas Shapefile. Downloaded from Natural Earth. Free vector and raster map data @ naturalearthdata.com.
License
The course material (text, images, presentation, videos) is licensed under a Creative Commons Attribution 4.0 International License.
The code (scripts, Jupyter notebooks) is licensed under the MIT License. For a copy, see https://opensource.org/licenses/MIT
Kindly give appropriate credit to the original author as below:
Copyright © 2021 Ujaval Gandhi www.spatialthoughts.com
If you would like to white-label these materials as part of a commercial offering, you can purchase a Trainer License for a small fee.
Please contact us for pricing and terms.
Citing and Referencing
You can cite the course materials as follows
- Gandhi, Ujaval, 2021. End-to-End Google Earth Engine Course. Spatial Thoughts. https://courses.spatialthoughts.com/end-to-end-gee.html
This course is offered as an instructor-led online class. Visit Spatial Thoughts to know details of upcoming sessions.