Entre las aplicaciones más importantes de los datos obtenidos de imágenes satelitales es la generación de mapas de cobertura y uso de suelo, las cuales pueden ser obtenidos a partir de técnicas de clasificacion de imágenes.

Las técnicas de clasificación de imagénes se basan en agrupar pixeles que poseen regiones del espectro similares y que pueden asociarse con diferentes clases de coberturas del suelo. En el caso particular de la clasificación supervisada, utiliza datos de entrenamiento que son variables predictoras medidas en cada unidad de muestreo, y asigna clases previas al muestreo. En siguiente diagrama muestra los pasos a realizar en una clasificación supervisada:

Los algoritmos de clasificación supervisada varían en la forma en que identifican y describen las regiones en el espacio espectral. Algunos pueden manejar clases definidas que se superponen entre sí, y otros generan límites firmes entre clases, entre los algoritmos más utilizados se encuentran : Random Forest, maxima verosimilitud, mínima distancia, paralelepipedos.

Machine learning ofrece la posibilidad de una clasificación eficaz y eficiente de imágenes satelites,entre los algoritmos más conocidos se encuentran árboles de decisiones, random forest, redes neuronales. Algunas de sus ventajas:

En este tutorial utilizaremos el algoritmo random forest para realizar una clasificación supervisada. R studio presenta diferentes paquetes como:

caret::varImp (Calculation of variable importance for regression and classification models)
randomForest::randomForest (Classification and Regression with Random Forest)

Random forest es un método de aprendizaje tanto de clasificación y regresión en donde se basa en la construcción de múltiples árboles de decisión durante la etapa de entrenamiento y como resultado brinda una clasificacion o el predicción promedio de los árboles individuales.

Random forest es uno de los clasificadores supervisados más utilizado en solución de problemas de reconocimiento de patrones, consiste en una combinación de clasificadores de árbol donde cada clasificador se genera usando un vector aleatorio muestreado independientemente del vector de entrada, y cada árbol arroja un voto unitario para la clase más popular para clasificar un vector de entrada. Entre las ventajas del uso de random forest:

El objetivo de este tutorial es realizar una evaluación preliminar del rendimiento de la clasificación de coberturas del suelo en una región del sur del Ecuador (cantón Cuenca), que presenta diversas categorías y un paisaje en ciertas zonas heterogéneo, mediante el uso del método clasificador random forest.

La zona de estudio es el cantón Cuenca, que se encuentra ubicado en la región interandina centro Sur del Ecuador. Se utilizó una imagen Landsat 8 con código LC80100622016325LGN01 capturada el 22 de octubre del 2016, esta fue recortada con la zona de estudio, se tomaron las bandas desde la 1 hasta la 7. Para la identificación de las 10 coberturas: nubes, agua, infraestructura, páramo, plantación forestal, bosque nativo, áreas agropecuarias (pastizales y cultivos), vegetación herbácea, área quemada y suelo descubierto; se utilizó como base el mapa de cobertura y uso de suelo del Ecuador a escala 1:25.000 del año 2012, de acuerdo a esto se generaron los poligonos de entrenamiento.

Código Cobertura
1 Nubes
2 Agua
3 Infraestructura
4 Paramo
5 Plantacion forestal
6 Bosque nativo
7 Pastizal
8 Vegetacion herbacea
9 Area quemada
10 Suelo descubierto

Este tutorial esta basado en el script realizado por João Gonçalves

Preparación de datos: Lo primero que realizamos es transformar los niveles digitales a valores de reflectancia utilizando la funcion >landsat8<. Despues de cortar las bandas con la zona de estudio visualizamos en diferentes combinaciones de bandas.

#librerias utilizadas
library(sp)
library(raster)
library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: D:/Lore/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: D:/Lore/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.4.4
library(maptools)
## Warning: package 'maptools' was built under R version 3.4.4
## Checking rgeos availability: TRUE
library(sp)
library(rgeos)
## rgeos version: 0.3-26, (SVN revision 560)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE
library(landsat8)
## Warning: package 'landsat8' was built under R version 3.4.4
library(randomForest)
## Warning: package 'randomForest' was built under R version 3.4.4
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
#Directorio de trabajo
setwd("D:/maestria/percepcion remota avanzada/taller1/datos")

#datos a utilizar. Imagen Landsat8 bandas y area de estudio
band1 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B1.TIF")
band2 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B2.TIF")
band3 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B3.TIF")
band4 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B4.TIF")
band5 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B5.TIF")
band6 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B6.TIF")
band7 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B7.TIF")
canton.cuenca <- readOGR("canton_cuenca.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "canton_cuenca.shp", layer: "canton_cuenca"
## with 1 features
## It has 6 fields
## Integer64 fields read as strings:  DPA_VALOR
#transformar a valores de  reflectancia 
#valores de reflectancia      reflconvS(x, Mp, Ap, sunelev)
#SUN_ELEVATION = 61.31854934 
#REFLECTANCE_MULT_BAND_1 = 2.0000E-05
#REFLECTANCE_ADD_BAND_1 = -0.100000

band1.dn<- as(band1, 'SpatialGridDataFrame')
band2.dn<- as(band2, 'SpatialGridDataFrame')
band3.dn<- as(band3, 'SpatialGridDataFrame')
band4.dn<- as(band4, 'SpatialGridDataFrame')
band5.dn<- as(band5, 'SpatialGridDataFrame')
band6.dn<- as(band6, 'SpatialGridDataFrame')
band7.dn<- as(band7, 'SpatialGridDataFrame')
b1.ref <- reflconvS(band1.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b2.ref <- reflconvS(band2.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b3.ref <- reflconvS(band3.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b4.ref <- reflconvS(band4.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b5.ref <- reflconvS(band5.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b6.ref <- reflconvS(band6.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b7.ref <- reflconvS(band7.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b1.ref <- as(b1.ref,'RasterLayer')
b2.ref <- as(b2.ref,'RasterLayer')
b3.ref <- as(b3.ref,'RasterLayer')
b4.ref <- as(b4.ref,'RasterLayer')
b5.ref <- as(b5.ref,'RasterLayer')
b6.ref <- as(b6.ref,'RasterLayer')
b7.ref <- as(b7.ref,'RasterLayer')

#unir las bandas
#cambiar los nombres a las bandas
bandas <- stack(b1.ref,b2.ref,b3.ref,b4.ref,b5.ref,b6.ref,b7.ref)
names(bandas) <- paste("Banda",1:7,sep="") 

#cortar con la zona de estudio y colocar en el sistema de proyecciones de la zona de estudio
cuenca.17n<-spTransform(canton.cuenca, CRS=crs(band1))
bandas <- crop(bandas, extent(cuenca.17n))
bandas <- mask(bandas, cuenca.17n)
cuenca.17s <- CRS("+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") 
bandas <- projectRaster(bandas, crs =cuenca.17s )
plot(bandas)

# Ploteamos la zona de estudio en una composicion RGB
plotRGB(bandas, r=4, g=3, b=2, stretch="hist",scale=2000) #Color natural

plotRGB(bandas, r=5, g=4, b=3, stretch="hist",scale=2000) #color infrarrojo

plotRGB(bandas, r=5, g=6, b=2, stretch="hist",scale=2000) #vegetacion saludable

Cargamos las muestras de entrenamiento usadas en la calibración, con estos datos el algoritmo clasificador aprende las firmas espectrales de cada cobertura. En la tabla se observa el numero de pixeles de entrenamiento por cada cobertura, las coberturas que tienes mayor area poseen mayor numero de pixeles; y las con menor area, menor numero de pixeles.

#cargar muestras de entrenamiento
entrenamiento <- readOGR("D:/maestria/percepcion remota avanzada/taller1/datos/entrenamiento.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:/maestria/percepcion remota avanzada/taller1/datos/entrenamiento.shp", layer: "entrenamiento"
## with 169 features
## It has 9 fields
## Integer64 fields read as strings:  OBJECTID Field1 ORIG_FID codigo
#Convertir las muestras de entrenamiento a formato raster
entren <- rasterize(entrenamiento, bandas, field='codigo')

plot(canton.cuenca)
plot(entren, add=TRUE)

# El numero de muestras de entrenamiento por categoria 
tab <- table(values(entren))
print(tab)
## 
##    1    2    3    4    5    6    7    8    9   10 
## 1500  417  863 1300 4450  661 3980 2000  300  670

El clasificador random forest utiliza los datos de entrenamiento en un formato data.frame, es por eso que transformamos a este formato nuestro rasterLayer, obteniendo valores de pixeles que van a ser utilizados para la calibración. Eliminamos los NA y convertimos la clase de referencia en una variable de categoría.

entren.df <- na.omit(values(stack(entren, bandas)))

entren.df[,"layer"] <- as.factor(as.character(entren.df[,"layer"]))

Establecemos algunos parámetros del modelo:

Creamos una matriz en donde se desplegarán los resultados: Precisión, kAPPA y PSS

nEvalRounds <- 8
pTrain <- 0.5

# Objetos que permiten almacenar informacion del procedimiento y validacion de la clasificación
n <- nrow(entren.df)
nClasses <- length(unique(entren.df[,"layer"]))

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)

Para la calibración y evaluación del clasificador Random forest utilizamos el código desarrollado por João Gonçalves que se ecuentra en el siguiente link: https://raw.githubusercontent.com/joaofgoncalves/Evaluation/master/eval.R

## Funcion de Evaluacion

Evaluate = function(actual=NULL, predicted=NULL, cm=NULL){
  
  if(is.null(cm)) {
    naVals = union(which(is.na(actual)), which(is.na(predicted)))
    if(length(naVals) > 0) {
      actual = actual[-naVals]
      predicted = predicted[-naVals]
    }
    f = factor(union(unique(actual), unique(predicted)))
    actual = factor(actual, levels = levels(f))
    predicted = factor(predicted, levels = levels(f))
    cm = as.matrix(table(Actual=actual, Predicted=predicted))
  }
  
  n = sum(cm) # numero de instancias
  nc = nrow(cm) # numero de clases
  diag = diag(cm) # número de instancias correctamente clasificadas por clase
  rowsums = apply(cm, 1, sum) # numero de instancias por calses
  colsums = apply(cm, 2, sum) # numero de predicciones por clases
  p = rowsums / n # distribución de instancias sobre las clases
  q = colsums / n # distribución de instancias sobre las clases predecidas
  
  # Exactitud
  accuracy = sum(diag) / n
  
  # per class prf
  recall = diag / rowsums
  precision = diag / colsums
  f1 = 2 * precision * recall / (precision + recall)
  
  # macro prf
  macroPrecision = mean(precision)
  macroRecall = mean(recall)
  macroF1 = mean(f1)
  
  # 1-vs-all matrix
  oneVsAll = lapply(1 : nc,
                    function(i){
                      v = c(cm[i,i],
                            rowsums[i] - cm[i,i],
                            colsums[i] - cm[i,i],
                            n-rowsums[i] - colsums[i] + cm[i,i]);
                      return(matrix(v, nrow = 2, byrow = T))})
  
  s = matrix(0, nrow=2, ncol=2)
  
  for(i in 1:nc){
    s = s + oneVsAll[[i]]
  }
  
  # promedio exactitud
  avgAccuracy = sum(diag(s))/sum(s)
  
  # micro prf
  microPrf = (diag(s) / apply(s,1, sum))[1];
  
  # majority class
  mcIndex = which(rowsums==max(rowsums))[1] # majority-class index
  mcAccuracy = as.numeric(p[mcIndex]) 
  
  mcRecall = 0*p
  mcRecall[mcIndex] = 1
  
  mcPrecision = 0*p
  mcPrecision[mcIndex] = p[mcIndex]
  
  mcF1 = 0*p
  mcF1[mcIndex] = 2 * mcPrecision[mcIndex] / (mcPrecision[mcIndex] + 1)
  
  # random/exactitud esperada
  expAccuracy = sum(p*q)
  
  # kappa
  kappa = (accuracy - expAccuracy) / (1 - expAccuracy)
  
  # puntaje de Peirce
  a = sum(colsums * rowsums) / n^2
  b = sum(rowsums^2) / n^2
  pss <- (accuracy - a) / (1 - b)
  
  # random guess
  rgAccuracy = 1 / nc
  rgPrecision = p
  rgRecall = 0*p + 1 / nc
  rgF1 = 2 * p / (nc * p + 1)
  
  # random weighted guess
  rwgAccurcy = sum(p^2)
  rwgPrecision = p
  rwgRecall = p
  rwgF1 = p
  
  classNames = names(diag)
  if(is.null(classNames)) classNames = paste("C",(1:nc),sep="")
  
  metrics = rbind(
    Accuracy = accuracy,
    Precision = precision,
    Recall = recall,
    F1 = f1,
    MacroAvgPrecision = macroPrecision,
    MacroAvgRecall = macroRecall,
    MacroAvgF1 = macroF1,
    AvgAccuracy = avgAccuracy,
    MicroAvgPrecision = microPrf,
    MicroAvgRecall = microPrf,
    MicroAvgF1 = microPrf,
    MajorityClassAccuracy = mcAccuracy,
    MajorityClassPrecision = mcPrecision,
    MajorityClassRecall = mcRecall,
    MajorityClassF1 = mcF1,
    Kappa = kappa,
    PSS = pss,
    RandomGuessAccuracy = rgAccuracy,
    RandomGuessPrecision = rgPrecision,
    RandomGuessRecall = rgRecall,
    RandomGuessF1 = rgF1,
    RandomWeightedGuessAccuracy = rwgAccurcy,
    RandomWeightedGuessPrecision = rwgPrecision,
    RandomWeightedGuessRecall = rwgRecall,
    RandomWeightedGuessF1 = rwgF1)
  
  colnames(metrics) = classNames
  
  return(list(ConfusionMatrix = cm, Metrics = metrics))
}

Corremos el clasificador. Utilizaremos 100 numero de árboles.

for(i in 1:nEvalRounds){
  
  # Crear índice aleatorio para la selección de filas en cada ronda
  sampIdx <- sample(1:n, size = round(n*pTrain))
  
  # Calibrar el clasificador random foreste
 rf  <- randomForest(y = as.factor(entren.df[sampIdx, "layer"]), 
                     x = entren.df[sampIdx, -1], 
                     ntree = 100)
  
  # Predecir las clases cone el st de prueba 
  testSetPred <- predict(rf, newdata = entren.df[-sampIdx,], type = "response")
  

  testSetObs <- entren.df[-sampIdx,"layer"]
  
  # Evaluar 
  evalData <- Evaluate(testSetObs, testSetPred)
  
  evalMatrix[i,] <- c(evalData$Metrics["Accuracy",1],
                      evalData$Metrics["Kappa",1],
                      evalData$Metrics["PSS",1])
  
  #Matriz de confusion por ronda evaluadad
  confMats[,,i] <- evalData$ConfusionMatrix
  
  # clasificar toda la imagen 
  rstPredClassTMP <- predict(bandas, model = rf, 
                             factors = levels(entren.df[,"layer"]))
  
  if(i==1){
    # raster predicho
    rstPredClass <- rstPredClassTMP
    
    # obtener la precision para cada clase
    Precision <- evalData$Metrics["Precision",,drop=FALSE]
    Recall <- evalData$Metrics["Recall",,drop=FALSE]
    
  }else{
    # Apilar los rásteres predichos
    rstPredClass <- stack(rstPredClass, rstPredClassTMP)
    
    #Obtenga precisión y recuperación para cada clase
    Precision <- rbind(Precision,evalData$Metrics["Precision",,drop=FALSE])
    Recall <- rbind(Recall,evalData$Metrics["Recall",,drop=FALSE])
  }
  
  setTxtProgressBar(pb,i)
  
}
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |===================                                              |  29%
  |                                                                       
  |============================                                     |  43%
  |                                                                       
  |=====================================                            |  57%
  |                                                                       
  |==============================================                   |  71%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |=================================================================| 100%

La tabla presenta las tres medidas usadas para la evaluación: precisión, indice KAPPA y puntaje de PEIRCE, para cada una de las rondas

knitr::kable(evalMatrix, digits = 3)
Accuracy Kappa PSS
R_1 0.954 0.944 0.940
R_2 0.956 0.947 0.942
R_3 0.954 0.943 0.940
R_4 0.957 0.948 0.944
R_5 0.954 0.944 0.940
R_6 0.954 0.944 0.939
R_7 0.953 0.943 0.938
R_8 0.954 0.945 0.941

Calculamos el promedio de las 20 rondas para cada una de las medidas y la desviación estandar para conocer que tanta separabilidad hay en las rondas de calibración

round(apply(evalMatrix,2,FUN = function(x,...) c(mean(x,...), sd(x,...))), 3)
##      Accuracy Kappa   PSS
## [1,]    0.954 0.945 0.940
## [2,]    0.001 0.002 0.002

Matriz de confusión con la mejor ronda

# Best round for Kappa
evalMatrix[which.max(evalMatrix[,"Kappa"]), , drop=FALSE]
##      Accuracy     Kappa       PSS
## R_4 0.9571305 0.9477819 0.9441173
# Show confusion matrix for the best kappa
cm <- as.data.frame(confMats[,,which.max(evalMatrix[,"Kappa"])])

# Change row/col names
colnames(cm) <- rownames(cm) <- paste("c",levels(as.factor(entren.df[,"layer"])),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_1 738 0 0 0 0 0 0 0 0 0
c_2 0 66 0 3 0 2 49 1 9 20
c_3 0 0 334 0 0 0 0 0 0 0
c_4 0 3 0 168 0 34 3 0 0 4
c_5 0 0 0 0 416 0 0 0 0 0
c_6 0 0 0 14 0 630 0 0 0 3
c_7 0 0 2 0 1 1 2182 0 5 16
c_8 0 0 0 0 0 0 1 274 39 10
c_9 0 0 3 0 1 0 0 19 2008 7
c_10 0 9 0 0 0 4 60 21 2 909

Como obtuvimos un mapa clasificado para cada ronda, podemos obtener toda esa información reuniéndola mediante un voto mayoritario (por ejemplo, calculando la clase de modelo). La función modal del paquete de ráster hace que sea muy fácil calcular esto

rstModalClass <- modal(rstPredClass)

rstModalClassFreq <- modal(rstPredClass, freq=TRUE)

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

Mediante la frecuencia del modelo de las 20 rondas de clasificación,observamos las clases con la mayor incertidumbre:

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

medFreq[order(medFreq[,2],decreasing = TRUE),]
##       ClassCode MedianModalFrequency
##  [1,]         1                    8
##  [2,]         2                    8
##  [3,]         3                    8
##  [4,]         4                    8
##  [5,]         5                    8
##  [6,]         6                    8
##  [7,]         7                    8
##  [8,]         8                    8
##  [9,]        10                    8
## [10,]         9                    6

Observamos la importancia de las variables dentro de la clasificación

varImpPlot(rf)

par(mfrow=c(1,2), cex.main=0.8, cex.axis=0.8)
cols <- c("aliceblue", "navy", "gray47","chartreuse1", "purple4", "darkgreen", "darkseagreen1", "lightcoral", "salmon4", "peru")
plot(rstModalClass,col=cols,main = "RF modal land cover class", legend=TRUE)
plot(rstModalClassFreq, main = "Modal frequency")

RESULTADOS

La matriz de confusión muestra el resultado de los pixeles clasificados de las diferentes coberturas, como relevante se observa que las coberturas de suelo descubierto (cod 10) es la que presenta mayor error de clasificación de las coberturas, confundiendose con coberturas de agua y vegetación herbacea principalmente; la coberturas de nubes presenta 0 error de clasificación.

En general las medidas evaluación precisión, indice KAPPA y puntaje de PEIRCE promedio presentan muy valores muy buenos mayores a 0.94 y con desviaciones estandares en relacion a las rondas muy bajas, lo cual nos dice que no hay variación en las rondas de calibración.

La grafica de disminución en el coeficiente Gini expresa la importancia de estas variables dentro de la clasificación random forest. En este caso la grafica muestra que las bandas 5 y 3 presentan los valores más altos, por lo cual debería ser evaluadas su uso para próximas clasificaciones en la zona de estudio.

DISCUSIONES Y CONCLUSIONES

Este trabajo puede ser visto como una prueba inicial del uso del clasificador de random forest en la zona de estudio, el método en general tuvo buenos resultados en la clasificación, obteniendo valores de porcentaje mayores al 0,94% de las medidas planteadas. Las coberturas que presentan mayor error de clasificación de pixeles, puede ser debido a que el tamaño de la muestra puede llegar a ser más grande que el elemento mapeado en la cobertura, causando una mezcla de clases, por lo cual es recomendable generar nuevos polígonos de entrenamiento que contengan muestras de pixeles más homogéneos. También es recomendable analizar el uso de las bandas 5 y 3 según los resultados obtenidos del coeficiente de impurezas Gini Es importante recalcar que el uso de imágenes satélites y software libre, permite generar información de gran calidad a bajo costo, muy importante en zonas donde los presupuestas de investigación y generación de información base es muy reducido