1 Introducción

Este [R Markdown] notebook ilustra como hacer la clasificación de coberturas a partir de imágenes Landsat8 utilizando las herramientas de R y RStudio. Este documento contiene 5 secciones principales: contextualización, métodos, resultados, discusión y conclusiones.

Este trabajo ha sido realizado con base en tres referencias: a- Evaluation of Water Indices for Surface Water Extraction in a Landsat 8 Scene of Nepal (Acharya, 2018). b- Identification of Water Bodies in a Landsat 8 OLI Image Using a J48 Decision Tree (Acharya, 2016). c- Notas y tutoriales de la clase de percepción remota alojados en https://rpubs.com/ials2un

Con el apoyo en la refencia a se aprendió la utilidad de aplicar tresholding sobre índices espectrales, de la referencia b se aprendió de la metodología de árboles de desición y la refencia c sirvió de guia para aplicar estos conceptos en el software R.

1.1 Contexto

Gracias a la teledetección se ha facilitado la tarea de identificar superficies de agua de forma exacta. La práctica más común es la estimación de agua superficial y aquí uno de los principales desafíos es descartar áreas que pueden ser confundidos con agua dentro de la clasificación, tales como nubes, nieve y sombras. Son numerosos los algoritmos que han sido desarrollados para estos fines: técnicas que reconocen patrones estadísticos, clasificación que utiliza datos reales o métodos de clasificación no supervisada. En este contexto, se plantea aplican algunos métodos para identificar espejos de agua en la zona de la Mojana - Colombia.

1.2 Problema

La dinámica de los espejos de agua en la región cenagosa de la Mojana se ve directamente influenciada por la actividad hidrológica de la zona. Lo anterior se debe a la función reguladora que las ciénagas realizan. En temporada húmeda los caudales extremos que fluyen en los ríos son almacenados para después ser liberados gradualmente en la época seca. Está dinámica previene parte de las inundaciones y desbordes, así como también se disminuyen los efectos negativos de la falta de precipitación.

La dinámica de las ciénagas se traduce en una ampliación y reducción periódica de la lámina de agua. Medir estos cambios con instrumentos in situ es una tarea costosa ya que la Mojana se compone de cientos de pequeñas ciénagas interconectadas y cada una de estas puede tener su propio comportamiento y conformación hidráulica. El problema que genera la dificultad de delimitación de este comportamiento es que en la Mojana se localizan pueblos pequeños que en algunas épocas del año pueden verse afectados ante una temporada de lluvias.

1.3 Objetivo

El objetivo de este informe es reportar y analizar los resultados de clasificación digital de la cobertura del suelo en la zona norte de la Mojana, el principal interés de la clasificación es identificar espejos de agua a partir de imágenes de satélite Landsat 8 utilizando Kmeans y el algoritmo Random Forest.

rm(list=ls())
library(terra)
library(sp)
library(raster)
library(dplyr)
library(corrplot)
library(PerformanceAnalytics)
library(rpart)
library(RColorBrewer)
library(sf)
library(ranger)
library(bookdown)
library(knitr)

1.4 Conceptos teóricos

1.4.1 Cobertura del suelo

El término cobertura del suelo hace referencia a la cobertura biofísica que tiene la tierra. Dicho de otro modo, este término se refiere a la descripción del material físico que cubre la tierra. En este orden de ideas, este término no solo pretende describir vegetación sino también elementos antrópicos y superficies como rocas y cuerpos de agua.

Conocer la cobertura del suelo toma importancia en la elaboración de análisis climáticos, biológicos e hidrológicos. Es muy común representar la cobertura en forma de mapas y así poder analizar los cambios en el territorio y ecosistemas.

En el aspecto global, la variable cobertura del suelo se considera como una variable esencial dentro del monitoreo y seguimiento del clima, en el sentido que el clima moldea la cobertura y viceversa.

1.4.2 Índices espectrales

Los índices espectrales han sido utilizados en el campo de los sensores remotos desde los inicios de esta tecnología, estos han surgido como una forma indirecta de estimar variables biofísicas de la vegetación, del agua, de las áreas impermeabilizadas entre otras. Estos índices con base en datos captados por sensores remotos explotan el alto contraste entre las diferentes bandas para caracterizarla y diferenciarla de otros objetos terrestres.

1.4.3 Índices espectrales para la identificación de agua

Los índices utilizados corresponden a los recomendados por Acharya et al. (2018), quienes realizaron un trabajo muy similar al que se propone en este estudio. En la tabla 1.1 se resumen las principales características de los índices utilizados para identificar superficies de agua.

tab_ind_for <- c("Índice espectral", "Ecuación", "Comentario","Referencia",
                 "Índice de vegetación de diferencia normalizada", "NDVI = (NIR-Rojo)/(NIR+Rojo)", "El agua toma valores negativos", "Rouse et al. (1973)",
                 
                 "Índice diferencial normalizado de agua", "NDWI = (Verde-NIR)/(Verde+NIR)", "El agua toma valores positivos","McFeeters (1996)",
           "Índice diferencial normalizado modificado de agua", "MNDWI1 = (Verde-SWIR1)/(Verde+SWIR1)", "El agua toma valores positivos", "Xu (2006)")

tab_ind_for <- matrix(tab_ind_for, nrow = 4, byrow=T)
colnames(tab_ind_for) <- tab_ind_for[1,]
tab_ind_for <- tab_ind_for[-1,]

kable(tab_ind_for, caption = "Índices multibanda utilizados para la extracción de características de agua.")
Table 1.1: Índices multibanda utilizados para la extracción de características de agua.
Índice espectral Ecuación Comentario Referencia
Índice de vegetación de diferencia normalizada NDVI = (NIR-Rojo)/(NIR+Rojo) El agua toma valores negativos Rouse et al. (1973)
Índice diferencial normalizado de agua NDWI = (Verde-NIR)/(Verde+NIR) El agua toma valores positivos McFeeters (1996)
Índice diferencial normalizado modificado de agua MNDWI1 = (Verde-SWIR1)/(Verde+SWIR1) El agua toma valores positivos Xu (2006)

1.4.4 Algoritmos de clasificación

Dentro del aprendizaje de máquina (Machine learning) se estudia el aprendizaje automático a partir de conjuntos de datos con el fin de lograr predicciones precisas. Con el aprendizaje de máquina se puede realizar clasificación de objetos o conjuntos de datos. A continuación, se describe los dos tipos de algoritmos que se aplicarán en este ejercicio.

Clasificación no supervisada: Los datos no presentan etiquetas o clasificación conocida, por tanto se clasifican a partir de relaciones generadas con su estructura interna.

Clasificación supervisada: Se puede definir esta clasificación como un proceso en que pixeles de identidad conocida, ubicados dentro de las áreas de entrenamiento, se utilizan para clasificar pixeles de identidad desconocida. La clasificación supervisada involucra las siguientes etapas: entrenamiento, selección del algoritmo de clasificación y validación.

2 Métodos

2.1 Zona de estudio

La zona de estudio se ubica en la región de la Mojana tal como lo muestra la Figura 2.1, mas exactamente entre las coordenadas 506578.04, 527181.24, 990220.75, 1013472.21 (xmin, xmax, ymin, ymax) (EPSG 32618).

El área de estudio se seleccionó bajo dos criterios. El primero, relacionado con la capacidad instalada del computador disponible para realizar el análisis y clasificación de las imágenes y el segundo, incluyendo al menos dos asentamientos urbanos dentro del área. Como resultado, la Zona de estudio posee un área de 268 km^2 y en ella se localizan seis veredas Santiago Apostol, Sancasaluma, Pinalito, Punta Franco, Punta nueva, Punta Blanco. Los datos de población encontrados no son oficiales, por tanto se calculó a partir de información publicada en noticias y estudios de entes privados que en estas veredas habitan cerca de 4.000 personas.

path <- "G:/Mi unidad/UNIVERSIDAD PC/MAESTRIA/IV SEMESTRE/PERCEPCION REMOTA/PROYECTO1/"
library(leaflet)
area_int <- shapefile("./SHP/SJ_MG_WGS84_32618.shp")
area_wgs84 <- spTransform(area_int, CRS("+proj=longlat +datum=WGS84 +no_defs"))
m <- leaflet() %>% #setView(lng = -122.091, lat = 37.910, zoom = 10)
addTiles() %>%
#addMarkers(lng = -74.84, lat = 9.06, popup = "Área de estudio")  %>%
addPolygons(data=area_wgs84,weight=5,col = 'red', popup = "Área de estudio")
m

Figure 2.1: Localización del área de estudio

2.2 Datos

Para este estudio se utilizaron imágenes Landsat 8 correspondientes al WRS path 099 y WRS row 054, el porcentaje de nubes de la imagen seleccionada es menor al 40%. La imagen fue adquirida el 31-mar-2021 a través del portal Earth Explorer desarrollado por el Servicio Geológico de los Estados Unidos (USGS) (https://earthexplorer.usgs.gov/). La imagen utilizada es de nivel 1 corrección de terreno (LT1) y previamente georefenciadas usando WGS84(EPSG 32618) como origen. En la Tabla 2.1 se presentan las características de las bandas de interés en el presente estudio.

Table 2.1: Especificaciones de las bandas utilizadas - Imagen Landsat8
Banda Longitud de onda (µm) Nombre Resolución (m)
1 0.435–0.451 Coastal Aerosol (CA) 30
2 0.452–0.512 Azul 30
3 0.533–0.59 Verde 30
4 0.636–0.673 Rojo 30
5 0.851–0.879 Infrarrojo cercano (NIR) 30
6 1.566–1.651 Infrarrojo onda corta (SWIR1) 30
7 2.107–2.294 Infrarrojo onda corta (SWIR2) 30

En la 2.2 se presenta la composición en color verdadero (4-3-2) para la imagen completa, el polígono rojo corresponde a la zona de estudio. Cabe anotarse que esta imagen fue seleccionada de forma estratégica para evitar la presencia de nubes dentro del área de interés, de esta forma se facilita el proceso de clasificación de espejos de agua.

path_rast <- "./RASTER/LC08_L1TP_009054_20180220_20200902_02_T1/LC08_L1TP_009054_20180220_20200902_02_T1_B" 
filenames <- paste0(path_rast, 1:7, ".tif")
landsat <- rast(filenames)
names(landsat) <- c('Ultra azul', 'Azul', 'Verde', 'Rojo', 'NIR', 'SWIR1', 'SWIR2')

landsatRGB <- c(landsat$Rojo, landsat$Verde,   landsat$Azul)

plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Composición color verdadero Landsat8")
plot(area_int, add=T, border="red")
Composiciones 4-3-2 Color verdadero

Figure 2.2: Composiciones 4-3-2 Color verdadero

La Figura 2.3 corresponde a la composición 4-3-2 (color verdadero) y 5-4-3 (Color falso) de la imagen a analizar después de recortarla a la ventana del área de estudio. Este proceso se realizó con miras a optimizar la capacidad computacional requerida para los procesos posteriores.

e <- ext(area_int)
print(paste0("Extensión del área de estudio: ", round(e,1)))
## [1] "Extensión del área de estudio: ext(506578, 527181.2, 990220.7, 1013472.2)"
landsatcrop <- crop(landsat, e)
landsatcrop[ landsatcrop == 0 ] <- NA

landsatRGB <- c(landsatcrop$Rojo, landsatcrop$Verde, landsatcrop$Azul)
landsatFCC <- c(landsatcrop$NIR, landsatcrop$Rojo, landsatcrop$Verde)

par(mfrow=c(1, 2))
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Composición 4-3-2")
plot(area_int, add=T, border="yellow", lwd=2)

plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Composición 5-4-3")
plot(area_int, add=T, border="yellow", lwd=2)
Composiciones 4-3-2 Color verdadero y 5-4-3 Color Falso

Figure 2.3: Composiciones 4-3-2 Color verdadero y 5-4-3 Color Falso

La información de coberturas de la tierra para la zona de estudio se tomó de IDEAM (2014) quien produjo una capa de coberturas para los años 2010-2012 (Figura ??). Para fines prácticos se tomó solo el primer nivel de la clasificación, se destaca que dentro del área de interés se encuentran 5 coberturas que se muestran en la 2.2.

#shape de clasificación corine land cover
landcover <- shapefile("./SHP/CORINE_LC_MOJANA_SJ_MG.shp")
landcover@data$CODIGO <- as.numeric(substr(as.character(landcover@data$CODIGO),1,1))
#http://documentacion.ideam.gov.co/openbiblio/bvirtual/021521/LIBROCORINEFINAL.pdf
landcover@data$cod <- ifelse(landcover$CODIGO==1, "Territorios artificializados",
                             ifelse(landcover$CODIGO==2, "Territorios agrícolas",
                                    ifelse(landcover$CODIGO==3,"Bosques y áreas seminaturales",
                                           ifelse(landcover$CODIGO==4,"Áreas húmedas","Superficies de agua"))))
landcover_clasi <- c(1,"Territorios artificializados",2,"Territorios agrícolas",3,
                     "Bosques y áreas seminaturales", 4, "Áreas húmedas", 5, "Superficies de agua")
                        
landcover_clasi <- matrix(landcover_clasi, ncol = 2, byrow = TRUE)

landcover_plot <- landcover
mycolor <- c("black", "darkgreen", "green", "lightblue", "navyblue")
landcover_plot$col = ifelse(landcover$CODIGO==1, mycolor[1], ifelse(landcover$CODIGO==2, mycolor[2],
                                    ifelse(landcover$CODIGO==3,mycolor[3], ifelse(landcover$CODIGO==4,mycolor[4], mycolor[5]))))

plot(landcover_plot, col=landcover_plot$col)
legend("left", fill = mycolor , legend = landcover_clasi[,2], ncol = 1, cex = 0.6 )

Debido a que el objetivo de este trabajo es identificar espejos de agua, se realizó una reclasificación de las coberturas Corine Land Cover ??. Como resultado se tienen dos clasificaciones, agua y no agua. La de agua corresponde a la cobertura 5 de la Figura ?? y no agua son las coberturas restantes de la misma figura. De igual forma se realizó la clasificación para las clases de la Figura ??, los resutlados se encuentran alojados en https://rpubs.com/cagarciae/ieaMOJANA.

#reclasificación

landcover@data$CODIGO <- ifelse(landcover@data$CODIGO==5,2,1)

landcover@data$cod <- ifelse(landcover$CODIGO==1, "No agua",
                                  "Superficies de Agua")
landcover_clasi <- c(1,"No agua",2, "Agua")

                        
landcover_clasi <- matrix(landcover_clasi, ncol = 2, byrow = TRUE)

landcover_plot <- landcover
mycolor <- c("brown", "lightblue")
landcover_plot$col = ifelse(landcover$CODIGO==1, mycolor[1], ifelse(landcover$CODIGO==2, mycolor[2],
                                    ifelse(landcover$CODIGO==3,mycolor[3], ifelse(landcover$CODIGO==4,mycolor[4], mycolor[5]))))

plot(landcover_plot, col=landcover_plot$col)
legend("left", fill = mycolor , legend = landcover_clasi[,2], ncol = 1, cex = 0.6 )

colnames(landcover_clasi) <- c("Código", "Leyenda")
kable(landcover_clasi, caption="Clasificación Corine Land Cover Nivel 1 para la zona de interés")
Table 2.2: Clasificación Corine Land Cover Nivel 1 para la zona de interés
Código Leyenda
1 No agua
2 Agua

2.3 Estadísticas

La aplicación de estadísticas dentro del procesamiento de imágenes satelitales permitirá conocer, cuantificar y mejorar las características de la imagen, bien sea para ajustar la visualización o para orientar cualquier proceso posterior. A continuación, se describen las estadísticas utilizadas en este informe (Salvatierra, 2004).

2.3.0.1 Histograma

Con ayuda de un histograma se puede cuantificar la frecuencia de niveles digitales de los pixeles de una banda. En el histograma se indica el número de pixeles presentes en cada rango y puede representarse de forma numérica y gráfica. En la forma gráfica el rango de los valores discretos que puede tener un pixel se muestra en el eje x, mientras que su frecuencia en el eje y. Con este gráfico se puede evaluar la distribución de los datos, se puede extraer una primera aproximación de las posibles clases temáticas y definir la calidad de la imagen por medio del contraste. Una imagen que genere un histograma muy estrecho corresponde a una imagen con un pobre contraste visual.

2.3.0.2 Medidas de tendencia central

Con este tipo de estadísticas se puede analizar y calcular el nivel digital al cual tienden los datos de cada banda. Las siguientes medidas se pueden aplicar a cada una de las bandas de una imagen multiespectral.

Media o promedio: Este valor expresa el promedio de los niveles digitales de todos los pixeles de la imagen (Ecuación (2.1)).

\[\begin{equation} \bar{\mu} = \frac{\sum_{i=1}^{m} Xi}{m} \tag{2.1} \end{equation}\]

Donde, Xi: Nivel digital del pixel m:Número total de pixeles en la imagen : Media aritmética de la imagen

Varianza: Esta medida calcula la dispersión de los niveles digitales de los pixeles alrededor de la media (Ecuación (2.2)). Si los valores se concentran alrededor de la media la varianza es baja y su histograma estrecho, en el caso contrario la varianza es alta.

\[\begin{equation} \bar{Var} = \frac{\sum_{i=1}^{m} (Xi - \mu)ˆ2}{m-1} \tag{2.2} \end{equation}\]

Donde, Xi: Nivel digital del pixel m: Número total de pixeles en la imagen : Media aritmética de la imagen

Desviación estándar: Se le denomina así a la raíz cuadrada de la varianza (Ecuación (2.3)).

\[\begin{equation} \bar{S} = \sqrt{Var} \tag{2.3} \end{equation}\]

2.3.0.3 Medidas de relación

Por medio de algunos estadísticos es posible conocer la relación entre las bandas de una imagen, tales como la covarianza y el índice de correlación.

Covarianza: Esta estadístico mide la variación conjunta de los niveles digitales de los pixeles de dos bandas alrededor de una media común (Ecuación (2.4)). Esta métrica toma valores positivos cuando ambas variables tienden a aumentar o disminuir a la vez y valores negativos cuando una variable tiene a incrementar mientras la otra disminuye o viceversa.

\[\begin{equation} \bar{COVxy} = \frac{\sum_{i=1}^{m}\sum_{j=1}^{m} (Xj - \mu x)(Yi-\mu y))}{m-1} \tag{2.4} \end{equation}\]

Donde, Xj: Nivel digital del pixel en la banda X Yi: Nivel digital del pixel en la banda Y x: Media aritmética de la banda X y: Media aritmética de la banda Y m:Número total de pixeles en la imagen

Coeficiente de correlación: Corresponde a una medida de relación existente entre los niveles digitales de dos bandas (Ecuación (2.5)). Esta métrica toma valores entre -1 y 1, siendo rxy=1 una relación perfecta entre las dos bandas, rxy=-1 una relación inversa entre las dos bandas y rxy=0 no existe relación entre las dos bandas.

\[\begin{equation} \bar{rxy} = \frac{COVxy}{Sx Sy} \tag{2.5} \end{equation}\]

Donde, COVxy: Covarianza de las bandas X y Y Sx: Desviación estándar de la banda X Sy: Desviación estándar de la banda Y

Matriz de covarianza: Corresponde a una matriz cuadrada que contiene las varianzas y covarianzas asociadas con diferentes variables. Los elementos de la diagonal de la matriz contienen las varianzas de las variables, mientras que los elementos que se encuentran fuera de la diagonal contienen las covarianzas entre todos los pares posibles de variables.

2.3.1 Componentes principales

Este es un método estadístico propuesto por Pearson (1901) y profundizado por Hotelling (1933) que permite simplificar espacios muestrales complejos de muchas dimensiones mientras se conserva la información. Para estudiar las relaciones que presentan las p variables correlacionadas se transforma el conjunto original en un nuevo conjunto de variables interrelacionadas entre sí, es decir que no tengan repetición o redundancia, a este nuevo conjunto se le llama conjunto de componentes principales.

Las nuevas variables se generan a partir de combinaciones lineales de las anteriores y se van construyendo por orden de importancia, este orden se refiere a la variabilidad total que recoge de la muestra, cabe resaltarse que a mayor variabilidad se considera que se agrupa más información.

2.3.2 Algoritmos de clasificación utilizados

Para realizar la clasificación de forma no supervisada se utilizará Kmeans (MacQueen, 1967), el cual tiene como objetivo realizar la agrupación de las características de una imagen minimizando las sumas de las distancias entre cada objeto y el centroide de su grupo (cluster). Algunos ejemplos se pueden ver en (https://estrategiastrading.com/k-means/).

De forma general el algoritmo Kmeans consta de tres fases: la primera, establecimiento del número de centroides, la segunda, asignación de objetos a los centroides y la tercera, actualización de centroides. Existe un proceso iterativo entre la segunda y tercera fase, el cual finaliza cuando los centroides no cambien de posición.

Para el caso de la aplicación de clasificación supervisada, se utilizarán los árboles de clasificación y regresión (CART). La clasificación a través de algoritmos CART son basados en el aprendizaje de máquina. De forma general, este método se basa en particiones sucesivas del conjunto de datos y su mayor ventaja es su carácter descriptivo, el cual permite visualizar las decisiones tomadas por el modelo.

Los árboles de decisión se dividen de la siguiente forma: nodos, ramas y hojas. Las variables de entrada son los nodos, los posibles valores de las variables de entrada son las ramas y los posibles valores de la variable de salida son las hojas. La variable de mayor relevancia en el proceso de clasificación es denominada nodo raíz. En el siguiente enlace se detalla con mayor precisión este método (https://bookdown.org/content/2274/metodos-de-clasificacion.html#arboles-de-clasificacion)).

Un segundo algoritmo utilizado dentro de la clasificación supervisada es Random Forest (Breiman, 2001). Esta es otra técnica de aprendizaje de máquina. Con este método se ejecutan varios algoritmos de árbol de decisión en lugar de uno solo. Para clasificar un nuevo objeto basado en atributos, cada árbol de decisión da una clasificación para después calcular la predicción del bosque como la media de la predicción de cada uno de los árboles. La proporción de árboles que toman una misma respuesta se interpreta como la probabilidad de esta.

2.4 Software

El presente reporte ha sido desarrollado en su totalidad utilizando las herramientas de R (R Core Team, 2020) y RStudio (RStudio Team, 2020)). Cabe acotarse que parte de este trabajo fue realizado con la ayuda de R Notebooks publicados en https://rpubs.com/ials2un/ (Lizarazo, 2021).

3 Resultados

La sección de resultados se compone del cálculo de las estadísticas sobre cada una de las bandas, cálculo de índices espectrales, análisis de componentes principales, clasificación no supervisada y supervisada.

3.1 Estadísticas

En primer lugar, se plotearon los histogramas de cada una de las bandas. Para este caso se convirtieron los niveles digitales a reflectancia TOA mediante la Ecuación (3.1).

\[\begin{equation} \bar{TOA} = FactorAditivo + Banda * FactorMultiplicativo \tag{3.1} \end{equation}\]

Donde, los factores aditivos y multiplicados se tomaron del METADATA ofrecido por el proveedor de la imagen.

La Figura 3.1 evidencia que para las bandas Ultra azul, Azul, Verde y Rojo siguen una distribución normal que se evidencia por la forma en campana de Gauss, según el análisis de estas bandas se identifica que la visualización es muy poco contrastada ( Figura 3.2) y que parecieran identificar una sola cobertura.

Ahora para el caso de las bandas restantes NIR, SWIR1 y SWIR2, los histogramas son bimodales, es decir, distinguen dos zonas la primera muy cercana a TOA=0 y la segunda varía con respecto a la banda NIR TOA=0.25, SWIR1 TOA=0.25 y SWIR2 TOA=0.15, en este caso cada una de las poblaciones puede a una cubierta diferente. En estas bandas los histogramas revelan que las bandas son bien contrastadas, lo cual concuerda con la visualización mostrada en la Figura 3.2.

#calcular reflectancia toa
toa13 <-  -0.1 + landsatcrop*2.0E-5

library(dplyr)
toa13_stack <- toa13 %>% stack()

par(mfrow=c(3, 3))
hist(toa13_stack, 1, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="darkblue", xlab="TOA reflectance" )
hist(toa13_stack, 2, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="lightblue",xlab="TOA reflectance" )
hist(toa13_stack, 3, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="green", xlab="TOA reflectance" )
hist(toa13_stack, 4, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="red", xlab="TOA reflectance")
hist(toa13_stack, 5, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="magenta", xlab="TOA reflectance" )
hist(toa13_stack, 6, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="maroon", xlab="TOA reflectance" )
hist(toa13_stack, 7, maxpixels=1000000, plot=TRUE, ylab="Frequency", col="salmon", xlab="TOA reflectance" )
Histogramas de las bandas recortadas al área de estudio

Figure 3.1: Histogramas de las bandas recortadas al área de estudio 

par(mfrow=c(3, 3))
plot(toa13_stack, c(1:7))
Histogramas de las bandas recortadas al área de estudio

Figure 3.2: Histogramas de las bandas recortadas al área de estudio

Las medidas de tendencia central son presentadas en la Tabla 3.1, se calculó la media, la varianza y la desviación estándar para cada una de las bandas, de igual forma se obtuvieron los valores mínimos, máximos y medios de estos estadísticos. Las desviaciones estándar son consistentemente mayores para las bandas NIR, SWIR1 y SWIR2 identificadas en los histogramas como bien contrastadas. El hecho de que la desviación estándar sea mayor indica que los datos están muy dispersos alrededor de la media, lo anterior se puede relacionar con imágenes bien contrastadas.

toa13_mean <- cellStats(toa13_stack, 'mean')
toa13_var <- cellStats(toa13_stack, 'var')
toa13_sd <- cellStats(toa13_stack, 'sd')


res_metrics <- rbind(toa13_mean, toa13_var, toa13_sd)
res_metrics <- cbind(res_metrics, apply(res_metrics,1,min), apply(res_metrics,1,max), apply(res_metrics,1,mean))
rownames(res_metrics) <- c("Media", "Varianza" , "Desviación estandar")
colnames(res_metrics) <- c(colnames(res_metrics)[1:7], "Mínimo", "Máximo", "Medio")

kable(round(res_metrics,3), caption="Resumen estadísticas de las bandas")
Table 3.1: Resumen estadísticas de las bandas
Ultra.azul Azul Verde Rojo NIR SWIR1 SWIR2 Mínimo Máximo Medio
Media 0.110 0.096 0.091 0.089 0.194 0.181 0.105 0.089 0.194 0.124
Varianza 0.000 0.000 0.000 0.000 0.007 0.009 0.004 0.000 0.009 0.003
Desviación estandar 0.008 0.010 0.014 0.022 0.083 0.097 0.060 0.008 0.097 0.042

La Figura 3.3 corresponde a un gráfico tradicional de correlaciones de Pearson, la diagonal de este gráfico corresponde a los histogramas de cada una de las bandas (discutidos en una sección anterior), bajo la diagonal se encuentran los gráficos de dispersión generados entre bandas y arriba de la diagonal la correlación de Pearson calculada. Debe resaltarse que entre más lineal se vea la línea de tendencia de los gráficos de dispersión mayor es la correlación.

tmp1 <- as.matrix(sample(toa13_stack,500000))
chart.Correlation(tmp1, histogram=TRUE, pch=19)
Correlación entre las bandas

Figure 3.3: Correlación entre las bandas

La Figura 3.4 evidencia que las bandas Rojo, Verde, Azul y Ultra azul tienen una alta correlación lineal entre sí, tomando valores cercanos 1. Lo anterior indica que si estas cuatro bandas aumentan y disminuyen de reflectancia en una misma proporción. Estas cuatro bandas no se correlacionan con las bandas NIR y SWIR1 (comportamiento representado por círculos pequeños). Otro grupo de bandas monotónicamente relacionadas es el compuesto por NIR, SWIR1 y SWIR, lo cual es congruente con lo encontrado en los estadísticos calculados anteriormente. Cabe resaltarse que no existe ningún caso donde la relación sea inversamente proporcional. Los resultados numéricos se pueden ver en la ??.

cmat <- cor(as.matrix(toa13_stack))
# 
ncmat <- round(cmat, 2)
corrplot(ncmat, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)
Correlación entre las bandas - gráfico moderno

Figure 3.4: Correlación entre las bandas - gráfico moderno 

round(cmat, 2)
##            Ultra.azul  Azul Verde  Rojo   NIR SWIR1 SWIR2
## Ultra.azul       1.00  0.99  0.94  0.92 -0.22  0.12  0.30
## Azul             0.99  1.00  0.96  0.95 -0.23  0.12  0.30
## Verde            0.94  0.96  1.00  0.94 -0.19  0.08  0.24
## Rojo             0.92  0.95  0.94  1.00 -0.09  0.27  0.43
## NIR             -0.22 -0.23 -0.19 -0.09  1.00  0.84  0.72
## SWIR1            0.12  0.12  0.08  0.27  0.84  1.00  0.97
## SWIR2            0.30  0.30  0.24  0.43  0.72  0.97  1.00

3.2 Análisis de componenetes principales

La Figura 3.5 corresponde a un Screeplot de los componentes principales identificados por el método PCA, en el eje de las ordenadas se presenta la varianza que es representada por cada componente principal. Visualmente se puede verificar que los dos primeros componentes son los que mayor parte de la varianza representan, por tanto, se trabajará con estos dos.

set.seed(12)
landsat <- rast(filenames)
e <- ext(area_int)
landsatcrop <- crop(landsat, e)
landsatcrop[ landsatcrop == 0 ] <- NA

pca_sample <- spatSample(landsatcrop, 10000)
#plot(pca_sample[,c(6,4)], main = "NIR-Red plot")

#linea de suelo
pca <- prcomp(pca_sample, scale = TRUE )
screeplot(pca, col="aquamarine2", main = "Componentes principales")
Número de componentes principales identificados

Figure 3.5: Número de componentes principales identificados

Se generaron las capas ráster resultantes por los componentes principales PC1 y PC2 (Figura 3.6).

# función que restringe la predicción a los dos componentes principales identificados
pca_predict <- function(model, data, ...) {
  predict(model, data, ...)[,1:2]
}
pci <- predict(landsatcrop, pca, fun=pca_predict)
plot(pci)
Componentes principales identificados

Figure 3.6: Componentes principales identificados

Después del análisis visual es posible definir que el mejor representa las áreas cubiertas por agua. En la Figura 3.7 se plotea la composición 5-6-4 en conjunto a la imagen PC2 generada.

m <- c(-Inf,-3,NA, -3,-2,0, -2,-1,1, -1,0,2, 0,1,3, 1,2,4, 2,4,5,4,Inf,NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)

pcClass <- classify(pci[[2]], rcl)
pcClass <- crop(pcClass , ext(landsatcrop))
par(mfrow=c(1,2))
plotRGB(landsatcrop, r = 5, g = 6, b = 4, stretch = "lin", axes=TRUE, main = "Composición 5-6-4")
plot(pcClass, main="PCA2")
Comparación del componente principal 2 vs falso color

Figure 3.7: Comparación del componente principal 2 vs falso color

3.3 Cálculo de índices espectrales

Para el cálculo de los índices seleccionados se hace necesario definir dos funciones que ayudarán a realizar el cálculo (clickear en code para ver las funciones) 3.8.

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}

vi2 <- function(img, m, k, i) {
  bm <- img[[m]]
  bk <- img[[k]]
  bi <- img[[i]]
  vi2 <- 2.5 * ((bm - bk) / (bm + 6* bk - 7.5* bi + 1))
  return(vi2)
}
# Para Landsat8 NIR = 5, red = 4.
ndvi <- vi(toa13, 5, 4)

# Para Landsat8 verde=3 ,NIR = 5.
ndwi <- vi(toa13, 3,5)

# Para Landsat8 verde=3 ,SWIR1= 6.
mndwi <- vi(toa13, 3,6)

par(mfrow=c(2, 2))
plot(ndvi, col=rev(terrain.colors(10)), main = "Landsat8-NDVI")
# plot(evi, col=rev(terrain.colors(10)), main = "Landsat8- EVI")
plot(seq(1:5), type="n", axes=F, yaxt="n", xaxt= "n" )
plot(ndwi, col=rev(topo.colors(10)), main = "Landsat8-NDWI")
plot(mndwi, col=rev(topo.colors(10)), main = "Landsat8-MNDWI")
Índices espectrales calculados en el área de estudio

Figure 3.8: Índices espectrales calculados en el área de estudio

De forma sencilla utilizando límites sobre los índices espectrales se puede hacer una primera aproximación a la clasificación del agua. En este sentido, los límites se establecieron de forma visual a partir de la Figuras 2.3 y 3.8. Para el caso del NDVI se tomó el valor límite como 0.0 (Figura 3.9).

#El NDVI mayor a 0.1 se toma como  agua
agua_ndwi <- classify(ndwi, cbind(-Inf, 0.0, NA))
#El MNDWI mayor a 0.4 se toma como agua
agua_mndwi <- classify(mndwi, cbind(-Inf, 0.0, NA))

par(mfrow=c(1, 2))
plot(agua_ndwi, main='Agua detectada por NDVI', col="blue", legend=F)
plot(agua_mndwi, main='Agua detectada por MNDWI', col="blue", legend=F)
Identificación de agua a partir de límites aplicados a los índices espectrales

Figure 3.9: Identificación de agua a partir de límites aplicados a los índices espectrales

Se decidió no aplicar métricas de exactitid a estas clasificaciones porque de entrada de forma visual se observa que la clasificación no es lo suficientemente buena, esta decisión se soportó tambien en los resultados del artículo publicado por Acharya et al (2016), quienes encontraron métricas de exactitud kappa de mas o menos 0.5.

3.4 Clasificación no supervisada

Esta clasificación se realizó con la ayuda del algoritmo Kmeans, se establecieron 2 clústeres y se corrió el algoritmo. El resultado de este algoritmo se presenta en la Figura 3.10.

#kmeans algoritmo
nr <- values(toa13)
set.seed(15)
# We want to create 10 clusters, allow 500 iterations, start with 5 random sets using "Lloyd" method
kmncluster <- kmeans(na.omit(nr), centers = 2, iter.max = 500, nstart = 5, algorithm="Lloyd")


# Use the ndvi object to set the cluster values to a new raster
knr <- ndwi

values(knr) <- round(as.numeric(kmncluster$cluster),1)

# par(omi = c(0, 0, 0, 2.3))
plot(knr, main = 'Clasificación no supervisada', col = mycolor, legend=F )
# par(omi = c(0, 0, 0, 0), mai = c(0,0,0,0))
legend("right", legend = seq(1:2),
       fill = mycolor)
Clasificación no supervisada con Kmeans

Figure 3.10: Clasificación no supervisada con Kmeans

De manera comparativa se ploteó el índice NDWI en una misma escala de color que los resultados de Kmeans (3.11).

# Use the ndvi object to set the cluster values to a new raster
knr <- ndwi

values(knr) <- round(as.numeric(kmncluster$cluster),1)


par(mfrow = c(1,2))
plot(ndwi, col = mycolor, main = "Landsat8-NDWI", legend=F)
legend("bottom", legend = seq(1:2), fill = mycolor, horiz = T)
plot(knr, main = 'Clasificación Kmeans', col = rev(mycolor), legend=F )
legend("bottom", legend = seq(1:2),   fill = rev(mycolor), horiz = T)
NDWI vs Clasificación Kmeans

Figure 3.11: NDWI vs Clasificación Kmeans 

print(landcover_clasi)
##      Código Leyenda  
## [1,] "1"    "No agua"
## [2,] "2"    "Agua"

3.5 Clasificación supervisada

Antes de hacer la clasificación por alguno de los dos métodos mencionados se hace necesario definir cuales serán las variables a incluir en la clasificación, como resultado se decidió incluir todas las bandas y los cuatro índices espectrales calculados en una sección anterior.

#agregar indices a la imagen
toa13 <- c(toa13, ndvi, ndwi, mndwi)
names(toa13) <- c("Ultra-blue", "Blue", "Green",      
                  "Red","NIR","SWIR1","SWIR2","NDVI", "NDWI", "MNDWI")

print(paste("Variables a incluir en la clasificación: ", "Ultra-blue,", "Blue,", "Green,",      
                  "Red,","NIR,","SWIR1,","SWIR2,","NDVI,", "NDWI,", "MNDWI"))
## [1] "Variables a incluir en la clasificación:  Ultra-blue, Blue, Green, Red, NIR, SWIR1, SWIR2, NDVI, NDWI, MNDWI"

Adicionalmente, con miras de hacer más comparables los resultados obtenidos con CART y Random Forest se utilizó el mismo valor semilla para la extracción de los muestreos, tanto en entrenamiento como en validación. En consecuencia el área para la validación cuantitativa se presenta en la Figura 3.12.

#obtener datos de entrenamiento
set.seed(612)
# generate stratified random samples from the polygons
# change size parameter as wanted
landcover <- landcover %>% vect()
ptsamp <- spatSample(landcover, size=40, method="random", strata= "CODIGO", chess="")
#extraer las coordenadas de los puntos de muestreo
xy <- as.matrix(geom(ptsamp)[,c('x','y')])

nptsamp <- as(ptsamp, "Spatial")
nref_data <- as(landcover, "Spatial")
ptsamp$class1 <- over(nptsamp, nref_data)$CODIGO

df <- extract(toa13, xy)
sampdata1 <- data.frame(class = as.factor(ptsamp$class1), df)

kable(head(sampdata1), caption= "Ejemplo de puntos de muestreo para entrenamiento")
Table 3.2: Ejemplo de puntos de muestreo para entrenamiento
class Ultra.blue Blue Green Red NIR SWIR1 SWIR2 NDVI NDWI MNDWI
1 0.11422 0.10224 0.09312 0.09584 0.21212 0.25386 0.15908 0.3775815 -0.3898572 -0.4632544
1 0.10506 0.08790 0.07984 0.06880 0.23778 0.21646 0.11634 0.5511775 -0.4972609 -0.4610867
1 0.12236 0.11260 0.11618 0.12180 0.14664 0.15724 0.11302 0.0925346 -0.1158968 -0.1501719
1 0.11578 0.10400 0.10266 0.12248 0.27316 0.35798 0.21542 0.3808513 -0.4536746 -0.5542723
1 0.11232 0.10102 0.09446 0.09558 0.26098 0.24106 0.13848 0.4638770 -0.4684898 -0.4369337
1 0.11228 0.09818 0.09486 0.09958 0.25172 0.30648 0.17916 0.4330771 -0.4525939 -0.5272836 
area_val <- shapefile("./SHP/AREA_VALIDACION.shp")
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Composición color verdadero Landsat8")
plot(area_val, add=T, border="green", lwd=2)
Área para validación cuantitativa

Figure 3.12: Área para validación cuantitativa 

ext_area_val <- ext(area_val)
print(paste("Extensión del área de validación visual:", round(ext_area_val,3)))
## [1] "Extensión del área de validación visual: ext(520808.956, 525266.533, 1003158.342, 1009684.747)"

3.5.1 CART

En la Figura 3.13 se presenta el árbol de clasificación generado después de entrenar el clasificador con los datos de la Tabla (3.2).

#entrenando el clasificador
cartmodel <- rpart(as.factor(class)~., data = sampdata1, method = 'class', minsplit = 4)
#print(cartmodel)

# imprimir el arbol de clasificación entrenado
par(xpd = NA) # otherwise on some devices the text is clipped
plot(cartmodel, uniform=TRUE, main="Árbol de clasificación")
text(cartmodel, cex=0.8 ,digits=3, col="blue")
Árbol de clasificación entrenado

Figure 3.13: Árbol de clasificación entrenado

Los resultados de la clasificación realizada con CART son mostrados en la Figura 3.14, cabe anotarse que hay 5 capas y que cada una representa la probabilidad de una clase de cobertura. El color verde representa los pixeles con mayor probabilidad para esa cobertura, ver por ejemplo la banda X5 - Superficies de agua donde se ve de forma clara la superficie con espejo de agua ya identificada con métodos anteriores.

## Hacer la clasificacion
classified_cart <- predict(toa13, cartmodel, na.rm = TRUE)
# classified_cart@ptr$names <- landcover_clasi[,2]
plot(classified_cart, legend=F)
legend("bottomright", legend = paste0("X",seq(1:2), " ",landcover_clasi[,2]), cex=0.7)
Resultados de la clasificación de coberturas por el método CART

Figure 3.14: Resultados de la clasificación de coberturas por el método CART

Una vez se identifica la probabilidad de clasificación de cada banda se hace un proceso de agrupamiento, aquí se toma la cobertura con mayor probabilidad por pixel y este sería el resultado de la clasificación con CART (3.15).

#crear una banda para almacenar valores NA
classified_cart$X3 <- classified_cart$X2
classified_cart$X3[is.na(classified_cart)] <- 1.0
#classified[[6]]

lulc <- app(classified_cart, fun = which.max)
#mycolor <- topo.colors(5)
#mycolor <- c(mycolor[2:5], mycolor[2])
plot(lulc,  legend=F, axes=TRUE,  maxcell=5000000, col=mycolor)
legend("left", legend = landcover_clasi[,2], cex=0.7, fill=mycolor)
Clasificación de coberturas - Método CART

Figure 3.15: Clasificación de coberturas - Método CART

Una vez realizada la clasificación se debe validar la calidad de esta. Esta validación se hará de dos formas, la primera, una validación cuantitativa que corresponde a una comparación visual. Es necesario resaltar que el objetivo de este trabajo es identificar agua, y esta cobertura es fácilmente identificable ante la visión humana y por consecuencia, es válido hacer este tipo de validación. La segunda corresponde a una validación cuantitativa donde se usarán índices de exactitud para medir el éxito de la clasificación.

Desde una perspectiva cualitativa es notorio que el resultado de este algoritmo ha sido bueno (Figura 3.16).Se identifica pérdida de información en del caño que cruza la ciénaga de este a norte.

#extend area de validacion

region.brick <- crop(toa13, ext_area_val)
region.lc <- crop(lulc, ext_area_val)

par(mfrow=c(1,2))
plotRGB(region.brick, r = 5, g = 6, b = 4, stretch = "lin")
plot(region.lc,  legend=FALSE, axes=FALSE,  col=mycolor)
Composición 5-6-4 vs Clasificación de coberturas método CART

Figure 3.16: Composición 5-6-4 vs Clasificación de coberturas método CART

Con el fin de explorar diferentes métodos para validar la clasificación de coberturas, se optó por aplicar una validación cruzada de 5 veces. Lo que se busca con esta validación es dividir los datos utilizados en 5 grupos. Después uno de los conjuntos se utilizará como validación del modelo y los otros para entrenamiento, así sucesivamente hasta que se hayan cruzado los 5 grupos. Una vez realizado este proceso se toman los resultados de cada grupo y se unen para hacer el cálculo de exactitud.

set.seed(2502)
# number of folds
k <- 10
j <- sample(rep(1:k, each = round(nrow(sampdata1))/k))
table(j)
## j
##  1  2  3  4  5  6  7  8  9 10 
##  8  8  8  8  8  8  8  8  8  8
#testear el modelo 5 veces
x <- list()
for (k in 1:5) {
  train <- sampdata1[j!= k, ]
  test <- sampdata1[j == k, ]
  cart <- rpart(as.factor(class)~., data=train, method = 'class',
                minsplit = 5)
  pclass <- predict(cart, test, na.rm = TRUE)
  # assign class to maximum probablity
  pclass <- apply(pclass, 1, which.max)
  # create a data.frame using the reference and prediction
  x[[k]] <- cbind(test$class, as.integer(pclass))
}

A continuación se presenta la tabla de contingencia de los resultados de la validación cruzada, esta tabla cuenta las ocurrencia de pares observados y predichos posibles.

classdf <- data.frame(classvalue = landcover_clasi[,1],
                      classnames = landcover_clasi[,2])
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
# confusion matrix
conmat1 <- table(y)
# 
print(conmat1)
##         predicted
## observed  1  2
##        1 16  4
##        2  3 17

En la Tabla 3.3 se presenta el resumen de las métricas de exactitud que se calcularon al clasificar con el método CART. Las métricas de exactitud de productor y de usuario se muestran en la Tabla 3.4.

# number of total cases/samples
n1 <- sum(conmat1)

# number of correctly classified cases per class
diag <- diag(conmat1)
# Overall Accuracy
OA1 <- sum(diag) / n1


#kappa index
# observed (true) cases per class
rowsums <- apply(conmat1, 1, sum)
p <- rowsums / n1
# predicted cases per class
colsums <- apply(conmat1, 2, sum)
q <- colsums / n1
expAccuracy <- sum(p*q)
kappa1 <- (OA1 - expAccuracy) / (1 - expAccuracy)

#producer and user accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc_1 <- data.frame(producerAccuracy = PA, userAccuracy = UA)

acc_1 <- c("Número de casos totales/muestras", as.character(n1),
           "Número de clasificados corectamente por clase", as.character(OA1),
           "Índice kappa", as.character(kappa1))
acc_1 <- matrix(acc_1, nrow = 3, byrow = T)
colnames(acc_1) <- c("Métrica","Valor")
acc_1 <- as.data.frame(acc_1)
acc_1[,2] <- round(as.numeric(acc_1[,2]),3)

kable(acc_1,caption = "Métricas de exactitud - clasificación CART")
Table 3.3: Métricas de exactitud - clasificación CART
Métrica Valor
Número de casos totales/muestras 40.000
Número de clasificados corectamente por clase 0.825
Índice kappa 0.650

La tabla 3.4 presetna las métricas de exactitud del productor y del usuario, la cobertura con valor más alto es la última (75.9% - 0.82%), esta corresponde a la 5- Superficies de agua.

kable(outAcc_1, caption="Exactitud de productor y de usuario - Clasificación CART")
Table 3.4: Exactitud de productor y de usuario - Clasificación CART
producerAccuracy userAccuracy
0.8421053 0.80
0.8095238 0.85

#Random Forest

Para aplicar este método se hace necesario definir la fórmula de predicción. En este caso, la formula incluye todos las bandas e índices.

coordenadas <- st_coordinates(as(as(ptsamp, "Spatial"),"sf"))

sampdata1$x <- coordenadas[,1]
sampdata1$y <- coordenadas[,2]
sampdata1$class <- as.factor(sampdata1$class)

#remover filas con valores vacios
sampdata1 <- na.omit(sampdata1)

#construir la formula de prediccion
predictors <- colnames(sampdata1)[2:11]
(fo <- as.formula(paste("class ~", paste(predictors, collapse = "+"))))
## class ~ Ultra.blue + Blue + Green + Red + NIR + SWIR1 + SWIR2 + 
##     NDVI + NDWI + MNDWI
#construir la funcion para evaluar la exactitud
ovAcc <- function(conmat){
  # number of total cases/samples
  n = sum(conmat)
  # number of correctly classified cases per class
  diag = diag(conmat)
  # Overall Accuracy
  OA = sum(diag) / n
  #
  # observed (true) cases per class
  rowsums = apply(conmat, 1, sum)
  p = rowsums / n
  # predicted cases per class
  colsums = apply(conmat, 2, sum)
  q = colsums / n
  expAccuracy = sum(p*q)
  kappa = (OA - expAccuracy) / (1 - expAccuracy)
  #
  # Producer accuracy
  PA <- diag / colsums
  # User accuracy
  UA <- diag / rowsums
  outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
  #print(outAcc)
  #
  global_acc = data.frame(overallAccuracy=OA, overallKappa=kappa)
  #print(global_acc)
  return(list(outAcc,global_acc))
}

La configuración del bosque definido para este ejercicio de clasificación con Random Forest se muestra a continuación.

#modelando landcover randomForest
#library(ranger)
fit <- ranger(fo, data = sampdata1, importance="impurity")
fit
## Ranger result
## 
## Call:
##  ranger(fo, data = sampdata1, importance = "impurity") 
## 
## Type:                             Classification 
## Number of trees:                  500 
## Sample size:                      80 
## Number of independent variables:  10 
## Mtry:                             3 
## Target node size:                 1 
## Variable importance mode:         impurity 
## Splitrule:                        gini 
## OOB prediction error:             12.50 %
fit$variable.importance
## Ultra.blue       Blue      Green        Red        NIR      SWIR1      SWIR2 
##   1.527374   1.629010   1.694291   1.258329   5.819480   3.033085   2.760781 
##       NDVI       NDWI      MNDWI 
##   6.823326   7.227287   7.776636

Se verifica que el entrenamiento hay sido exitoso.

pred <- predict(fit, data =sampdata1 , type = "response")
conmat <- table(pred = pred$predictions, obs = sampdata1$class)
conmat
##     obs
## pred  1  2
##    1 40  0
##    2  0 40
cart_acc <- ovAcc(conmat)
kable(cart_acc[[2]], caption = "Métricas de exactitud global")
Table 3.5: Métricas de exactitud global
overallAccuracy overallKappa
1 1 
cart_acc <- ovAcc(conmat)
kable(cart_acc[[1]], caption = "Métricas de exactitud de productor y usuario")
Table 3.6: Métricas de exactitud de productor y usuario
producerAccuracy userAccuracy
1 1
1 1

Una vez con el modelo entrenado se pueden hacer predicciones, esto significa clasificar todos los pixeles de la imagen y no solo en lo puntos de muestreo.

#clasificaion de las imagenes multiespectrales

classified <- predict(toa13, fit, 
                      fun = function(model, ...) predict(model,   ...)$predictions)

#classified
classified_plot <- classified %>% raster()
#plot(classified)
lc_colors <- c('#fc4e2a','#f7fcb9' , '#99d8c9','#006837',
               '#41ae76', '#d9f0a3',
               '#addd8e', '#2b8cbe') %>%
  colorRampPalette()

# plot(classified_plot,  legend=TRUE, axes=TRUE, col=lc_colors(9))
plot(classified_plot,  legend=F, axes=TRUE, col=mycolor)
legend("left", legend = landcover_clasi[,2], cex=0.7, fill=mycolor)

Tal y como se hizo con el algoritmo CART, se validará la calidad de la predicción. La Figura 3.17 compara la Composición 5-6-4 vs clasificación Random Forest en un área menor a la de la clasificación. El espejo de agua de la ciénaga se diferencia adecuadamente en este punto, y al igual que el método CART se pierde la representación de entidades más finos como es el caso del caño que se aprecia en la imagen de la izquierda.

ext.region <- extent(area_val)
region.brick <- crop(landsatcrop, ext.region)

region.lc <- crop(classified, ext.region)


par(mfrow=c(1,2))
plotRGB(region.brick, r = 5, g = 6, b = 4, stretch = "lin", axes=F)
plot(region.lc, legend=FALSE, axes=FALSE, col=mycolor)
Composición 5-6-4 vs clasificación Random Forest

Figure 3.17: Composición 5-6-4 vs clasificación Random Forest

Para evaluar la calidad de la predicción de forma cualitativa se presentan los índices kappa y el de exactitud global.

set.seed(15)
# Generar aleatoreamente los puntos de validación
ptsamp_val <- spatSample(landcover, size=40, method="random", strata="CODIGO", chess="")
ptsamp_val$reference <-  as.factor(ptsamp$pol.id)


prediction <- terra::extract(classified, ptsamp_val)

ptsamp_val$prediction <- as.factor(prediction$lyr1)

#ptsamp_val$prediction[1:10]
#ptsamp_val$reference[1:10]

#error matrix
(conmat <- table(pred= ptsamp_val$prediction, obs= ptsamp_val$reference))
##     obs
## pred  1  2
##    1 35  5
##    2  5 35
cart_acc <- ovAcc(conmat)
kable(cart_acc[[2]], caption = "Métricas de exactitud global - Validación Random Forest")
Table 3.7: Métricas de exactitud global - Validación Random Forest
overallAccuracy overallKappa
0.875 0.75

Con este algoritmo de clasificación los mejores resultados de las métricas de exactitud productor/usuario se dan en las categorías territorios artificializados y superficies de agua.

cart_acc <- ovAcc(conmat)
kable(cart_acc[[1]], caption = "Métricas de exactitud de productor y usuario - Validación Random Forest")
Table 3.8: Métricas de exactitud de productor y usuario - Validación Random Forest
producerAccuracy userAccuracy
0.875 0.875
0.875 0.875

4 Discusión

Teniendo en cuenta que se busca evaluar la clasificación de agua utilizando como base una imagen Landsat8, se considera interesante comparar los diferentes resultados encontrados, es decir, se incluyen los límites calculados a los índices espectrales, componentes principales, clasificación no supervisada, clasificación supervisada CART y Random Forest. De la Figura 4.1 se detalla que a medida que se complejiza el algoritmo de clasificación se identifican espejos de agua mas amplios. En el caso de los índices espectrales, cuando identifican agua lo hacen con exactitud; no obstante, estos métodos tomaron como no agua muchas zonas del área de estudio, las cuales a simple vista se ven cubiertas por agua.En el caso de Kmeans, la cantidad de pixeles identificados como agua es mayor, lo cual es mejor teniendo en cuenta que la escena representada es en época húmeda y las ciénagas están inundadas, la limitación de este método es que pierde detalle en la representación de cuerpos de agua diferentes a ciénagas, por ejemplo, ríos y caños. Los resultados obtenidos con el método PCA son bastante buenos, de hecho las manchas de agua se asemejan a las producidas por las clasificaciones supervisadas.

#El NDVI mayor a 0.1 se toma como  agua
agua_ndwi <- classify(ndwi, cbind(-Inf, 0.00, NA))
#El MNDWI mayor a 0.4 se toma como agua
agua_mndwi <- classify(mndwi, cbind(-Inf, 0.00, NA))
#PC2 mayor a 0.15 se toma como agua
agua_pc2 <- classify(pcClass, cbind(-Inf, 2.5, NA) )
#nosup menor a 3 se toma como agua
nosup_agua <- classify(knr, cbind(0.2, 1.9, NA) )
#CART 5 se toma como agua
cart_agua <- classify(lulc, cbind(-Inf, 1.2, NA) )
#Random Forest 5 se toma como agua
classified_spat <- as(classified_plot, "SpatRaster")
rf_agua <- classify(classified_spat, cbind(-Inf, 1.2, NA) )



par(mfrow=c(3,2))
plot(agua_ndwi, main='NDVI', col="blue", legend=F)
plot(area_int, add=T, border="red")
plot(agua_mndwi, main='MNDWI', col="blue", legend=F)
plot(area_int, add=T, border="red")
plot(agua_pc2, main="PCA2", col="blue", legend=F)
plot(area_int, add=T, border="red")
plot(nosup_agua, main="KMeans", col="blue", legend=F)
plot(area_int, add=T, border="red")
plot(cart_agua, main="CART", col="blue", legend=F)
plot(area_int, add=T, border="red")
plot(rf_agua, main="Random Forest", col="blue", legend=F)
plot(area_int, add=T, border="red")
Comparación de clasificaciones realizadas por diferentes métodos

Figure 4.1: Comparación de clasificaciones realizadas por diferentes métodos

Los mejores resultados se obtuvieron al utilizar los algoritmos de clasificación supervisada. En la Figura 4.3 se presentan los resultados generado por ambos métodos junto con las composiciones de bandas de color verdadero y 5-6-4. Los resultados por los dos métodos son similares, con la diferencia de que la clasificación con Random Forest respeta mejor los detalles, el río que fluye de sur a norte se ve considerablemente mejor con este algoritmo. Por otro lado, ninguo de los métodos logró representar de forma adecuada los caños.

par(mfrow=c(1,2))
plotRGB(landsatcrop, r = 4, g = 3, b = 2, stretch = "lin", main= "Color verdadero")
plotRGB(landsatcrop, r = 5, g = 6, b = 4, stretch = "lin", main = "Color falso")
Comparación las clasificaciones realizadas usando CART y Random Forest

Figure 4.2: Comparación las clasificaciones realizadas usando CART y Random Forest 

par(mfrow=c(1,2))
plot(lulc, main="CART", col=mycolor, legend=F, axes=F)
plot(classified_plot, main="Random Forest", col=mycolor, legend=F, axes=F ,box=F)
Comparación las clasificaciones realizadas usando CART y Random Forest

Figure 4.3: Comparación las clasificaciones realizadas usando CART y Random Forest

Ahora, para abordar los resultados desde el punto de vista numérico en la Tabla 4.1 se presenta un resumen de las métricas globales de exactitud. Los valores de Random Forest son mejores que los de CART, lo que se explica porque Random Forest corresponde a una versión mejorada de los árboles de desición.

res_res <- c("Métrica de exactitud", "Valor CART", "Valor Random Forest",
             "Exactitud global",as.character(acc_1[2,2]),as.character(cart_acc[[2]][1]),
             "Índice kappa", as.character(acc_1[3,2]),as.character(cart_acc[[2]][2] ))
res_res <- matrix(res_res, nrow=3, byrow=T)
colnames(res_res) <- res_res[1,]
res_res <- res_res[-1,]
res_res[,2] <- round(as.numeric(res_res[,2]),3)
res_res[,3] <- round(as.numeric(res_res[,3]),3)

kable(res_res, caption = "Métricas de exactitud global de las clasificaciones realizadas con CART y Radom Forest")
Table 4.1: Métricas de exactitud global de las clasificaciones realizadas con CART y Radom Forest
Métrica de exactitud Valor CART Valor Random Forest
Exactitud global 0.825 0.875
Índice kappa 0.65 0.75

A continuación, se comparan los resultados de exactitud de productor y de usuario para la clase de interés que es la número 2 (Tabla 4.1). Nuevamente se encontró que la mejor clasificación la hizo el algoritmo Random Forest. Aquí cabe resaltarse que la diferencia de los resutlados obtenidos con ambos algoritmos no es muy amplia.

#[5,]
#

res_2 <- c("Métrica de exactitud", "Valor CART", "Valor Random Forest",
             "Exactitud productor",as.character(outAcc_1[2,1]),as.character(cart_acc[[1]][2,1]),
             "Exactitud usuario", as.character(outAcc_1[2,2]),as.character(cart_acc[[1]][2,2] ))
res_2 <- matrix(res_2, nrow=3, byrow=T)
colnames(res_2) <- res_2[1,]
res_2 <- res_2[-1,]
res_2[,2] <- round(as.numeric(res_2[,2]),3)
res_2[,3] <- round(as.numeric(res_2[,3]),3)

kable(res_2, caption = "Métricas de exactitud de productor y de usuario calculadas para la clase 5 - Superficies de agua")
Table 4.2: Métricas de exactitud de productor y de usuario calculadas para la clase 5 - Superficies de agua
Métrica de exactitud Valor CART Valor Random Forest
Exactitud productor 0.81 0.875
Exactitud usuario 0.85 0.875

5 Conclusiones

Del trabajo realizado para la identificación de espejos de agua en la región de la Mojana se encontró que existe una relación entre la complejidad del algoritmo y los resultados obtenidos, en otras palabras, la clasificación de agua mejoró a medida que los métodos utilizados requerían de un mayor procesamiento. Después del análisis visual se puede concluye que si se debieran ordenar los algoritmos el orden descendente (de mejor a menor representación), el orden sería el siguiente:

  1. Random Forest
  2. CART
  3. PCA
  4. Kmeans
  5. Índices espectrales

Cabe resaltarse que cada metodología tiene potencialidad de ser aplicada frente a diferentes necesidades. Por ejemplo, para aplicar índices espectrales no se requiere ser programador, basta con usar una herramienta de SIG para computar estos índices, el proceso es sencillo y los resultados son adecuados para hacerse una idea general de la cobertura del agua.

Se destaca que evaluar la calidad de la clasificación mediante métricas de exactitud es un proceso importante, en algunas oportunidades a simple vista la imagen parecía estar bien clasificada, pero al momento de analizar estas métricas los valores eran muy bajos. Esto se debe a que el ojo humano puede ser engañado al momento de analizar imágenes con pixeles pequeños y aquí es donde entran las variables cuantitativas.

Al comparar los resultados obtenidos en este ejercicio con los de los artículos de referencia se encontró una gran diferencia en las métricas de exactitud. Las métricas de los artículo que utilizaba árboles de desición rondaban valores de 0.95, la razón al menor desempeño es que la capa de información utilizada para entrenar los modelos es una capa a una resolución muy gruesa. Con seguridad, si se mejora la resolución de la información de entrenamiento las métricas de exactitud podrían aumentar.

Para finalizar, las herramientas que existen para realizar todo el proceso de clasificación se encuentran bien documentadas y son bastante diversas. Tanto las plataformas para descargar las imágenes como las herramientas ofrecidas por R son producto de años de desarrollo, esto facilita que una persona no experta en el lenguaje logre el objetivo de clasificar una imagen y se inicie en el mundo de la teledetección.

6 Referencias

Acharya, Tri D.; Lee, Dong H.; Yang, In T.; Lee, Jae K. 2016. “Identification of Water Bodies in a Landsat 8 OLI Image Using a J48 Decision Tree” Sensors 16, no. 7: 1075. https://doi.org/10.3390/s16071075

Acharya TD, Subedi A, Lee DH. (2018) Evaluation of Water Indices for Surface Water Extraction in a Landsat 8 Scene of Nepal. Sensors (Basel). 2018 Aug 7;18(8):2580. doi: 10.3390/s18082580. PMID: 30087264; PMCID: PMC6111878.

Breiman, L. Random Forests. Machine Learning 45, 5–32 (2001). https://doi.org/10.1023/A:1010933404324

Hotelling H., (1933). Analysis of a Complex of Statistical Variables Into Principal Components, Journal of Educational Psychology, volume 24, pages 417-441 and 498-520.

Huete, A., K. Didan, T. Miura, E.P. Rodriguez, X. Gao and L.G. Ferreira, (2002): Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1): 195–213. DOI: 10.1016/S0034-4257(02)00096-2. (For more information on this paper, please contact the IDMP HelpDesk).

IDEAM (2014). Cobertura de la tierra metogología Corine Land Cover adaptada para Colombia 2010-2012.

Cite as: Lizarazo, I., 2021. A tutorial on pixel-based land cover classification using random forests in R. Available at https://rpubs.com/ials2un/rf_landcover.

MacQueen, J. B. (1967). Some Methods for Classification and Analysis of MultiVariate Observations. In L. M. L. Cam & J. Neyman (eds.), Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability (p./pp. 281-297), : University of California Press.

McFeeters, S.K. (1996) The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int. J. Remote Sens.

Pearson, Karl, (1901). On lines and planes of closest fit to systems of points in space, Philosophical Magazine, Series 6, vol. 2, no. 11, pp. 559-572.

RStudio Team (2020). RStudio: Integrated Development for R. RStudio, PBC, Boston, MA URL http://www.rstudio.com/.

R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. (1973) Monitoring vegetation systems in the Great Plains with ERTS (Earth Resources Technology Satellite). In Proceedings of the Third Earth Resources Technology Satellite Symposium, Greenbelt, ON, Canada, 10–14 December 1973; pp. 309–317

Salvatierra, Cristina. (2004). Fundamentos de procesamiento digital de imágenes. http://ffyl1.uncu.edu.ar/IMG/pdf/Notas_clase_Procesamiento_digital.pdf

Xu, H. (2006) Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens. 2006, 27, 3025–3033.

7 Reproducibilidad

Cite as: García-Echeverri, C. (2021). Identificación de espejos de agua en la zona baja de la Mojana a partir de imágenes Landsat 8. Available at https://rpubs.com/cagarciae/ieaMOJANA_Landsat8.

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] leaflet_2.0.3              knitr_1.30                
##  [3] bookdown_0.21              ranger_0.12.1             
##  [5] sf_0.9-7                   RColorBrewer_1.1-2        
##  [7] rpart_4.1-15               PerformanceAnalytics_2.0.4
##  [9] xts_0.12.1                 zoo_1.8-8                 
## [11] corrplot_0.84              dplyr_1.0.5               
## [13] raster_3.4-5               sp_1.4-5                  
## [15] terra_1.0-10              
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6         highr_0.8          pillar_1.4.7       compiler_4.0.3    
##  [5] class_7.3-17       tools_4.0.3        digest_0.6.27      jsonlite_1.7.2    
##  [9] evaluate_0.14      lifecycle_1.0.0    tibble_3.0.4       lattice_0.20-41   
## [13] pkgconfig_2.0.3    rlang_0.4.10       Matrix_1.2-18      DBI_1.1.0         
## [17] crosstalk_1.1.0.1  rgdal_1.5-19       yaml_2.2.1         xfun_0.20         
## [21] e1071_1.7-4        stringr_1.4.0      htmlwidgets_1.5.3  generics_0.1.0    
## [25] vctrs_0.3.6        classInt_0.4-3     grid_4.0.3         tidyselect_1.1.0  
## [29] glue_1.4.2         R6_2.5.0           rmarkdown_2.7      purrr_0.3.4       
## [33] magrittr_2.0.1     units_0.6-7        codetools_0.2-16   ellipsis_0.3.1    
## [37] htmltools_0.5.0    assertthat_0.2.1   quadprog_1.5-8     KernSmooth_2.23-17
## [41] stringi_1.5.3      crayon_1.3.4