📝 ACTIVIDAD DE EVALUACIÓN 1 CLASIFICACIÓN SUPERVISADA:

Esta actividad se realiza sobre la ciudad canadiense de Calgary. Los datos necesarios están incluidos en el fichero calgary que contiene:

Banda Denominación
Banda 2 Blue
Banda 3 Green
Banda 4 Red
Banda 5 NIR
Banda 6 SWIR 1
Banda 7 SWIR 2

Esta práctica tiene por objeto obtener un mapa temático aplicando el algoritmo CART a partir del uso de puntos de entrenamiento.

► Cargar los paquetes necesarios para la actividad

library(terra)
library(tidyverse)
library(RStoolbox)
library(sf)
library(rpart)
library(rpart.plot)

► Cargar la escena multibanda.

# Obtén la lista de archivos .tif en tu carpeta
archivos <- list.files(path = "D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/actividad_evaluacion/supervisada/ejercicio_calgary/", 
                       pattern = "\\.tif$", 
                       full.names = TRUE)
imagen <- rast(archivos)
print(imagen)                     # Opcional: Revisa que el número de capas sea correcto

► Transformar los valores de reflectancia superficial de números digitales a reflectancia verdadera.

imagen <- (imagen * 0.0000275) - 0.2

► Cargar el fichero de entrenamiento y verificar su contenido

vect_puntos_entrenamiento <- vect("D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/actividad_evaluacion/supervisada/ejercicio_calgary/calgary_puntos_entrenamiento.geojson")
print(vect_puntos_entrenamiento)                                            # Verifica el contenido

► Compruebar la coincidenciaentre el sistema de coordenadas de referencia del objeto ráster y del objeto vectorial.

vect_puntos_entrenamiento <- project(vect_puntos_entrenamiento, crs(imagen))

► Representar la escena multibanda en color verdadero y superpuestos, los puntos de entrenamiento en color amarillo.

# Representar en color verdadero.
plotRGB(imagen, r = 3, g = 2, b = 1, stretch = "lin")

# Superponer los puntos en amarillo
plot(vect_puntos_entrenamiento, add = TRUE, col = "yellow", cex = 0.8, pch = 16)

► Extraer los valores espectrales correspondientes a los puntos de entrenamiento

df_muestra <- terra:: extract(imagen,                                      # Imagen de satélite
                              vect_puntos_entrenamiento, 
                              ID = FALSE)    

df_muestra <- cbind(as.data.frame(vect_puntos_entrenamiento), df_muestra)    # Usamos 'as.data.frame' para obtener los atributos del vector y combinarlos
df_muestra$id <- NULL         # 'extract' añade una columna 'ID' que suele coincidir con el orden de las filas
head(df_muestra)

► Dividir el dataframe en dos subconjuntos, uno de entrenamiento (con el 75 % de los casos) y otro de verificación (con el 25 % restante).

df_entrenamiento <- df_muestra %>% 
  sample_frac(0.7)                                # Crear el set de entrenamiento (70%)

df_verificacion <- df_muestra %>% 
  anti_join(df_entrenamiento)          # Crear el set de verificación con el resto (30%). Usamos anti_join para asegurar que no se repitan filas

► Aplicar el algoritmo CART a los datos de entrenamiento.

modelo_cart <- rpart(as.factor(uso)~., 
                   data = df_entrenamiento, 
                   method = 'class', 
                   control = rpart.control(cp = 0.0, minsplit = 2))

► Representar gráficamente el árbol de decisión

rpart.plot(modelo_cart, 
           box.palette = 0, 
           main = "Árbol de clasificación")

► Verificar el modelo con el 25% restante

prediccion_test <- predict(modelo_cart, 
                           df_verificacion, 
                           type = "class")

► Evaluar el resultado de la clasificación creando una matriz de confusión en la que se compare la clasificación obtenida con los datos de entrenamiento y la obtenida con los datos de verificacion

matriz_confusion <- table(Real = df_verificacion$uso, 
                          Predicho = prediccion_test)
print(matriz_confusion)

► Calcular la precisión Global (Accuracy)

accuracy <- sum(diag(matriz_confusion)) / sum(matriz_confusion)
print(paste("Precisión Global:", round(accuracy, 4)))

n <- sum(matriz_confusion)                 # Primero se calcula el número total de casos/muestras
n

diag <- diag(matriz_confusion)             # Número de casos correctamente clasificados por clase
OA <- sum(diag) / n                        # Precisión general
OA

► Calcular el valor de Kappa:

rowsums <- apply(matriz_confusion, 1, sum)                   # Casos observados (verdaderos) por clase
p <- rowsums / n

colsums <- apply(matriz_confusion, 2, sum)                    # Casos predichos por clase
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa

► Clasificar todas las celdas de la imagen original.

raster_cart <- predict(imagen, 
                      modelo_cart, 
                      na.rm = TRUE)
raster_cart
plot(raster_cart)

► Representa gráficamente en un mapa temático el uso de suelo más probable.

raster_cart_final <- which.max(raster_cart)
raster_cart_final

raster_cart_final <- as.factor(raster_cart_final)                             # Conversión en factor

levels(raster_cart_final) <- data.frame(ID = 1:4,
                          clase = c("nubes","artificial","natural", "agua"))  # Creación de las etiquetas

colores_clasificacion_cart <- c("nubes" = "grey", "artificial" = "red", "natural" = "darkgreen", "agua" = "blue", "urbano" = "red")
levels(raster_cart_final)

ggR(raster_cart_final, geom_raster = TRUE, stretch = "lin", forceCat = TRUE) + 
  ggtitle("Clasificación no supervisada. Método CART") +                     # Título
  labs(x ="Longitud (m)", y ="Latitud (m)") +                                   # Etiquetas de los ejes
  theme(plot.title = element_text(hjust =0.5,                                   # Alineación del título 
                                  size =20),                                    # Tamaño del título
        axis.title = element_text(size =10),                                    # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),                                          # Tamaño de la leyenda
        legend.title = element_text(size = 20),                                 # Tamaño del título de la leyenda
        legend.text = element_text (size = 10)) +                               # Tamaño del texto de la leyenda
   # Dibuja los colores según un gradiente 
scale_fill_manual(values = colores_clasificacion_cart,  
                  name = "Clase") +
  theme_minimal()

► Crear un mapa representando sólo con las superficies naturales

solo_natural <- ifel(raster_cart_final == 3, 1, NA)
plotRGB(imagen, r=4, g=3, b=2, axes=TRUE, stretch="lin")           # Dibuja la imagen RGB
plot(solo_natural, col = "green", add = TRUE, legend = FALSE)            # Superponer la categoría 1

📝 ACTIVIDAD DE EVALUACIÓN 2 CLASIFICACIÓN SUPERVISADA:

Esta actividad consiste en la realización de un procedimiento de clasificación supervisada a partir de puntos de entrenamiento extraidos de una base de datos con usos de suelo (por ejemplo, CORINE, Urban Atlas etc… Los datos necesarios están incluidos en el fichero etiopia que contiene:

Banda Denominación
Banda 2 Blue
Banda 3 Green
Banda 4 Red
Banda 5 Red Edge 1
Banda 6 Red Edge 2
Banda 7 Red Edge 3
Banda 8A NIR
Banda 11 SWIR 1
Banda 12 SWIR 2

En este caso, las categorías de usos de suelo son los siguientes:

ID usos
1 cultivos
2 bambú
3 suelo desnudo
4 plantaciones
5 bosque
6 humedales

El objetivo, como en el caso anterior, es obtener un mapa de usos de suelo de naturaleza temática, pero en este caso, los puntos de entrenamiento se conseguirán a partir del fichero ráster.

► Cargar los paquetes necesarios para la actividad.

library(tidyverse)           # Para el manejo de bases de datos
library(terra)               # Para el procesamiento de imágenes en formato ráster
library(RStoolbox)           # Para la visualización de imágenes en formato ráster, con ggplot2
library(ggplot2)             # Para la visualización de imágenes en formato ráster
library(doParallel)          # Para el procesamiento en paralelo
library(caret)               # Para la clasificación
library(kernlab)             # Para implementar SVM

► Importar el fichero de usos de suelo, cambiando el nombre de la variable lulcGewata por el nombre categoria.

raster_usos <- rast("D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/Ejercicio_clasificacion_supervisada/wur/data/usos_lulc.tif")
raster_usos
names(raster_usos) <- "categoria"
unique(raster_usos)

► Asignar los correspondientes nombres a las categorías del fichero ráster (no es categórico aún).

categoria <- c("cultivos", "bambú", "suelo desnudo", "plantaciones", "bosque", "humedales")
df <- data.frame(value = c(1,2,3,4,5,6), 
                 names = categoria)
levels(raster_usos) <- df
is.factor(raster_usos)

► Proporcionar los siguientes colores (como códigos hexidecimales) para su representación gráfica

Categoria colores
cultivos #E6E600
bambú “#81AF32
suelo desnudo #D2B48C
plantaciones “#73FFDF
bosque #006400
humedales #4169E1 
colores_raster_usos <- c("cultivos" = "#E6E600", # Amarillo intenso (estándar para agricultura)
                         "bambú"  = "#81AF32", # Verde oliva claro (distinto al bosque denso)
                         "suelo desnudo" = "#D2B48C", # Canela o ocre (tan)
                         "plantaciones" = "#73FFDF", # Verde aguamarina o verde brillante
                         "bosque" = "#006400", # Verde oscuro (DarkGreen)
                         "humedales" = "#4169E1")  # Azul real (RoyalBlue para zonas inundables)

plot(raster_usos, col = colores_raster_usos)

► Generar una muestra estratificada de 25 puntos de entrenamiento en cada categoría. Representa gráficamente estos puntos sobre el mapa de usos de suelo.

set.seed(123)
puntos_entrenamiento_raster <- spatSample(raster_usos, 
                                          size = 25, 
                                          method="stratified", 
                                          as.points=TRUE)                 # Al trabajar sobre un ráster es necesario señalar que se crean puntos
plot(raster_usos, col = colores_raster_usos)
points(puntos_entrenamiento_raster, col = "red", cex = 1)

head(puntos_entrenamiento_raster)

► Grabar estos puntos de entrenamiento como fichero vectorial:

writeVector(puntos_entrenamiento_raster, "D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/Ejercicio_clasificacion_supervisada/wur/data/puntos_entrenamiento.geojson", overwrite = TRUE)

► Cargar el gránulo de Sentinel 2 denominado S2B2A_T36NZP_20201227T075239_20m.tif

imagen <- rast("D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/Ejercicio_clasificacion_supervisada/wur/data/S2B2A_T36NZP_20201227T075239_20m.tif")
names(imagen)
plot(imagen)

► Extracción de los valores de reflectancia de todas las bandas correspondientes a los puntos de entrenamiento y comprobación de la estructura del dataframe.

df_reflectividad_pixeles <- terra:: extract(imagen,                                      # Imagen de satélite
                                             puntos_entrenamiento_raster,                # Fichero con los puntos de entrenamiento
                                             ID=FALSE)    
head(df_reflectividad_pixeles)

► Creación de otro dataframe combinando

df_muestra <- data.frame(categoria = puntos_entrenamiento_raster$names, 
                         df_reflectividad_pixeles)
df_muestra

► Grabación de este dataframe con formato *.csv para su uso posterior

write.csv2(df_muestra, "D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/Ejercicio_clasificacion_supervisada/wur/data/df_muestra.csv")

► Cálculo de los valores medios de reflectancia para cada clase y cada banda.

df_muestra_promedios <- aggregate(df_muestra[    ,-1],                        # Omite la primera columna, con las categorías
                                  list(puntos_entrenamiento_raster$names), 
                                  mean)

► Renombrar la primera columna del dataframe como categoria.

df_muestra_promedios <- df_muestra_promedios %>% 
  rename(categoria = Group.1)
df_muestra_promedios

► Pasar a formato largo (necesario para ggplot)

df_long <- df_muestra_promedios %>%
  pivot_longer(cols = -categoria, names_to = "banda", values_to = "reflectancia")

► Asegurar el orden de las bandas en el eje X (si no, R las pondrá en orden alfabético)

df_long$banda <- factor(df_long$banda, 
                         levels = c("B02", "B03", "B04", "B05", "B06", "B07", "B8A", "B11", "B12"))

► Crear el gráfico

ggplot(df_long, aes(x = banda, y = reflectancia, 
                    color = factor(categoria), 
                    group = categoria)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  scale_color_manual(
    values = c("1" = "#E6E600", 
               "2" = "#81AF32", 
               "3" = "#D2B48C", 
               "4" = "#73FFDF", 
               "5" = "#006400", 
               "6" = "#4169E1"),
    labels = c("1" = "cultivos", 
               "2" = "bambú", 
               "3" = "suelo desnudo", 
               "4" = "plantaciones", 
               "5" = "bosque", 
               "6" = "humedales")
  ) +
  labs(title = "Perfiles Espectrales por Uso de Suelo",
       x = "Bandas", y = "Reflectancia",
       color = "Categoría") + 
  theme_bw() + 
  theme(
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 1.5),
    panel.grid.major = element_line(color = "black", linetype = "dashed", linewidth = 0.2),
    panel.grid.minor = element_blank(),
    axis.title.x = element_text(color = "black", face = "bold", size = 14),
    axis.title.y = element_text(color = "black", face = "bold", size = 14),
    axis.text.x = element_text(color = "black", face = "bold", size = 12),
    axis.text.y = element_text(color = "black", face = "bold", size = 12),
    axis.ticks = element_line(color = "black", linewidth = 1),
    axis.ticks.length = unit(0.3, "cm")
  )

► Cargar el dataframe con la muestra de puntos de entrenamiento

df_muestra <- read.csv2("D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/Ejercicio_clasificacion_supervisada/wur/data/df_muestra.csv", row.names = 1)
df_muestra

► División del dataframe con la muestra en dos subconjuntos, uno de entrenamiento (con el 70 % de los casos) y otro de verificación (con el 30% restante).

# Crear el set de entrenamiento (70%)
df_entrenamiento <- df_muestra %>% 
  sample_frac(0.7)

# Crear el set de verificación con el resto (30%). Usamos anti_join para asegurar que no se repitan filas
df_verificacion <- df_muestra %>% 
  anti_join(df_entrenamiento)

► A partir de aquí, cada alumno elegirá si utiliza el algoritmo Random Forest o el SVM y procederá a efectuar todos los cálculos necesarios para elaborar un mapa temático.

# Procesamiento en paralelo
mc <- makeCluster(detectCores())       # Pregunta al sistema operativo el número de núcleos y crea "copias" que se ejecutan en segundo plano.
registerDoParallel(mc)                # "conecta" el equipo que acabas de crear con las funciones que van a realizar el trabajo pesado

# Ajuste de los parámetros para el entrenamiento
tune_par <- trainControl(method = "repeatedcv", 
                         number = 10, 
                         repeats = 5, 
                         allowParallel=TRUE )


# Entrenamiento del modelo SVM

df_entrenamiento$categoria <- as.factor(df_entrenamiento$categoria)             # Convertir la categoría a factor

svm_model <- train(categoria ~ ., 
                   data = df_entrenamiento, 
                   method = "svmRadial", 
                   metric = "Accuracy", 
                   trControl = tune_par)
svm_model

# Detención del procesamiento en paralelo.
stopCluster(mc)

# Predicción de los valores sobre los datos de entrenamiento
predicciones_entrenamiento <- predict(svm_model, 
                                      df_entrenamiento)

# Evaluación de la precisión
# Asegurar que la referencia sea factor
df_entrenamiento$categoria <- as.factor(df_entrenamiento$categoria)

# Forzar a las predicciones a ser factor con los MISMOS niveles que la referencia
predicciones_entrenamiento <- factor(predicciones_entrenamiento, 
                                     levels = levels(df_entrenamiento$categoria))

confusionMatrix(predicciones_entrenamiento, df_entrenamiento$categoria)

# Validación del modelo SVM con el conjunto de datos de verificación 
predicciones_verificacion <- predict(svm_model, 
                                     df_verificacion)

# Para la evaluación de la precisión en el conjunto de datos de verificación
# Asegurar que la referencia sea factor
df_verificacion$categoria <- as.factor(df_verificacion$categoria)

# Forzar a las predicciones a ser factor con los MISMOS niveles que la referencia
predicciones_verificacion <- factor(predicciones_verificacion, 
                                     levels = levels(df_verificacion$categoria))

confusionMatrix(predicciones_verificacion, df_verificacion$categoria)

# Clasificación de la imagen de satélite
raster_svm <- predict(imagen, svm_model)
raster_svm                                             # Muestra los atributos

# Representación gráfica

# Conversión del ráster en factor para que reconozca las etiquetas (asumiendo que los niveles están en el mismo orden en `df_entrenamiento`.
clases_nombres <- levels(as.factor(df_entrenamiento$categoria))
levels(raster_svm) <- data.frame(ID = 1:length(clases_nombres), categoria = clases_nombres)

ggR(raster_svm, geom_raster = TRUE) + 
  ggtitle("Ejemplo de clasificación con SVM") + 
  labs(x = "Longitud", y = "Latitud") + 
  theme(plot.title = element_text(hjust = 0.5, size = 24),
        axis.title = element_text(size = 20),
        legend.key.size = unit(1, "cm"),
        legend.title = element_text(size = 20),
        legend.text = element_text(size = 15)) +
  # Usamos scale_fill_manual para asignar colores fijos a categorías de texto
  scale_fill_manual(values = c("grey", "blue", "wheat", "darkgreen", "red", "skyblue"), 
                    name = "Clase", 
                    labels = clases_nombres)