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:
Capacidad de manejar datos de alta dimensionalidad.
Mapear clases con características muy complejas.
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:
Se ejecuta de manera eficiente cuando se usa bases de datos grandes
Puede manejar miles de variables de entrada sin variables supresión
Da estimaciones de qué variables son importantes en la clasificación
Genera una estimación interna imparcial de la generalización error (oob error)
Calcula proximidades entre pares de casos que se pueden usar en la localización de valores atípicos
Es relativamente robusto a los valores atípicos y el ruido
Es computacionalmente más ligero que otros métodos de ensamblaje de árboles
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:
Para evaluar el rendimiento del clasificador RF: validación cruzada de holdout (HOCV), utiliza un conjunto de datos de entrenamiento y otro de prueba
Porcion de datos para entrenamiento (pTrain) = 0.5, es decir la mitad
Numero de repeticiones para la validacion cruzada holdout = 8
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")
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.
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