Demo 2: Land Cover Classification - Learn Python with BNHR

This demonstration utilizes scripts originally developed for the 4-Day Advanced Online Training on AI/ML for Earth Observation for Philippine EO Professionals, part of the CoPhil Programme (EU-Philippines Copernicus Capacity Support Programme), a flagship initiative under the European Union’s Global Gateway strategy (2025).

Learning Objectives¶

Set up and authenticate Google Earth Engine in Google Colab
Acquire and preprocess Sentinel-2 satellite imagery
Calculate spectral indices (NDVI, NDWI, NDBI, EVI)
Create training datasets using verified open data sources
Train a Random Forest classifier
Perform land cover classification
Assess accuracy and generate statistics

Supervised Classification Workflow¶

Key Steps¶

Preprocessing: Cloud masking, atmospheric correction
Feature Extraction: Spectral bands, indices (NDVI, NDWI)
Training Data: Collect representative samples for each class
Model Training: Train classifier on training data
Classification: Apply model to entire scene
Validation: Assess accuracy with independent test data

Why Random Forest for Earth Observation?¶

Random Forest is one of the most popular algorithms for land cover classification because:

Handles high-dimensional data: Works well with many spectral bands (Sentinel-2 has 13 bands)
Robust to overfitting: Ensemble approach reduces variance
Feature importance: Reveals which bands are most informative
No feature scaling required: Unlike neural networks
Fast training: Efficient even with large datasets
Interpretable: Can visualize decision rules\

Random Forest: The “Forest” Analogy¶

One tree (Decision Tree): One expert’s opinion

Can be very confident but sometimes wrong
Might overfit to specific training examples

Forest (Random Forest): Committee of experts

Each expert sees slightly different data
Each expert considers different features
Final decision: Majority vote
Wisdom of crowds: Group decision more reliable than individual

💡 If you ask 100 independent experts and 75 say “Forest”, you can be more confident than if only 1 expert says “Forest”.

A. Setup and Authentication¶

This code block initializes the environment by linking Google Earth Engine for processing and visualizing satellite imagery.

# Install required packages
pip install earthengine-api geemap pandas numpy matplotlib seaborn -q

# Import all libraries

# Geospatial Analysis Libraries
import ee                           # Google Earth Engine API
import geemap                       # Interactive mapping and GEE utility tools

# Data Manipulation & Math
import pandas as pd                 # DataFrames for tabular data analysis
import numpy as np                  # Numerical operations and array handling

# Data Visualization
import matplotlib.pyplot as plt     # Standard plotting engine
import seaborn as sns               # Statistical data visualization

# Utility Libraries
from datetime import datetime       # Date and time handling for time-series

# Set plot style
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("husl")

print("✓ Libraries imported successfully")

This code prompt you through an authentication process and will open a browser window for account selection and authorization. Make sure that you select the email you provided when you registered to this training. If not, inform your instructor/s if you use another e-mail account.

# Authenticate and initialize Earth Engine
ee.Authenticate()
ee.Initialize(project='test-1-459001')

print("✓ Earth Engine initialized successfully")

You may encounter an error: <EEException: ee.Initialize: no project found>, in case you cannot authenticate and initialize Earth Engine.

In this case, you may need to provide a project that you own or have permissions to use. Go to your GEE account (you may also check in your Google Cloud) and choose any of the projects that you have and then copy the selected Project ID.

B. Define Study Area¶

This code block defines the geographic bounding box for Palawan, Philippines and initializes an interactive map centered on its coordinates, displaying the defined boundary as a red polygon layer for visual verification.

# Define Palawan boundary
palawan_coords = [
    [117.0, 8.5],   # Southwest
    [117.0, 12.5],  # Northwest
    [120.5, 12.5],  # Northeast
    [120.5, 8.5],   # Southeast
    [117.0, 8.5]    # Close polygon
]

palawan = ee.Geometry.Polygon(palawan_coords)

# Calculate area
area_km2 = palawan.area().divide(1e6).getInfo()
print(f"Study Area: Palawan Province")
print(f"Approximate Area: {area_km2:,.0f} km²")

# Create interactive map
Map = geemap.Map(center=[10.5, 118.8], zoom=8, height='600px')
Map.addLayer(palawan, {'color': 'red'}, 'Palawan Boundary')
Map

C. Sentinel-2 Data Acquisition¶

This code block filters the Sentinel-2 Harmonized collection for low-cloud imagery within the 2024 date range, applies cloud masking, and generates a median composite clipped to Palawan for a clear, seamless RGB visualization.

# Define date range (using 2024 for better availability)
start_date = '2024-01-01'
end_date = '2024-12-31'

# Cloud masking function
def mask_s2_clouds(image):
    qa = image.select('QA60')
    cloud_bit_mask = 1 << 10
    cirrus_bit_mask = 1 << 11
    mask = (qa.bitwiseAnd(cloud_bit_mask).eq(0)
            .And(qa.bitwiseAnd(cirrus_bit_mask).eq(0)))
    return image.updateMask(mask).divide(10000)

print(f"Date range: {start_date} to {end_date}")

# Load Sentinel-2 collection
s2_collection = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                 .filterBounds(palawan)
                 .filterDate(start_date, end_date)
                 .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                 .map(mask_s2_clouds))

# Create median composite
s2_composite = s2_collection.median().clip(palawan)

# Select bands
bands = ['B2', 'B3', 'B4', 'B8', 'B11', 'B12']
s2_composite = s2_composite.select(bands)

print(f"✓ Sentinel-2 composite created")
print(f"Images in collection: {s2_collection.size().getInfo()}")

# Visualize RGB composite
rgb_vis = {'bands': ['B4', 'B3', 'B2'], 'min': 0.0, 'max': 0.3, 'gamma': 1.4}
Map.addLayer(s2_composite, rgb_vis, 'Sentinel-2 RGB')
Map

D. Calculate Spectral Indices¶

This code block calculates four key spectral indices—

NDVI (vegetation),
NDWI (water),
NDBI (built-up areas), and
EVI (enhanced vegetation)—

by performing mathematical transformations on specific Sentinel-2 bands to highlight different environmental features.

# Calculate all indices
ndvi = s2_composite.normalizedDifference(['B8', 'B4']).rename('NDVI')
ndwi = s2_composite.normalizedDifference(['B3', 'B8']).rename('NDWI')
ndbi = s2_composite.normalizedDifference(['B11', 'B8']).rename('NDBI')

evi = s2_composite.expression(
    '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))',
    {
        'NIR': s2_composite.select('B8'),
        'RED': s2_composite.select('B4'),
        'BLUE': s2_composite.select('B2')
    }
).rename('EVI')

print("✓ Spectral indices calculated")

E. Feature Stack Preparation¶

This code block combines the original Sentinel-2 spectral bands with the four (4) newly calculated indices into a single multi-band image stack, creating a comprehensive feature set for spatial analysis and machine learning classification.

# Create feature stack
feature_stack = (s2_composite
                 .addBands(ndvi)
                 .addBands(ndwi)
                 .addBands(ndbi)
                 .addBands(evi))

feature_names = feature_stack.bandNames().getInfo()
print(f"✓ Feature stack created with {len(feature_names)} features")
print(f"Features: {feature_names}")

F. Training Data from Verified Open Sources¶

Using verified datasets:

ESA WorldCover 2021: 10m resolution, 76.7% accuracy
Global Mangrove Watch: 25m resolution, 87.4% accuracy for mangroves
OpenStreetMap: Urban features

# Define class properties
class_info = {
    1: {'name': 'Forest', 'color': '006400'},
    2: {'name': 'Agriculture', 'color': 'FFFF00'},
    3: {'name': 'Water', 'color': '0000FF'},
    4: {'name': 'Urban', 'color': 'FF0000'},
    5: {'name': 'Mangrove', 'color': '008B8B'}
}

print("Land Cover Classes:")
for class_id, info in class_info.items():
    print(f"  {class_id}: {info['name']}")

# Load ESA WorldCover 2021
worldcover = ee.ImageCollection('ESA/WorldCover/v200') \
    .filterDate('2021-01-01', '2022-01-01') \
    .first() \
    .clip(palawan)

# Remap to our classes
worldcover_remapped = worldcover.remap(
    [10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100],
    [1,  1,  1,  2,  4,  0,  0,  3,  0,  5,  0],
    0
).rename('landcover')

print("✓ ESA WorldCover loaded (10m, 76.7% accuracy)")

# Generate stratified random training points
training_points = worldcover_remapped.stratifiedSample(
    numPoints=500,
    classBand='landcover',
    region=palawan,
    scale=30,
    seed=42,
    geometries=True
).filter(ee.Filter.neq('landcover', 0))

print(f"✓ Training points generated: {training_points.size().getInfo()} total")

# Add urban training points from known cities
urban_centers = [
    ee.Geometry.Point([118.7384, 9.7392]).buffer(5000),   # Puerto Princesa
    ee.Geometry.Point([119.507, 10.693]).buffer(2000),     # Taytay
    ee.Geometry.Point([119.408, 11.041]).buffer(2000),     # El Nido
    ee.Geometry.Point([118.3707, 9.3085]).buffer(2000),    # Bataraza
]

urban_areas = ee.FeatureCollection([ee.Feature(center, {'landcover': 4}) for center in urban_centers])

urban_points = ee.FeatureCollection.randomPoints(
    region=urban_areas,
    points=200,
    seed=42
).map(lambda feat: feat.set('landcover', 4))

# Combine all training points
all_training_points = training_points.merge(urban_points)
print(f"✓ Final training dataset: {all_training_points.size().getInfo()} points")

# Sample spectral values from training points
training_samples = feature_stack.sampleRegions(
    collection=all_training_points,
    properties=['landcover'],
    scale=10,
    geometries=False,
    tileScale=4
)

print(f"✓ Training samples extracted: {training_samples.size().getInfo()} pixels")

# Visualize training points
Map_training = geemap.Map(center=[10.5, 118.8], zoom=8, height='600px')
Map_training.addLayer(s2_composite, rgb_vis, 'Sentinel-2 RGB')

for class_id, info in class_info.items():
    class_points = all_training_points.filter(ee.Filter.eq('landcover', class_id))
    Map_training.addLayer(class_points, {'color': info['color']}, f"Training: {info['name']}")

Map_training

G. Random Forest Training¶

This code block configures and trains a Random Forest classifier using 100 decision trees to learn the relationship between the input spectral features and the corresponding land cover training labels.

Parameters¶

n_trees

Description: Number of trees in forest
Typical values: 50-500
Impact: More trees → better performance (diminishing returns)

max_depth

Description: Maximum depth of each tree
Typical values: 10-50 or None
Impact: Deeper → more complex, risk overfitting

min_samples_split

Description: Min samples to split node
Typical values: 2-10
Impact: Higher → simpler trees, less overfitting

max_features

Description: Features per split
Typical values: 𝑛√ or log2(𝑛)
Impact: Balance between accuracy and diversity

bootstrap

Description: Use bootstrap sampling
Typical values: True
Impact: Almost always True for RF

💡 GEE’s ee.Classifier.smileRandomForest() defaults:

numberOfTrees: 100
variablesPerSplit: 𝑛√ (automatic)
minLeafPopulation: 1

# Train Random Forest classifier
rf_classifier = ee.Classifier.smileRandomForest(
    numberOfTrees=100,                            # The number of decision trees to create.
    variablesPerSplit=None,                       # The number of variables per split. If unspecified, uses the √ of the number of variables.
    minLeafPopulation=1,                          # Only create nodes whose training set contains at least this many points
    bagFraction=0.632,                            # The fraction of input to bag per tree.
    seed=42                                       # The randomization seed
)

trained_classifier = rf_classifier.train(
    features=training_samples,
    classProperty='landcover',
    inputProperties=feature_names
)

print("✓ Random Forest training complete")

H. Image Classification¶

This code block applies the trained Random Forest model to the entire multi-band feature stack to generate a discrete land cover map, which is then displayed using a custom color palette for visual analysis.

# Apply classifier
classified_image = feature_stack.classify(trained_classifier)

# Define palette
class_palette = [class_info[i]['color'] for i in sorted(class_info.keys())]

# Visualize
Map_classified = geemap.Map(center=[10.5, 118.8], zoom=8, height='700px')
Map_classified.addLayer(s2_composite, rgb_vis, 'Sentinel-2 RGB', False)
Map_classified.addLayer(classified_image, {'min': 1, 'max': 5, 'palette': class_palette}, 'Land Cover')

print("✓ Classification complete")
Map_classified

I. Accuracy Assessment¶

This code block evaluates the model’s performance by generating an independent validation dataset and comparing the predicted land cover classes against known reference values. It calculates a confusion matrix to derive the Overall Accuracy and Kappa Coefficient, providing a statistical measure of how reliably the classifier identifies different land cover types.

Overall Accuracy

Definition: Percentage of correctly classified samples
Formula: Overall Accuracy = Correct Predictions / Total Samples
Interpretation: Simple, intuitive metric
Limitation: Can be misleading with imbalanced classes
Example: 95% accuracy sounds great, but if 95% of pixels are forest, a “classify everything as forest” model achieves 95%!

Kappa Coefficient

Definition: Measure of agreement beyond chance
Interpretation: κ < 0.4: Poor agreement κ = 0.4-0.6: Moderate agreement κ = 0.6-0.8: Good agreement κ > 0.8: Excellent agreement
Why use Kappa? Accounts for agreement by random chance, especially important for imbalanced datasets.

# Create validation dataset
validation_points = worldcover_remapped.stratifiedSample(
    numPoints=100,
    classBand='landcover',
    region=palawan,
    scale=30,
    seed=123,  # Different seed
    geometries=True
).filter(ee.Filter.neq('landcover', 0))

validation_samples = feature_stack.sampleRegions(
    collection=validation_points,
    properties=['landcover'],
    scale=10,
    geometries=False
)

# Classify validation set
validation_classified = validation_samples.classify(trained_classifier)

# Generate confusion matrix
confusion_matrix = validation_classified.errorMatrix('landcover', 'classification')
overall_accuracy = confusion_matrix.accuracy().getInfo()
kappa = confusion_matrix.kappa().getInfo()

print(f"Overall Accuracy: {overall_accuracy*100:.2f}%")
print(f"Kappa Coefficient: {kappa:.4f}")

Summary and Key Takeaways¶

Random Forest Ensemble¶

Combines many trees to reduce variance and improve stability
Uses bootstrap sampling (bagging) for training diversity
Uses random feature selection to decorrelate trees
Final prediction by majority voting (classification) or averaging (regression)
More robust than single trees, less prone to overfitting

Accuracy Metrics¶

Overall accuracy can hide class-specific problems
Other accuracy metrics includes: (1) Precision (user’s accuracy): Reliability of positive predictions, (2) Recall (producer’s accuracy): Completeness of detection, and (3) F1-score: Harmonic mean balancing precision and recall
Choice of metric depends on application cost (false positives vs. false negatives)

For Earth Observation¶

Random Forest works well with multi-spectral data
Feature importance reveals spectral signature insights
Confusion patterns often reflect spectral similarity
Fast training enables rapid iteration

References¶

CoPhil EO AI/ML Training Programme. (2025). CoPhil EO AI/ML Training: 4-Day Advanced Training on AI/ML for Earth Observation. https://skotsopoulos-del.github.io/cophil-training-v1.0/

End of Demo 2 Notebook