# 9d Land Cover Classification Practical
# Random Forest Classification

# First, we want to load all necessary libraries, before we start the actual work
library(terra)    # For raster data handling
library(sf)       # For spatial data handling
library(randomForest)  # Random Forest model
library(ggplot2) # nice visuals

# Step 1: Load Raster Data, make sure to change the file path on your PC!

# Set the path to the folder with your Sentinel-2 bands
folder_path <- "C:/Users/schul/Documents/Berufliches/Hiwi EOCap4Africa/Exercices Finalised/Geodata/9d Land Cover Classification Practical/Data/QGIS_data"  # change this to your actual path

# List only .tif files in that folder
band_files <- list.files(folder_path, pattern = "\\.tiff$", full.names = TRUE)

stacked <- rast(band_files)  # this works directly!

# Save to one multiband GeoTIFF
writeRaster(stacked, "sentinel2_stacked.tif", overwrite = TRUE)
sentinel_raster <- stacked  # keep the original raster object

# Visually inspect the Raster Image
plot(sentinel_raster)

# Step 2: Load Training Data (Shapefile), make sure to change the file path on your PC!
training_path <- "C:/Users/schul/Documents/Berufliches/Hiwi EOCap4Africa/Exercices Finalised/Geodata/9d Land Cover Classification Practical/Data/Training_data.shp"
training_data <- st_read(training_path) # This should output some information about the File


# Step 3: Transforming the Training Data
# Random Forest Models in R want Points as inputs, not Polygons, however we can simply generate points within Polygons to give us a larger training set

# First we need to change the Class attribute into a Factor
training_data$Class <- as.factor(training_data$Class)

# visually inspect out Data
ggplot(training_data) +
  geom_sf(aes(fill = Class)) +
  theme_minimal() +
  guides(fill = guide_legend(keywidth = 2))  # encountered a weird bug, this prevents it


# Convert Sf data into SpactVector
training_data_sv <- vect(training_data)

# Set seed for reproducibility
set.seed(1)

# Generate random points inside polygons (1000 per class, this leads to more accuracy)
training_points <- spatSample(training_data_sv, 1000, method="random")


# Step 4: Extract Training Samples
# Extract raster values at the sampled points
training_samples <- extract(sentinel_raster, training_points)

# Add land cover labels from polygons
training_samples$Class <- training_points$Class # make sure your Class column has the same name!


# Step 5: Prepare Data for Modeling
# Convert to dataframe
training_df <- as.data.frame(training_samples)

# Clean up the dataframe
colnames(training_df) <- make.names(colnames(training_df))

# Lets inspect the Dataframe!
print(training_df)

# We want to remove "ID" column from training data before training,  we do not need it and it can cause issues later on
training_df <- training_df[, !(colnames(training_df) %in% c("ID"))]


# Ensure labels are factors
training_df$Class <- as.factor(training_df$Class)


# Step 6: Training the Model
# Upping the number for ntree can yield more accurate results
rf_model <- randomForest(Class ~ ., data = training_df, ntree = 500, importance = TRUE)

# Check model output
print(rf_model)


# Step 7: Apply the trained model to the raster
# Convert raster to dataframe for prediction
raster_df <- as.data.frame(sentinel_raster, xy = TRUE)  # Keep X, Y coordinates

# remove NA values for accurate Prediction
raster_df <- na.omit(raster_df)
# Clean the column names of the raster dataframe to match the training data
colnames(raster_df) <- make.names(colnames(raster_df))

# Now, we can predict the land cover classes for each pixel
predicted_classes <- predict(rf_model, newdata = raster_df)


# Step 8: Apply Predictions to the Raster
# Add predictions back to the dataframe
raster_df$Class <- predicted_classes

# Convert predictions to raster format
classified_raster <- sentinel_raster[[1]]  # Use the original raster as a template
values(classified_raster) <- as.numeric(raster_df$Class)


# Step 9: Visual Insepction of Results
# Get the factor levels (class names), that way our map will have the proper labels and not just numbers for classes
class_names <- levels(training_df$Class)

# Print to check the order of class names
print(class_names)

# Convert raster values to factor
classified_raster <- as.factor(classified_raster)

# Assign class labels to raster categories
levels(classified_raster) <- data.frame(ID = 1:length(class_names), Name = class_names)

# Print raster levels to check mapping
levels(classified_raster)

# Plot raster with categorical legend
plot(classified_raster, main = "Random Forest Land Cover Classification")

# Step 10: Safe your Classified Raster
# Define file path of your choice
output_path <- "C:/Users/schul/Documents/Berufliches/Hiwi EOCap4Africa/Exercices Finalised/Geodata/9d Land Cover Classification Practical/LCC_labeled.tif"

# Save raster with specified path
writeRaster(classified_raster, output_path, overwrite = TRUE)

