Google Earth Engine (GEE) Production Tutorials
Planetary-scale remote sensing masterclass: Sentinel-2 Level-2A SCL cloud masking, Sealed Urban Index (SUI) time-series charting, Random Forest classification, and cloud batch exports.
Massively parallel distributed map-reduce computation over multi-petabyte public Earth Observation catalogs.
Bottom-of-Atmosphere (BOA) Surface Reflectance with Scene Classification Layer (SCL) quality flags.
Interactive web Code Editor scripts and headless Python automation via geemap and ee libraries.
Distributed non-parametric decision tree ensemble for multi-temporal land cover classification.
Tutorial 1: Sentinel-2 Level-2A SCL Cloud & Shadow Masking
Reliable surface monitoring requires removing cloud cover and cloud shadows from multi-temporal optical imagery. While raw Top-of-Atmosphere (TOA, Level-1C) scenes feature only a coarse 60-meter QA60 bitmask, the Copernicus Sentinel-2 Level-2A Bottom-of-Atmosphere (BOA) collection includes a 20-meter Scene Classification Layer (SCL) generated by the Sen2Cor processor.
SCL Quality Classification Values
The SCL band assigns each pixel an integer value from 0 to 11:
- Values to Mask Out:
0(No Data),1(Saturated/Defective),3(Cloud Shadows),8(Cloud Medium Probability),9(Cloud High Probability),10(Thin Cirrus),11(Snow/Ice). - Valid Land Cover Classes:
4(Vegetation),5(Bare Soil),6(Water Bodies),7(Unclassified / Built-Up Features).
Production GEE JavaScript Cloud Masking Code
/**
* Sentinel-2 Level-2A SCL Cloud Masking Pipeline
* Platform: Google Earth Engine JavaScript API
*/
// 1. Define Area of Interest (Mediterranean Coastal Basin)
var aoi = ee.Geometry.Point([16.86, 41.12]).buffer(25000);
// 2. Define SCL Cloud Masking Function
function maskS2sr(image) {
var scl = image.select("SCL");
// Create mask keeping only valid clear land pixels:
// 4 = Vegetation, 5 = Bare Soil, 6 = Water, 7 = Unclassified/Built
var clearMask = scl.eq(4)
.or(scl.eq(5))
.or(scl.eq(6))
.or(scl.eq(7));
// Surface reflectance values are scaled by 10,000 in Harmonized collection
return image.updateMask(clearMask)
.divide(10000)
.select(["B2", "B3", "B4", "B8", "B11", "B12"], ["blue", "green", "red", "nir", "swir1", "swir2"])
.copyProperties(image, ["system:time_start"]);
}
// 3. Load Harmonized Level-2A Collection
var s2Collection = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
.filterBounds(aoi)
.filterDate("2026-05-01", "2026-09-30")
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 30))
.map(maskS2sr);
print("Total Cloud-Masked Scenes Loaded:", s2Collection.size());
// 4. Generate Cloud-Free Summer Median Composite
var medianComposite = s2Collection.median().clip(aoi);
// 5. Visualize Natural Color (RGB) and False Color Infrared (CIR)
Map.centerObject(aoi, 11);
Map.addLayer(medianComposite, {bands: ["red", "green", "blue"], min: 0.0, max: 0.25}, "Summer True Color (RGB)");
Map.addLayer(medianComposite, {bands: ["nir", "red", "green"], min: 0.0, max: 0.40}, "Color Infrared (CIR)");
Tutorial 2: Sealed Urban Index (SUI) & Multi-Year Time-Series Charting
Accurately mapping artificial impervious soil sealing requires distinguishing stable built-up concrete from dry fallow agricultural soils that undergo seasonal plowing. The Sealed Urban Index (SUI) developed by Ahmad Asad et al. leverages four spectral bands (Green, Red, SWIR1, SWIR2) to suppress bare mineral soils while enhancing artificial structures.
Mathematical Formulation in Earth Engine
In GEE, implementing complex multi-band ratio expressions is optimized using ee.Image.expression():
// Compute Sealed Urban Index (SUI) with epsilon zero-division guard
var sui = image.expression(
"((GREEN - SWIR1) * (SWIR2 - RED)) / " +
"(max(GREEN + SWIR1, 1e-6) * max(SWIR2 + RED, 1e-6)) + 0.5",
{
"GREEN": image.select("green"),
"RED": image.select("red"),
"SWIR1": image.select("swir1"),
"SWIR2": image.select("swir2")
}
).rename("SUI");
Extracting Multi-Temporal Urban Phenology Charts
The following script maps the SUI calculation across a multi-year collection and plots monthly mean index values across an urban core vs an adjacent agricultural zone:
// Add SUI calculation to collection
var s2WithSui = s2Collection.map(function(img) {
var sui = img.expression(
"((g - s1) * (s2 - r)) / (max(g + s1, 1e-6) * max(s2 + r, 1e-6)) + 0.5",
{
"g": img.select("green"),
"r": img.select("red"),
"s1": img.select("swir1"),
"s2": img.select("swir2")
}
).rename("SUI");
return img.addBands(sui);
});
// Define inspection samples: Urban Industrial vs Cropland
var urbanSample = ee.Geometry.Point([16.84, 41.10]).buffer(100);
var cropSample = ee.Geometry.Point([16.92, 41.05]).buffer(100);
// Generate Time-Series Chart
var chart = ui.Chart.image.seriesByRegion({
imageCollection: s2WithSui.select("SUI"),
regions: ee.FeatureCollection([
ee.Feature(urbanSample, {label: "Industrial Sealed Roofs"}),
ee.Feature(cropSample, {label: "Rainfed Farmland"})
]),
reducer: ee.Reducer.mean(),
band: "SUI",
scale: 10,
xProperty: "system:time_start",
seriesProperty: "label"
}).setOptions({
title: "Multi-Temporal SUI Trajectory: Urban Sealing vs Agricultural Phenology",
hAxis: {title: "Date", format: "YYYY-MM"},
vAxis: {title: "Sealed Urban Index (SUI)", viewWindow: {min: 0.3, max: 0.7}},
colors: ["#ef4444", "#10b981"],
lineWidth: 2
});
print(chart);
Tutorial 3: Supervised Land Cover Classification with Random Forest
Supervised classification categorizes spectral pixels into discrete thematic land cover classes. In Google Earth Engine, the non-parametric Random Forest classifier (ee.Classifier.smileRandomForest) handles multi-collinear spectral bands and non-linear feature spaces with high generalization accuracy.
Step-by-Step Classification Pipeline
- Prepare Predictor Stack: Assemble Sentinel-2 spectral bands (Blue, Green, Red, NIR, SWIR1, SWIR2) and derived indices (NDVI, NDWI, SUI) into an 8-band composite image.
- Collect Ground Truth Polygons: Digitized reference geometries representing 4 classes: 0 = Water, 1 = Dense Canopy Forest, 2 = Agricultural Crops, 3 = Sealed Urban Infrastructure.
- Sample Training Pixels: Extract predictor pixel values across training polygons using
image.sampleRegions(). - Train Classifier: Instantiate
ee.Classifier.smileRandomForest(100)with 100 decision trees. - Classify & Assess Accuracy: Apply the trained model across the study area and compute the Confusion Matrix on independent validation polygons.
// 1. Assemble Predictor Bands Stack
var ndvi = medianComposite.normalizedDifference(["nir", "red"]).rename("NDVI");
var ndwi = medianComposite.normalizedDifference(["green", "nir"]).rename("NDWI");
var suiBand = s2WithSui.select("SUI").median();
var featureStack = medianComposite
.addBands(ndvi)
.addBands(ndwi)
.addBands(suiBand);
var bands = ["blue", "green", "red", "nir", "swir1", "swir2", "NDVI", "NDWI", "SUI"];
// 2. Training Sample Extraction (Assume trainingPoints feature collection exists)
var training = featureStack.select(bands).sampleRegions({
collection: trainingPolygons,
properties: ["landcover_class"],
scale: 10
});
// 3. Split 80% Training / 20% Validation
var withRandom = training.randomColumn("random");
var trainSet = withRandom.filter(ee.Filter.lt("random", 0.8));
var testSet = withRandom.filter(ee.Filter.gte("random", 0.8));
// 4. Train Random Forest Classifier
var classifier = ee.Classifier.smileRandomForest({
numberOfTrees: 100,
variablesPerSplit: 3,
minLeafPopulation: 2
}).train({
features: trainSet,
classProperty: "landcover_class",
inputProperties: bands
});
// 5. Classify the Regional Composite
var classified = featureStack.select(bands).classify(classifier);
// 6. Confusion Matrix and Accuracy Assessment
var validated = testSet.classify(classifier);
var testAccuracy = validated.errorMatrix("landcover_class", "classification");
print("Validation Error Matrix:", testAccuracy);
print("Overall Classification Accuracy:", testAccuracy.accuracy());
print("Kappa Coefficient:", testAccuracy.kappa());
// 7. Render Map
var palette = ["#0284c7", "#15803d", "#eab308", "#dc2626"];
Map.addLayer(classified, {min: 0, max: 3, palette: palette}, "Land Cover Classification");
Tutorial 4: Large-Scale Cloud Exports to Google Drive & GCS
Interactive browser maps in Earth Engine are capped at 32 megabytes of raw pixel data per request. To export high-resolution 10-meter production orthomosaics spanning whole provinces or river basins, jobs must be dispatched to backend batch export pipelines:
// Dispatch Asynchronous Batch Export Task to Google Drive
Export.image.toDrive({
image: classified.toInt8(),
description: "Regional_LandCover_Classification_2026",
folder: "GEE_Exports",
fileNamePrefix: "puglia_classified_10m",
region: aoi,
scale: 10,
crs: "EPSG:32633", // Target metric UTM Zone 33N
maxPixels: 1e13,
fileFormat: "GeoTIFF",
formatOptions: {
cloudOptimized: true // Generates Cloud-Optimized GeoTIFF with internal pyramids
}
});
Common Pitfalls and Field Troubleshooting in GEE
Developing in Earth Engine requires adhering to its distributed map-reduce model:
1. The Client-Side JavaScript Loop Anti-Pattern
Writing a standard JavaScript for (var i = 0; i < list.length; i++) to loop through images in an ee.ImageCollection causes client browser crashes. The client loop forces the browser to synchronously query the server on every iteration. Always use server-side mapping functions (collection.map(function(img) { ... })) so the operations execute concurrently across hundreds of Google Cloud worker nodes.
2. 'User Memory Limit Exceeded' in Reducers
Calling image.reduceRegion() on a vast continental polygon with thousands of vertices exceeds the 100 MB memory allocated to an interactive browser request. Solve this by:
- Increasing the
scaleparameter (e.g., from 10m to 30m or 60m for regional overviews). - Setting
tileScale: 4ortileScale: 8in the reducer parameters to divide the spatial geometry into smaller computational tiles. - Submitting the aggregation to
Export.table.toDrive()as a batch task rather than printing results interactively to the browser console.
3. False Cloud Detections Over High-Albedo Urban Roofs
In bright Mediterranean environments, large white industrial warehouse rooftops and solar panels exhibit high visible reflectance that automated cloud masks (such as SCL classes 8 and 9) occasionally misidentify as clouds. To prevent eroding urban structures, combine the SCL mask with thermal infrared criteria from Landsat 8-9 TIRS (clouds are cold, $< 280\text{ K}$, while summer tin roofs exceed $320\text{ K}$).
Frequently Asked Questions on Google Earth Engine Workflows
Why does using standard JavaScript for-loops on Earth Engine collections cause browser memory crashes?
Google Earth Engine operates on a client-server architecture. Variables declared with ee.Image or ee.ImageCollection represent server-side proxy objects executed on Google Cloud compute clusters. Standard client-side JavaScript for-loops run on the local browser thread, attempting to synchronously fetch server-side objects one by one. This causes severe network latency and browser out-of-memory errors. Developers must use server-side mapping functions like collection.map() to execute operations in parallel across Google Cloud workers.
What is the physical difference between Sentinel-2 SCL cloud masking and the older QA60 band?
The older QA60 band is a coarse 60-meter bitmask that flags opaque clouds and cirrus with frequent false negatives in littoral and mountain zones. The Scene Classification Layer (SCL) is a 20-meter product generated by the Sen2Cor atmospheric correction processor. SCL classifies each pixel into 12 discrete categories, including cloud shadows, medium/high probability clouds, thin cirrus, water, and snow, enabling much cleaner, higher-resolution cloud removal without losing cloud-free land pixels.
How does omitting the scale parameter in GEE reducers lead to inaccurate spatial statistics?
If the scale parameter is omitted in functions like image.reduceRegion(), Earth Engine defaults to the native pyramid resolution of the display canvas or a coarse 1000-meter grid. For 10-meter Sentinel-2 data, evaluating an urban polygon at 1000m scale aggregates thousands of diverse pixels into a single blended value, obliterating local land cover boundaries and producing wildly incorrect zonal averages.
How can 'User memory limit exceeded' errors be avoided when training Random Forest models on large regions?
Memory limit errors occur when the number of training sample points multiplied by the number of spectral bands exceeds the 100 MB memory allocated to an interactive worker thread. Mitigate this by: 1) Reducing training sample counts to 500 to 2,000 points per class; 2) Limiting predictor bands to essential features; 3) Using image.sample() with a defined scale; and 4) Exporting the classifier training job to Google Cloud Storage or Google Drive as an asynchronous batch task.
Why is ee.Image.expression() computationally superior to chained arithmetic methods in GEE?
Chaining individual arithmetic methods (such as image.select('B3').subtract(image.select('B11')).divide(...)) constructs deep, deeply nested Abstract Syntax Tree (AST) graphs on the server, increasing serialization overhead. The ee.Image.expression() method compiles mathematical band equations directly into an optimized C++ kernel on the backend, reducing memory footprints and speeding up computation by 30% to 50%.