library(terra)

# Common prefix
prefix <- "S2A_OPER_MSI_L1C_TL_MPS__20160629T132555_A005324_T34UEV_"

# Build filenames for the needed bands
f_b2 <- paste0(prefix, "B02.tif")   # Blue
f_b3 <- paste0(prefix, "B03.tif")   # Green
f_b4 <- paste0(prefix, "B04.tif")   # Red
f_b8 <- paste0(prefix, "B08.tif")   # NIR

# Read bands
b2 <- rast(f_b2)
b3 <- rast(f_b3)
b4 <- rast(f_b4)
b8 <- rast(f_b8)

# If needed, align to the first band
b3 <- resample(b3, b2)
b4 <- resample(b4, b2)
b8 <- resample(b8, b2)

# Create multiband raster
ms <- c(b2, b3, b4, b8)
names(ms) <- c("Blue", "Green", "Red", "NIR")

# Inspect
print(ms)

# Plot individual bands
par(mfrow = c(2, 2))
plot(ms[[1]], main = "Blue")
plot(ms[[2]], main = "Green")
plot(ms[[3]], main = "Red")
plot(ms[[4]], main = "NIR")

# True colour composite
plotRGB(ms, r = 3, g = 2, b = 1, stretch = "lin",
        main = "True colour composite")

# False colour composite
plotRGB(ms, r = 4, g = 3, b = 2, stretch = "lin",
        main = "False colour composite")

# NDVI
ndvi <- (ms[[4]] - ms[[3]]) / (ms[[4]] + ms[[3]])
plot(ndvi, main = "NDVI")

# Save outputs
writeRaster(ms, "sentinel_stack.tif", overwrite = TRUE)
writeRaster(ndvi, "ndvi.tif", overwrite = TRUE)

------------------------------------
27.03.2026

# =========================================================
# Random Forest classification from sample polygons
# Sentinel-2 bands: B2, B3, B4, B8 (10 m)
# Output: classified raster + confusion matrix
# =========================================================

# Install packages if needed:
# install.packages(c("terra", "sf", "randomForest", "caret", "dplyr"))

library(terra)
library(sf)
library(randomForest)
library(caret)
library(dplyr)

# -----------------------------
# 1. INPUTS
# -----------------------------
# Raster with 4 bands: B2, B3, B4, B8
raster_path <- "path/to/sentinel_B2_B3_B4_B8.tif"

# Training polygons
shp_path <- "path/to/UA_clip_sample_polygons.shp"

# Name of the attribute containing class labels
# Change this to the real field in your shapefile, e.g. "class", "Class", "landcover"
class_field <- "class"

# Output files
classified_raster_path <- "rf_classification.tif"
model_rds_path <- "rf_model.rds"

# -----------------------------
# 2. LOAD DATA
# -----------------------------
img <- rast(raster_path)
polygons <- st_read(shp_path, quiet = TRUE)

# Check geometry and CRS
print(img)
print(polygons)

# Reproject polygons if needed
if (st_crs(polygons)$wkt != crs(img)) {
  polygons <- st_transform(polygons, crs(img))
}

# Convert sf to terra vector
polygons_v <- vect(polygons)

# Check if class field exists
if (!(class_field %in% names(polygons_v))) {
  stop(
    paste0(
      "Field '", class_field, "' not found in shapefile.\n",
      "Available fields: ", paste(names(polygons_v), collapse = ", ")
    )
  )
}

# Make sure class labels are factors
polygons_v[[class_field]] <- as.factor(polygons_v[[class_field]])

# Optional: set band names if needed
if (nlyr(img) == 4) {
  names(img) <- c("B2", "B3", "B4", "B8")
}

# -----------------------------
# 3. EXTRACT TRAINING SAMPLES
# -----------------------------
# Extract pixel values under polygons
# touches=TRUE includes pixels touched by polygon borders
samples <- extract(img, polygons_v, df = TRUE, touches = TRUE)

# 'extract' returns ID linking rows to polygons
# Join the class label from polygons
poly_df <- as.data.frame(polygons_v)
poly_df$ID <- 1:nrow(poly_df)

samples <- samples %>%
  left_join(poly_df[, c("ID", class_field)], by = "ID")

# Remove missing values
samples <- na.omit(samples)

# Keep only predictor bands + class
samples <- samples %>%
  select(all_of(names(img)), all_of(class_field))

# Convert class to factor
samples[[class_field]] <- as.factor(samples[[class_field]])

# Inspect class balance
print(table(samples[[class_field]]))

# -----------------------------
# 4. TRAIN / TEST SPLIT
# -----------------------------
set.seed(123)

train_index <- createDataPartition(samples[[class_field]], p = 0.7, list = FALSE)
train_data <- samples[train_index, ]
test_data  <- samples[-train_index, ]

# -----------------------------
# 5. TRAIN RANDOM FOREST MODEL
# -----------------------------
rf_model <- randomForest(
  x = train_data[, names(img)],
  y = train_data[[class_field]],
  ntree = 500,
  importance = TRUE
)

print(rf_model)
print(importance(rf_model))
varImpPlot(rf_model)

saveRDS(rf_model, model_rds_path)

# -----------------------------
# 6. PREDICT TEST SET
# -----------------------------
test_pred <- predict(rf_model, newdata = test_data[, names(img)])

# Confusion matrix
cm <- confusionMatrix(test_pred, test_data[[class_field]])
print(cm)

# If you want only the matrix:
print(cm$table)

# Overall accuracy and kappa
cat("\nOverall Accuracy:", cm$overall["Accuracy"], "\n")
cat("Kappa:", cm$overall["Kappa"], "\n")

# -----------------------------
# 7. APPLY MODEL TO FULL RASTER
# -----------------------------
classified <- predict(img, rf_model, type = "response", na.rm = TRUE)

# Save output
writeRaster(classified, classified_raster_path, overwrite = TRUE)

# Plot result
plot(classified, main = "Random Forest Classification")