Lesson 4: Supervised Classification

Open slides

Important
  • Try this without AI first — struggle builds understanding.
  • Ask if you are stuck. Asking builds understanding too!

Data: Download the Ces aerial imagery (TIF files for 1961–1995, 2021, 2024) from Moodle and place them in data-large/ces/. The ground truth polygons (data/ces_polygons.gpkg) are already in the repository.

Task 1

Load the 1961 aerial image (data-large/ces/1961.tif). The raw pixel values are in the range [0, 255] — normalize them to [0, 1] by dividing by 255, then rename the layer "luminosity".

Load the ground truth polygons from data/ces_polygons.gpkg. Inspect the available layers with st_layers(). Load both the ces_polygons layer (forest / non-forest polygons) and the study_area layer. Inspect the ces_polygons attribute table — what classes are present?

Plot the image and overlay the polygons with plot(st_geometry(...), add = TRUE).

Task 2

Use st_sample() to draw 1000 random points from the forest/non-forest polygons. Convert the result to an sf object with st_sf(), then use st_join() to attach the class column from ces_polygons. Drop the jahr column and convert class to a factor.

Split the points into a 70 % training set and a 30 % test set. Use set.seed(42) and sample() on the row indices.

How many points ended up in each class in the training set?

Task 3

Extract the pixel values from ces1961 at the training points using terra::extract() (set ID = FALSE). Combine the class column with the extracted values using cbind(), then drop the geometry with st_drop_geometry().

Train a CART classifier with rpart(class ~ luminosity, data = ...). Plot the decision tree with rpart.plot(). What luminosity threshold separates forest from non-forest?

Task 4

Predict land cover for the full raster using predict(ces1961, cart_model, type = "prob"). This returns one probability layer per class. Use which.max() to get a single-layer raster with the most probable class index.

Set the categorical levels on the result using levels():

classes_df <- data.frame(id = 1:2, class = levels(data_train2$class))
levels(pred_class) <- classes_df

Plot the predicted classes. Then evaluate your model:

  1. Extract predictions at the test points with terra::extract()
  2. Build a confusion matrix using table(predicted = ..., actual = ...)
  3. Compute overall accuracy as sum(diag(...)) / sum(...)

Task 5

The single-band classifier uses only brightness to distinguish forest from non-forest. Add a texture feature using a focal standard deviation filter over a 5×5 window:

i <- 5
focal_window <- focalMat(ces1961, d = (i-1)/2 * res(ces1961)[1], type = "rectangle")
ces1961_focal_sd <- focal(ces1961, focal_window, fun = "sd")

Before combining features, normalize the focal layer to [0, 1] using min-max scaling (see slides). Stack it with the original raster and name the layers c("luminosity", "focal_sd").

Extract features from the multi-band raster at the training points, retrain the CART model on both features, predict on the full raster, and compute accuracy on the test set.

How much does the texture feature improve accuracy? Look at the decision tree — does the model actually use focal_sd?

Part 2: RGB Classification (2021)

The 2021 swissimage data covers the same study area in colour (red, green, blue bands at 2 m resolution). Your task is to repeat the classification workflow — this time using RGB instead of grayscale.

The ground truth polygons for 2021 are already prepared in data/ces_polygons.gpkg under the layer ces_polygons_2021 (same Forest / Non-forest structure as before).

Task 6

Load data-large/ces/2021_rgb.tif. The raw values are in [0, 255] — normalise to [0, 1] and name the three bands c("red", "green", "blue").

Load the ces_polygons_2021 layer from data/ces_polygons.gpkg. Sample 1000 points, attach class labels, drop jahr and area, convert class to a factor, then do a 70/30 train/test split with set.seed(42).

How are Forest and Non-forest distributed in the training set? Is the class balance similar to the 1961 split?

Task 7

Extract pixel values from all three bands at the training points. Train a CART model using all available features (class ~ .). Plot the decision tree — which band(s) does the model use?

Predict on the full 2021 raster, set categorical levels, and plot the result. Evaluate on the test set with a confusion matrix and overall accuracy.

Tip

The classes_df lookup should list levels in the same order R assigns them — alphabetical by default. Use levels(train_df$class) to check.

Task 8

Compare the two classifications to quantify forest change between 1961 and 2021.

Step 1 — Compute total forest area (in hectares) for each year. The pixel area in m² is prod(res(raster)).

Step 2 — Create a change map. The two rasters have different resolutions (0.5 m vs 2 m), so resample the 1961 result to 2 m first with resample(..., method = "near"). Then classify each pixel as: stable forest, stable non-forest, forest gained, or forest lost.

Tip

terra::ifel(condition, yes, no) is the raster equivalent of ifelse(). You can nest calls to encode multiple categories.

What does the change map tell you? Where did forest gain or loss occur?