Module 8: Land Use & Land Cover Classification

Mapping land use and land cover (LULC) from satellite imagery is one of the most impactful applications of Earth observation. This module covers supervised classification with Random Forest, accuracy assessment using confusion matrices and the Kappa statistic, and change detection for monitoring urban expansion, deforestation, and agricultural conversion.

1. Supervised Classification & Accuracy Metrics

Supervised classification assigns each pixel to a predefined land cover class based on its spectral, textural, and contextual features. The classifier learns the mapping from labeled training samples. The confusion matrix (error matrix) is the standard tool for accuracy assessment, from which we derive:

$$OA = \frac{\sum_i N_{ii}}{N_{total}}, \quad \kappa = \frac{OA - P_e}{1 - P_e}$$

where $N_{ii}$ are the diagonal elements of the confusion matrix (correctly classified pixels), $N_{total}$ is the total number of validation pixels, and $P_e$ is the expected agreement by chance:

$$P_e = \frac{1}{N_{total}^2} \sum_i N_{i+} \cdot N_{+i}$$

where $N_{i+}$ and $N_{+i}$ are the row and column marginal totals, respectively. Kappa values above 0.80 indicate strong agreement, 0.60–0.80 substantial agreement, and below 0.40 poor agreement. Per-class metrics include:

Producer's Accuracy (Recall)

The probability that a reference pixel of class $i$ is correctly classified: $PA_i = N_{ii} / N_{+i}$. Measures omission error ($1 - PA_i$).

User's Accuracy (Precision)

The probability that a pixel classified as class $i$ actually belongs to that class: $UA_i = N_{ii} / N_{i+}$. Measures commission error ($1 - UA_i$).

Sample Size Requirements

A common rule of thumb is 50–100 validation samples per class for reliable accuracy assessment. The samples should be independent of the training data (spatially separated) and collected using a probability-based design (stratified random sampling is recommended by Olofsson et al., 2014) to enable unbiased area estimation and confidence intervals.

2. Feature Stack for Random Forest Classification

A rich feature stack improves classification accuracy by providing the classifier with diverse spectral, index, and textural information. For Sentinel-2 based LULC classification, a typical 14-feature stack includes:

14-Feature Stack

Spectral Bands (10)

  • B2 (Blue, 490 nm) — 10 m
  • B3 (Green, 560 nm) — 10 m
  • B4 (Red, 665 nm) — 10 m
  • B5 (Red Edge 1, 705 nm) — 20 m
  • B6 (Red Edge 2, 740 nm) — 20 m
  • B7 (Red Edge 3, 783 nm) — 20 m
  • B8 (NIR, 842 nm) — 10 m
  • B8A (NIR narrow, 865 nm) — 20 m
  • B11 (SWIR 1, 1610 nm) — 20 m
  • B12 (SWIR 2, 2190 nm) — 20 m

Spectral Indices (4)

  • NDVI = (B8 − B4) / (B8 + B4) — vegetation vigor
  • NDWI = (B3 − B8) / (B3 + B8) — water bodies
  • NDBI = (B11 − B8) / (B11 + B8) — built-up areas
  • NBR = (B8A − B12) / (B8A + B12) — burn scars

Optional Extras

  • GLCM texture (contrast, homogeneity)
  • DEM-derived slope and aspect
  • SAR backscatter (VV, VH)
  • Multi-temporal composites

Random Forest Advantages for LULC

  • ●Handles high-dimensional feature spaces without feature selection (ensemble of decorrelated trees).
  • ●Provides built-in feature importance rankings via mean decrease in impurity (Gini importance) or permutation importance.
  • ●Robust to noisy features and outliers; does not require feature normalization.
  • ●Out-of-bag (OOB) error provides an unbiased estimate of generalization accuracy without a separate validation set.

Hyperparameter Tuning

Key parameters: n_estimators (100–500 trees typically sufficient), max_depth (None for full trees or limit to prevent overfitting), min_samples_leaf (5–10 for pixel-based classification), and max_features($\sqrt{n_{features}}$ is the default for classification).

3. Random Forest Classification Pipeline

This pipeline generates a synthetic multi-class dataset simulating a Sentinel-2 feature stack, trains a Random Forest classifier, computes the full confusion matrix with OA and Kappa, and visualizes feature importances and the classified map.

Random Forest LULC Classification with Accuracy Assessment

Python
script.py193 lines

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Land Cover Change Detection

Change detection identifies pixels that transition between land cover classes over time. The two main approaches are:

Post-Classification Comparison

Classify two dates independently, then compute the transition matrix. Simple but errors in either classification propagate. The change matrix entry $T_{ij}$gives the area that transitioned from class $i$ to class $j$. Total change for class $i$: $\Delta A_i = \sum_j T_{ij} - \sum_j T_{ji}$.

Spectral Change Detection

Compute pixel-wise differences in spectral values or indices between dates. Change Vector Analysis (CVA) uses the magnitude and direction of the spectral change vector: $|\Delta\vec{\rho}| = \sqrt{\sum_b (\rho_{b,t2} - \rho_{b,t1})^2}$. Pixels exceeding a magnitude threshold are flagged as changed.

Area-Adjusted Accuracy (Olofsson et al., 2014)

Standard accuracy metrics assume equal sampling probability across classes. When class proportions are unequal (e.g., 5% urban, 60% forest), area-adjusted estimators correct for sampling bias:

$$\hat{A}_i = A_{total} \sum_j W_j \frac{n_{ji}}{n_{j\cdot}}$$

where $W_j$ is the area proportion of mapped class $j$,$n_{ji}$ is the number of validation samples mapped as class $j$but reference class $i$, and $n_{j\cdot}$ is the total validation samples in mapped class $j$. This provides unbiased area estimates with confidence intervals.

BFAST & Continuous Change Detection

For dense time series (e.g., Sentinel-2 every 5 days), algorithms like BFAST (Breaks For Additive Seasonal and Trend) and CCDC (Continuous Change Detection and Classification) fit harmonic models to the time series and detect breakpoints where the observed values deviate significantly from the model. This enables near-real-time deforestation alerts and gradual change detection.

5. GEE Dynamic World: Urban Expansion (Illustrative)

Google's Dynamic World is a near-real-time LULC product at 10 m resolution, generated for every Sentinel-2 scene using a deep learning model. The following GEE script demonstrates how to detect urban expansion between two periods using Dynamic World classifications.

// ── GEE Dynamic World Urban Expansion Detection ──
// Illustrative code for Google Earth Engine Code Editor

// 1. Define region of interest
var city = ee.Geometry.Point([77.2, 28.6]); // Delhi
var region = city.buffer(30000); // 30km radius

// 2. Load Dynamic World for two periods
var dw_early = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
  .filterBounds(region)
  .filter(ee.Filter.date('2018-01-01', '2018-12-31'))
  .select('label');

var dw_late = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
  .filterBounds(region)
  .filter(ee.Filter.date('2023-01-01', '2023-12-31'))
  .select('label');

// 3. Compute mode (most frequent class) for each period
var mode_early = dw_early.reduce(ee.Reducer.mode()).rename('label');
var mode_late = dw_late.reduce(ee.Reducer.mode()).rename('label');

// Dynamic World classes: 0=water, 1=trees, 2=grass, 3=flooded_veg,
// 4=crops, 5=shrub, 6=built, 7=bare, 8=snow_ice

// 4. Detect urban expansion (non-urban -> urban)
var built_class = 6;
var was_not_built = mode_early.neq(built_class);
var is_built = mode_late.eq(built_class);
var urban_expansion = was_not_built.and(is_built);

// 5. Compute areas by previous land cover
var previous_lulc = mode_early.updateMask(urban_expansion);
var area_by_class = ee.Image.pixelArea()
  .addBands(previous_lulc)
  .reduceRegion({
    reducer: ee.Reducer.sum().group({
      groupField: 1, groupName: 'class'
    }),
    geometry: region,
    scale: 10,
    maxPixels: 1e10
  });
print('Urban expansion by previous class:', area_by_class);

// 6. Compute total urban expansion area
var expansion_area = urban_expansion.selfMask()
  .multiply(ee.Image.pixelArea()).divide(1e6)
  .reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: region,
    scale: 10,
    maxPixels: 1e10
  });
print('Total urban expansion (km2):', expansion_area);

// 7. Visualization
var dwVis = {
  min: 0, max: 8,
  palette: ['419bdf','397d49','88b053','7a87c6','e49635',
            'dfc35a','c4281b','a59b8f','b39fe1']
};

Map.centerObject(region, 11);
Map.addLayer(mode_early.clip(region), dwVis, 'LULC 2018');
Map.addLayer(mode_late.clip(region), dwVis, 'LULC 2023');
Map.addLayer(urban_expansion.selfMask().clip(region),
  {palette: ['ff0000']}, 'Urban Expansion 2018-2023');

// 8. Transition matrix (sampled)
var transition = mode_early.addBands(mode_late)
  .reduceRegion({
    reducer: ee.Reducer.frequencyHistogram(),
    geometry: region,
    scale: 30,
    maxPixels: 1e9
  });
print('Transition frequencies:', transition);

Pipeline Explanation

  • ●Steps 1–2: Load Dynamic World classifications for 2018 and 2023 over the Delhi metropolitan area.
  • ●Step 3: Compute the mode (most frequent class) per pixel across all scenes in each year, reducing temporal noise.
  • ●Step 4: Boolean change detection: pixels that were not built-up in 2018 but are built-up in 2023.
  • ●Steps 5–6: Compute expansion area in km², broken down by what land cover was replaced (cropland, trees, etc.).

6. Global LULC Products

Several global land cover products are available, each with different characteristics:

ESA WorldCover (10 m)

Based on Sentinel-1 and Sentinel-2 data. 11 classes. Global coverage for 2020 and 2021. OA ~75%. The highest resolution global product available.

Google Dynamic World (10 m)

Near-real-time, per-Sentinel-2-scene classification. 9 classes with probability maps. Deep learning model. Available from 2015 to present. Enables temporal analysis.

Copernicus Global Land Cover (100 m)

Discrete and fractional cover maps. 23 classes. Annual products from 2015. Provides fractional cover of trees, shrubs, herbaceous, and bare soil per pixel.

MODIS MCD12Q1 (500 m)

Annual land cover from 2001 to present. Multiple classification schemes (IGBP, UMD, LAI/fPAR). The longest consistent moderate-resolution record available.

7. Key Takeaways

✓Random Forest classification with a 14-feature Sentinel-2 stack (10 bands + NDVI + NDWI + NDBI + NBR) routinely achieves OA > 85% for 6–10 class LULC maps.
✓The Kappa coefficient accounts for chance agreement and provides a more conservative accuracy measure than OA alone.
✓Feature importance from Random Forest reveals which bands and indices contribute most to class separability, guiding feature engineering.
✓Change detection via post-classification comparison or spectral change vectors enables monitoring of urbanization, deforestation, and agricultural expansion.
✓Google Dynamic World provides near-real-time 10 m LULC classification for every Sentinel-2 scene, enabling planetary-scale monitoring.