Clasificación Supervisada

Introducción

Con el siguiente ejercicio se pretende aplicar algoritmos de clasificación supervisada a una imagen raster, con la finalidad de entender cómo funcionan estos algoritmos.

En un ejercicio anterior usé algoritmos de clasificación no supervisada como K-Means y Clara, a diferencia de ese tipo de clasificación, en la clasificación supervisada se determinan unas características iniciales que serán útiles a la hora de elaborar mi modelo. Otra característica importante es que en con K-Means y Clara se busca hacer agrupación de datos (clustering), con la clasificación supervisada se va a obtener que característica tiene más peso para resolver mi problema.

Metodología

Para la clasificación supervisada se hace un flujograma de trabajo:

  • Se definen unas clases temáticas
  • Clasificación de áreas de entrenamiento adecuadas
  • Calibración del algoritmo de clasificación
  • Clasificación de la imagen
  • Evaluación y verificación de resultados

Usualmente con la clasificación supervisada, se puede obtener una señal espectral con las bandas de imágenes multi-espectrales aplicando una serie de reglas para cada clase como el radio de las bandas, los índices espectrales, las características texturales, las variaciones temporales entre otras características. Estos datos sirven para entrenar un modelo que permite clasificar toda la imagen más allá del espacio de trabajo.

En este taller usaré Random Forest (RF), un algoritmo incluido como un paquete de R, este algoritmo permite hacer una clasificación mediante la construcción de árboles de decisión (de ahí su nombre) durante la etapa de entrenamiento, lo que nos permitirá llegar a una clasificación que además nos permita identificar que variables tienen mayor peso a la hora de realizar la clasificación.

Este ejercicio se realizó con imágenes de Sentinel-2a (S2)

legBerlin <- read.csv(url("https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/legend_berlin.csv"))

knitr::kable(legBerlin)
Type Code Land.cover.class.description
Artificial 1 Compact high-rise
Artificial 2 Compact midrise
Artificial 3 Compact low-rise
Artificial 4 Open high-rise
Artificial 5 Open midrise
Artificial 6 Open low-rise
Artificial 7 Lightweight low-rise
Artificial 8 Large low-rise
Artificial 9 Sparsely built
Artificial 10 Heavy industry
Vegetated 11 Dense trees
Vegetated 12 Scattered trees
Vegetated 13 Bush and scrub
Vegetated 14 Low plants
Vegetated 15 Bare rock or paved
Vegetated 16 Bare soil or sand
Vegetated 17 Water

Para este ejercicio igual que en el anterior vamos a definir un directorio de trabajo para eso usamos el código dentro del siguiente Chunk:

getwd()
## [1] "D:/R-Scrips"
setwd('D:/R-Scrips/') ## Establece un directorio de trabajo

## Posteriormente se crea una carpeta llamada "data-raw" dentro del directorio de trabajo, aquí descargaré la imagen para mi ejercicio

if(!dir.exists("./data-raw")) dir.create("./data-raw")

## Ahora se descarga la imagen con la que voy a trabajar y se descomprimirá en el directorio que creé

download.file("https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/berlin.zip", "./data-raw/berlin.zip", method = "auto") #Descarga un archivo de una dirección específica y lo pone en la carpeta que se indica

unzip("./data-raw/berlin.zip", exdir = "./data-raw") #Descomprime un archivo en, la ruta especificada

Ahora que ya tengo los datos con los que voy a trabajar instalo las librerías necesarias para poder trabajar y se cargan para poder trabajar con estas

#install.packages('raster')
#install.packages('randomForest')
library(raster)
## Loading required package: sp
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.

Ahora que tengo los archivos y las librerías necesarias voy a crear una rasterStack (un apilamiento), además voy a cambiar el nombre de las capas por un nombre más conveniente:

fl <- list.files("./data-raw/berlin/S2", pattern = ".tif$", full.names = TRUE)

rst <- stack(fl)
names(rst) <- c(paste("b",2:8,sep=""),"b8a","b11","b12")

Ahora voy a usar la función plotRGB para visualizar la imagen que descargué, con esta función puedo combinar las capas de un raster en un orden diferente, dependiendo de lo que deseo resaltar, en esta caso quiero ver la imagen en su color natural para esto usamos la combinación 4-3-2, donde 4 es el rojo, el 3 la banda verde y la 2 la azul, a la hora de hacer la combinación hay que tener en cuanta el orden en que se hace el stack y como se renombramos las capas.

plotRGB(rst, r=3, g=2, b=1, scale=1E5, stretch="lin")

Otro ejemplo de combinación es el infrarrojo para resaltar la vegetación, para eso combinamos las bandas 8-4-3 en el orden: r7m g3 y b2

plotRGB(rst, r=7, g=4, b=2, scale=1E5, stretch="lin")

Por último radiación atmosférica (12-11-8A):

plotRGB(rst, r=10, g=9, b=8, scale=1E5, stretch="lin")

Ahora cargamos los ejemplos de entrenamiento usados en la calibración, estos datos sirven como un ejemplo de cómo el algoritmo de clasificación “aprenderá” de la respuesta espectral de cada una de las clases

rstTrain <- raster("./data-raw/berlin/train/berlin_lcz_GT.tif")

rstTrain[rstTrain==0] <- NA # Remover valores nulos del set de entrenamiento (background NA)

rstTrain <- ratify(rstTrain) # Convertir a una representación del tipo discreto

names(rstTrain) <- "trainClass" # Cambiar el nombre de la capa

plot(rstTrain, main="Áreas de Entrenamiento por Clases") # Mostrar el raster de entrenamiento

Ahora veremos cuantos pixeles de entrenamiento hay por cada una de las 12 clases

tab <- table(values(rstTrain))
print(tab)
## 
##    2    4    5    6    8    9   11   12   13   14   16   17 
## 1534  577 2448 4010 1654  761 4960 1028 1050 4424  359 1732

Aunque no siempre sea la mejor aproximación, voy a convertir mi raster data-set en un data.frame, para esto podemos usar RF. Sin embargo, hay que tener en cuenta que en algunos casos, dependiendo del tamaño del RasterStack y la cantidad de memoria disponible podría no ser posible la conversión. Una posible solución a este inconveniente podría ser convertir el raster de entrenamiento en un objeto SpatialPoints, luego correr la función extract.

A continuación para obtener los valores de pixel dentro del dataframe de calibración usamos el siguiente código

rstDF <- na.omit(values(stack(rstTrain, rst)))
rstDF[,"trainClass"] <- as.factor(as.character(rstDF[,"trainClass"]))

Ahora, se configurarán algunos parámetros. En el ejemplo se usará una validación cruzada (HOCV) para evaluar el rendimiento del clasificador RF Esto significa que se usará un conjunto de datos que se repitan con set de prueba y entrenamiento, para esto es necesario definir el número de instancias para el entrenamiento (con la función pTrain).

Hay que tomar en consideración el echo que RF tiende a tomar algo de tiempo para calibrar un extenso número de observaciones. Así mismo se necesita definir el número de repeticiones en el HOVC (nEvalRounds)

nEvalRounds <- 20 # Número de Repeticiones en el HOVC

pTrain <- 0.5 # Proporción de los datos usados por el algoritmo

Ahora se prepararán algunos objetos que nos permitirán almacenar alguna información del rendimiento de la clasificación

n <- nrow(rstDF)
nClasses <- length(unique(rstDF[,"trainClass"]))


# Preparar Objetos:

confMats <- array(NA, dim = c(nClasses,nClasses,nEvalRounds))

evalMatrix<-matrix(NA, nrow=nEvalRounds, ncol=3, 
                   dimnames=list(paste("R_",1:nEvalRounds,sep=""), 
                                 c("Accuracy","Kappa","PSS")))

pb <- txtProgressBar(1, nEvalRounds, style = 3)

Ahora con todos los sets se calibrará y evaluará el clasificador RF

# Función Evaluación
source(url("https://raw.githubusercontent.com/joaofgoncalves/Evaluation/master/eval.R"))

# Correr el algoritmo
for(i in 1:nEvalRounds){
 
  # Crear un índice que es aleatorio que permite la selección de las filas 
  sampIdx <- sample(1:n, size = round(n*pTrain))

  # Se calibra el RF 
  rf <- randomForest(y = as.factor(rstDF[sampIdx, "trainClass"]),
                     x = rstDF[sampIdx, -1], 
                     ntree = 200)

  # vector que ayudará a predecir el conjunto de prueba
  testSetPred <- predict(rf, newdata = rstDF[-sampIdx,], type = "response")

  # Obtener el vector de las clases observadas de la muestra realizada
  testSetObs <- rstDF[-sampIdx,"trainClass"]
 
  # Se evalúa la clasificación realizada de la muestra
  evalData <- Evaluate(testSetObs, testSetPred)
 
  evalMatrix[i,] <- c(evalData$Metrics["Accuracy",1], 
                    evalData$Metrics["Kappa",1], 
                    evalData$Metrics["PSS",1])
 
  # Se almacena la matriz y se realiza su evaluación
  confMats[,,i] <- evalData$ConfusionMatrix
 
  # Se clasifica toda la imagen completa
  rstPredClassTMP <- predict(rst, model = rf,factors = levels(rstDF[,"trainClass"]))
 
  if(i==1){
    # Se realiza el raster predicho
    rstPredClass <- rstPredClassTMP
   
    # Se Obtiene la precisión para cada clase.
    Precision <- evalData$Metrics["Precision",,drop=FALSE]
    Recall <- evalData$Metrics["Recall",,drop=FALSE]
   
  }else{
    # Se unen todos los raster que se predijeron
    rstPredClass <- stack(rstPredClass, rstPredClassTMP)
   
    # Obtenga precisión que se obtuvo para cada clase
    Precision <- rbind(Precision,evalData$Metrics["Precision",,drop=FALSE])
    Recall <- rbind(Recall,evalData$Metrics["Recall",,drop=FALSE])
  }
 
  setTxtProgressBar(pb,i)
 
}
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |===                                                              |   5%
  |                                                                       
  |=======                                                          |  11%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |==============                                                   |  21%
  |                                                                       
  |=================                                                |  26%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |========================                                         |  37%
  |                                                                       
  |===========================                                      |  42%
  |                                                                       
  |===============================                                  |  47%
  |                                                                       
  |==================================                               |  53%
  |                                                                       
  |======================================                           |  58%
  |                                                                       
  |=========================================                        |  63%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |================================================                 |  74%
  |                                                                       
  |===================================================              |  79%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |==========================================================       |  89%
  |                                                                       
  |==============================================================   |  95%
  |                                                                       
  |=================================================================| 100%
#save.image(file = "./data-raw/P8-session.RData")

Para evaluar la calidad de la clasificación se usarán 3 variables, la exactitud, exhaustividad y la precisión

knitr::kable(evalMatrix, digits = 3)
Accuracy Kappa PSS
R_1 0.789 0.756 0.750
R_2 0.787 0.754 0.748
R_3 0.785 0.752 0.746
R_4 0.790 0.757 0.751
R_5 0.791 0.758 0.752
R_6 0.791 0.758 0.752
R_7 0.788 0.755 0.749
R_8 0.789 0.756 0.750
R_9 0.786 0.752 0.745
R_10 0.795 0.762 0.756
R_11 0.794 0.762 0.757
R_12 0.787 0.753 0.747
R_13 0.786 0.753 0.746
R_14 0.787 0.754 0.747
R_15 0.785 0.750 0.744
R_16 0.782 0.748 0.742
R_17 0.786 0.752 0.747
R_18 0.787 0.754 0.748
R_19 0.786 0.753 0.747
R_20 0.786 0.753 0.748

Ahora se calculasn las estadísticas de las variables anteriores:

round(apply(evalMatrix,2,FUN = function(x,...) c(mean(x,...), sd(x,...))), 3)
##      Accuracy Kappa   PSS
## [1,]    0.788 0.755 0.749
## [2,]    0.003 0.003 0.004

Estas medidas parecen indicar que los resultados son aceptables con muy poca variación entre las rodas de calibración.

Ahora verificamos la exactitud y la exhaustividad:

avgPrecision <- apply(Precision,2,mean)
print(round(avgPrecision, 3))
##     1    10    11    12     2     3     4     5     6     7     8     9 
## 0.965 0.698 0.660 0.489 0.717 0.723 0.907 0.768 0.996 0.653 0.325 0.515
avgRecall <- apply(Recall,2,mean)
print(round(avgRecall, 3))
##     1    10    11    12     2     3     4     5     6     7     8     9 
## 0.978 0.869 0.655 0.269 0.689 0.652 0.964 0.441 1.000 0.621 0.090 0.474
avgF1 <- (2 * avgPrecision * avgRecall) / (avgPrecision+avgRecall)
print(round(avgF1, 3))
##     1    10    11    12     2     3     4     5     6     7     8     9 
## 0.971 0.774 0.657 0.347 0.703 0.686 0.934 0.560 0.998 0.637 0.141 0.494

Ahora se verificará la matriz de confusión para los mejores resultados

# Mejores resultados de Kappa
evalMatrix[which.max(evalMatrix[,"Kappa"]), , drop=FALSE]
##       Accuracy     Kappa       PSS
## R_11 0.7942783 0.7621452 0.7565946
##      Accuracy     Kappa       PSS
## R_4 0.7946858 0.7626567 0.7569426

# esta sentencia muestra los mejores resultados de la variable kappa
cm <- as.data.frame(confMats[,,which.max(evalMatrix[,"Kappa"])])

# Change row/col names
colnames(cm) <- rownames(cm) <- paste(as.factor("c"),levels(as.factor(rstDF[,"trainClass"])),sep="_")

knitr::kable(cm)
c_1 c_2 c_3 c_4 c_5 c_6 c_7 c_8 c_9 c_10 c_11 c_12
c_1 2446 12 0 5 22 4 0 0 0 0 1 0
c_2 16 1699 13 58 17 9 4 1 0 9 11 120
c_3 1 35 550 11 4 8 17 11 3 77 11 115
c_4 5 199 4 119 5 12 16 1 0 0 1 11
c_5 46 29 1 10 369 32 22 0 0 0 1 1
c_6 7 7 4 6 66 337 101 0 0 0 0 1
c_7 0 8 3 11 15 37 2151 5 0 0 0 1
c_8 0 4 45 6 3 6 31 84 0 0 1 4
c_9 0 0 0 0 0 0 0 0 871 0 0 0
c_10 0 26 70 0 0 0 4 0 0 499 3 170
c_11 1 94 20 15 2 4 6 3 0 8 31 92
c_12 3 316 93 15 5 4 12 1 0 162 32 589
Como ob tuvimos una cl asific ación para c ada ro nda, po dremos reuni r toda es info rmación mediante una mejor clasificación: 
rstModalClass <- modal(rstPredClass)

rstModalClassFreq <- modal(rstPredClass, freq=TRUE)

medFreq <- zonal(rstModalClassFreq, rstTrain, fun=median)

Luego usando el modelo de la frecuencia de las 20 clasificaciones, revisamos cuales clases obtuvieron la mayor incertidumbre:

colnames(medFreq) <- c("ClassCode","MedianModalFrequency")

medFreq[order(medFreq[,2],decreasing = TRUE),]
##       ClassCode MedianModalFrequency
##  [1,]         6                   20
##  [2,]        11                   20
##  [3,]        14                   20
##  [4,]        17                   20
##  [5,]         2                   19
##  [6,]         8                   19
##  [7,]        12                   19
##  [8,]        13                   19
##  [9,]         5                   15
## [10,]        16                   14
## [11,]         9                   13
## [12,]         4                   11

Estos valores confirman los valores de la precisión. Finalmente ploteamos nuestros resultados

par(mfrow=c(1,2), cex.main=0.8, cex.axis=0.8)

plot(rstModalClass, main = "Clasificación Supervisada por RF")
plot(rstModalClassFreq, main = "Frecuencia Modal")

Conclusiones:

En la práctica pudimos evaluar la exactitud, la exhaustividad y la precisión de las clasificaciones con resultados aceptables. Sin embargo, se presentaron algunos valores anómalos a la hora de identificar suelo desnudo, o tipos urbanos, quizás debido a la perdida de información a la hora de hacer el resampleo. Si bien en los mapas productos permiten identificar algunas areas, estos podrían mejorarse.