Comparing different scikit-learn classifiers for land cover mapping

Introduction¶

This lab introduces participants to the application of machine learning algorithms for land cover classification using satellite imagery. Land cover classification is a fundamental task in remote sensing, supporting a wide range of applications including environmental monitoring, urban planning, agriculture, and natural resource management.

In this lab, we explore and compare the performance of four commonly used machine learning classifiers: k-nearest neighbor (KNN), support vector machines (SVM), decision trees, and random forest. These classifiers are implemented using the Scikit-Learn library.

Dataset:

Sentinel-2 image: S2_JFM_2025.tif (Jan–Mar 2025 composite)
Training areas: training_polygons.geojson (7 land cover classes)

Goal¶

By the end of this lab, participants will have a foundational understanding of how to apply and evaluate basic machine learning algorithms for classifying land cover types from Earth observation data.

Imports and Setup¶

Install libraries¶

Install any additional libraries not installed by default.

!pip install rasterio earthpy

Import libraries¶

Import the necessary libraries (pandas, numpy, scikit-learn, rasterio, geopandas, etc.).

import os
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
# import earthpy.plot as ep
from rasterio.mask import mask as rio_mask
from shapely.geometry import mapping
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
import matplotlib.pyplot as plt
from matplotlib.colors import from_levels_and_colors
from matplotlib.patches import Patch
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

OUTPUT_DIR = '../outputs'
os.makedirs(OUTPUT_DIR, exist_ok=True)

Define paths and variables¶

Define the paths to the Sentinel-2 image and training polygon GeoJSON. We also define the spectral bands used as input features, the class labels, and the visualization palette.

# File paths
Image_Path    = 'data/S2_JFM_2025.tif'
Training_Path = 'data/training_polygons.geojson'

# Spectral bands (feature columns) — assumes image bands follow this order
Bands = ['B2', 'B3', 'B4', 'B8']

# Classes use Cl_Id (0-indexed) from training_polygons.geojson
Classes   = [0, 1, 2, 3, 4, 5, 6]
N_Classes = 7
Names     = ['Water', 'Built-up', 'Forest', 'Clouds', 'Bare Soil', 'Pond', 'Crops']

# Distinct palette — each class uses a perceptually separate color
Palette = [
    '#1565C0',  # Water      — dark blue
    '#E53935',  # Built-up   — red
    '#1B5E20',  # Forest     — dark green
    '#CE93D8',  # Clouds     — lavender
    '#8D6E63',  # Bare Soil  — brown
    '#00BCD4',  # Pond       — cyan
    '#F9A825',  # Crops      — amber
]

print('Classes:', list(zip(Classes, Names)))
print('Palette:', Palette)

Load the Sentinel-2 Imagery¶

We use rasterio to open the .tif file and inspect its properties. A true-color composite (B4, B3, B2) is displayed for a quick visual check.

# Open the multiband Sentinel-2 image
image = rasterio.open(Image_Path)

band_count = image.count
height     = image.height
width      = image.width
crs        = image.crs
transform  = image.transform

print(f'Image has {band_count} bands.')
print(f'Image dimensions: {height} rows x {width} columns')
print(f'Coordinate Reference System (CRS): {crs}')

# Validate band count matches expected Bands list
if band_count != len(Bands):
    print(f'Warning: image has {band_count} bands but Bands list has {len(Bands)} entries.')
    print('Update the Bands list above to match the image.')
    Bands = [f'B{i+1}' for i in range(band_count)]  # fallback generic names

# True-color composite (B4=band3, B3=band2, B2=band1 in 1-based indexing)
image_vis = np.stack([image.read(b) for b in [3, 2, 1]])

# plt.figure(figsize=(8, 8))
# ep.plot_rgb(image_vis, stretch=True)
# plt.title('Sentinel-2 True Color Composite (B4, B3, B2)')
# plt.savefig(os.path.join(OUTPUT_DIR, 'S2_true_color.png'), dpi=150, bbox_inches='tight')
# plt.show()

Load and Prepare Training Data¶

Training data is provided as polygon geometries in training_polygons.geojson. We use rasterio.mask to sample pixel values from the Sentinel-2 image within each polygon. Each sampled pixel becomes one training sample, with its spectral band values as features and the polygon’s Cl_Id as the class label.

# Load training polygons
gdf = gpd.read_file(Training_Path)
print(f'Loaded {len(gdf)} training polygons')
print('\nClasses in training data:')
print(gdf[['Cl_Id', 'class_name']].drop_duplicates().sort_values('Cl_Id').to_string(index=False))

# Reproject to match image CRS if needed
if gdf.crs is None or gdf.crs != image.crs:
    gdf = gdf.to_crs(image.crs)
    print(f'\nReprojected training polygons to: {image.crs}')

# Extract pixel values from each polygon
records = []
for _, row in gdf.iterrows():
    geom = [mapping(row.geometry)]
    try:
        out_image, _ = rio_mask(image, geom, crop=True, nodata=0)
        # Reshape from (bands, rows, cols) → (pixels, bands)
        pixels = out_image.reshape(band_count, -1).T
        # Remove nodata pixels (all-zero rows)
        valid = np.any(pixels != 0, axis=1)
        pixels = pixels[valid]
        for pixel in pixels:
            rec = {Bands[i]: float(pixel[i]) for i in range(len(Bands))}
            rec['class'] = int(row['Cl_Id'])
            records.append(rec)
    except Exception as e:
        print(f'Warning: skipped polygon {row.name} — {e}')

df = pd.DataFrame(records)
print(f'Total training samples extracted (before cap): {len(df)}')

# Limit to a maximum number of samples per class
MAX_SAMPLES_PER_CLASS = 200
sampled = [
    group.sample(min(len(group), MAX_SAMPLES_PER_CLASS), random_state=42)
    for _, group in df.groupby('class')
]
df = pd.concat(sampled).reset_index(drop=True)
print(f'Total training samples after capping at {MAX_SAMPLES_PER_CLASS}/class: {len(df)}')
print('\nSamples per class:')
print(df['class'].value_counts().sort_index().rename(lambda x: f'Class {x} ({Names[x]})'))

# Separate features (X) and label (y)
X = df[Bands]
y = df['class']

print(f'\nMissing values in features: {X.isnull().sum().sum()}')
print(f'Missing values in label:    {y.isnull().sum()}')

# Split into training and testing subsets
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
)
print(f'\nTraining samples: {X_train.shape[0]}')
print(f'Testing samples:  {X_test.shape[0]}')

Exploratory Data Analysis (EDA)¶

EDA is a critical initial step in the geospatial ML workflow. It uses statistical summaries and visualizations to gain insights into the dataset before applying any classification algorithm.

Class distribution¶

First, create a count plot to show the distribution of training samples per land cover class.

df_train = X_train.copy()
df_train['class'] = y_train.values

plt.figure(figsize=(10, 6))
sns.countplot(x='class', hue='class', data=df_train, palette=Palette, legend=False)
plt.title('Distribution of Samples Across Classes')
plt.xlabel('Class')
plt.ylabel('Number of Samples')
plt.xticks(ticks=Classes, labels=Names, rotation=15, ha='right')
plt.tight_layout()
plt.savefig(os.path.join(OUTPUT_DIR, 'eda_class_distribution.png'), dpi=150, bbox_inches='tight')
plt.show()

Scatter plot¶

Create a scatter plot of two spectral bands colored by class. This helps identify how well different classes separate in feature space.

plt.figure(figsize=(10, 6))
sns.scatterplot(x='B4', y='B8', hue='class', data=df_train, palette=Palette)
plt.title('Scatter Plot of Band B4 vs. Band B8')
plt.xlabel('Red (B4) reflectance')
plt.ylabel('NIR (B8) reflectance')
plt.legend(title='Class', labels=Names)
plt.tight_layout()
plt.savefig(os.path.join(OUTPUT_DIR, 'eda_scatter_B4_B8.png'), dpi=150, bbox_inches='tight')
plt.show()

Pairwise scatter plots¶

Pairwise scatter plots visualize how each pair of spectral bands relates to one another across classes, helping to spot patterns and assess feature separability.

g = sns.pairplot(df_train, vars=Bands, hue='class', palette=Palette)
g.fig.suptitle('Pairwise Scatter Plots for All Bands', y=1.02)
g.fig.savefig(os.path.join(OUTPUT_DIR, 'eda_pairplot.png'), dpi=100, bbox_inches='tight')
plt.show()

3D scatter plot¶

A 3D scatter plot visualizes three bands simultaneously (Blue B2, Red B4, NIR B8), revealing cluster separation that may not be visible in 2D.

color_map = {cls: Palette[i] for i, cls in enumerate(Classes)}
colors = df_train['class'].map(color_map)

fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(df_train['B2'], df_train['B4'], df_train['B8'],
           c=colors, s=15)
ax.set_xlabel('Blue Band (B2)', labelpad=15)
ax.set_ylabel('Red Band (B4)', labelpad=15)
ax.set_zlabel('NIR Band (B8)', labelpad=15)
ax.dist = 11

legend_elements = [
    Patch(facecolor=Palette[i], edgecolor='k', label=f'{Names[i]}')
    for i in range(N_Classes)
]
ax.legend(handles=legend_elements, title='Classes', loc='upper left',
          bbox_to_anchor=(1.05, 1))
plt.title('3D Scatter Plot: Blue, Red, and NIR Bands', pad=20)
plt.savefig(os.path.join(OUTPUT_DIR, 'eda_3d_scatter.png'), dpi=150, bbox_inches='tight')
plt.show()

t-SNE projection¶

t-SNE reduces the high-dimensional spectral feature space to 2D, helping to visualize cluster structure and class separability across all bands at once.

from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, perplexity=30, random_state=42, max_iter=1000)
X_tsne = tsne.fit_transform(X)

colors_tsne = [color_map[label] for label in y]

plt.figure(figsize=(10, 8))
plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=colors_tsne, s=10, alpha=0.7)
legend_elements = [
    Patch(facecolor=Palette[i], edgecolor='k', label=Names[i])
    for i in range(N_Classes)
]
plt.legend(handles=legend_elements, title='Classes', bbox_to_anchor=(1.05, 1), loc='upper left')
plt.title('t-SNE Projection of Land Cover Samples')
plt.xlabel('t-SNE Dimension 1')
plt.ylabel('t-SNE Dimension 2')
plt.tight_layout()
plt.savefig(os.path.join(OUTPUT_DIR, 'eda_tsne.png'), dpi=150, bbox_inches='tight')
plt.show()

Box plots¶

Box plots show the distribution of reflectance values per band for each land cover class, helping to identify which bands are most discriminative.

df_melt = df_train.melt(
    id_vars='class',
    value_vars=Bands,
    var_name='Band',
    value_name='Reflectance'
)

plt.figure(figsize=(14, 8))
sns.boxplot(x='Band', y='Reflectance', hue='class', data=df_melt, palette=Palette)
plt.title('Box Plots of Reflectance Values by Band and Class')
plt.xlabel('Spectral Band')
plt.ylabel('Reflectance Value')
plt.legend(title='Class', labels=Names, loc='upper right')
plt.tight_layout()
plt.savefig(os.path.join(OUTPUT_DIR, 'eda_boxplots.png'), dpi=150, bbox_inches='tight')
plt.show()

Train and Evaluate Multiple Classifiers¶

We train four classifiers: KNN, SVM, Decision Tree, and Random Forest, then compare their performance using accuracy, confusion matrices, and classification reports.

# 1) KNN
knn = KNeighborsClassifier(n_neighbors=5)
knn.fit(X_train, y_train)
knn_preds = knn.predict(X_test)

# 2) SVM
svm = SVC(kernel='rbf', C=3)
svm.fit(X_train, y_train)
svm_preds = svm.predict(X_test)

# 3) Decision Tree
dt = DecisionTreeClassifier(max_depth=10, random_state=42)
dt.fit(X_train, y_train)
dt_preds = dt.predict(X_test)

# 4) Random Forest
rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)
rf_preds = rf.predict(X_test)

Model validation¶

Compare each model’s performance using accuracy, classification report, and confusion matrix.

models = {
    'KNN':           knn_preds,
    'SVM':           svm_preds,
    'Decision Tree': dt_preds,
    'Random Forest': rf_preds,
}

for name, preds in models.items():
    print(f'\n** {name} **')
    print('Accuracy:', (preds == y_test).mean())
    print(classification_report(y_test, preds, target_names=Names))

    cm = confusion_matrix(y_test, preds, labels=Classes)
    disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=Names)
    fig, ax = plt.subplots(figsize=(9, 7))
    disp.plot(cmap='Blues', ax=ax)
    ax.set_title(f'{name} — Confusion Matrix')
    plt.xticks(rotation=30, ha='right')
    plt.tight_layout()
    plt.savefig(os.path.join(OUTPUT_DIR, f'cm_{name.lower().replace(" ","_")}.png'),
                dpi=150, bbox_inches='tight')
    plt.show()

Land Cover Classification¶

We use the Random Forest classifier (typically the best performer) to predict a class label for every pixel in the Sentinel-2 image.

i. Read the full image stack¶

# Read all bands as a NumPy array — shape: (band_count, height, width)
img_array = image.read()
print('Image array shape:', img_array.shape)

ii. Reshape for prediction¶

Scikit-Learn expects a 2D array where each row is a pixel and each column is a band. Reshape from (bands, height, width) → (height×width, bands).

img_reshaped = img_array.reshape(band_count, -1).T  # (rows*cols, bands)
print('Reshaped array for prediction:', img_reshaped.shape)

iii. Predict and reshape back¶

Run the Random Forest classifier on every pixel and reshape the result back to the original image dimensions.

img_reshaped_df = pd.DataFrame(img_reshaped, columns=Bands)
prediction = rf.predict(img_reshaped_df).astype(np.uint8)
prediction_map = prediction.reshape(height, width)
print('Predicted map shape:', prediction_map.shape)

img_reshaped_df = pd.DataFrame(img_reshaped, columns=Bands)
prediction = svm.predict(img_reshaped_df).astype(np.uint8)
prediction_map = prediction.reshape(height, width)
print('Predicted map shape:', prediction_map.shape)

Visualize the land cover map¶

Display the resulting map with the class color palette.

levels = Classes + [max(Classes) + 1]
cmap, norm = from_levels_and_colors(levels, Palette)

plt.figure(figsize=(12, 10))
im = plt.imshow(prediction_map, cmap=cmap, norm=norm)
plt.title('Land Cover Classification Map\n(Sentinel-2 JFM 2025 — Random Forest)')

cbar = plt.colorbar(im, shrink=0.7)
tick_positions = [i + 0.5 for i in Classes]
cbar.set_ticks(tick_positions)
cbar.set_ticklabels(Names)

plt.tight_layout()
plt.savefig(os.path.join(OUTPUT_DIR, 'LC_map_RF_S2_JFM_2025.png'), dpi=150, bbox_inches='tight')
plt.show()

# Build a valid-pixel mask — skip pixels where all bands are 0 (nodata)
valid_mask = np.any(img_reshaped != 0, axis=1)  # shape: (height*width,)
print(f'Total pixels:  {len(valid_mask):,}')
print(f'Valid pixels:  {valid_mask.sum():,}  ({valid_mask.mean()*100:.1f}%)')

# Predict only on valid pixels (much faster than predicting on the full image)
valid_pixels = pd.DataFrame(img_reshaped[valid_mask], columns=Bands)
valid_preds  = rf.predict(valid_pixels).astype(np.uint8)

# Place predictions back into a full-size array (255 = nodata)
prediction = np.full(len(valid_mask), 255, dtype=np.uint8)
prediction[valid_mask] = valid_preds

prediction_map = prediction.reshape(height, width)
print('Predicted map shape:', prediction_map.shape)

Export the Land Cover Map¶

Save the classification result as a GeoTIFF, preserving the original CRS and spatial transform.

output_tif = os.path.join(OUTPUT_DIR, 'LC_RF_S2_JFM_2025.tif')

with rasterio.open(
    output_tif,
    'w',
    driver='GTiff',
    height=height,
    width=width,
    count=1,
    dtype=prediction_map.dtype,
    crs=crs,
    transform=transform
) as dst:
    dst.write(prediction_map, 1)

print('Classification result exported to:', output_tif)

👩🏽‍💻 TEST IT OUT

What are ways that you can improve the results of your model?