Clasificaciones no supervisada y supervisada aplicada a una imagen Satelital Landsat-8 de la zona del PNN Sierra Nevada de Santa Marta
Miguel Ángel Ramirez Gutierrez
14/10/2019
Introducción
El empleo de imágenes satelitales para monitoreo de la tierra y las fotografías aéreas han sido empleados históricamente para la delimitación de coberturas de la tierra (Hansen & Loveland, 2012). Con el desarrollo de las ciencias de la computación se han investigado y aplicado métodos para la clasificación de imágenes (Younes Cárdenas, Joyce, & Maier, 2017)con el propósito de delimitar coberturas de la tierra a partir de su respuesta espectral (Goldblatt et al., 2018). Los métodos de clasificación han sido supervisados y no supervisados (Ali, Qazi, & Aslam, 2018), dependiendo sus resultados del propósito y la naturaleza de los datos. También el desarrollo de Machine Learning (Suleiman, Tight, & Quinn, 2019) han motivado la aplicación de estos métodos a imágenes para la detección de fenómenos de interés.
El presente informe muestra de forma detallada el análisis y la obtención de las coberturas presentes en la zona alta del Parque Nacional Natural Sierra Nevada de Santa Marta. La importancia radica porque el monitoreio y conservación de este sistema montañoso es crucial para el aseguramiento del recurso hidrico en la región, además de cuidar la vivienda de las comunidades indigenas arhuacos (o ikas), los wiwas, los kogis y los kankuamos. Esto es congruente con los objetivos de desarrollo sostenible de la agenda 2030 relacionados con hambre cero y la vida en ecosistemas terrestres.
Metodología
1. Zona de Estudio.
El parque Natural Nacional (PNN) Sierra Nevada de Santa Marta se encuentra ubicado en la zona norte de Colombia, específicamente entre los departamentos de Magdalena, Guajira y Cesar. Corresponde a un macizo montañoso que se eleva abruptamente desde las costas del Mar Caribe hasta alcanzar una altura de 5.775 metros en sus picos nevados Bolívar y Colón; ubicados a tan sólo 42 kilómetros del mar. Por sus condiciones topográficas especiales, es lugar de nacimiento de los drenajes más importantes de la región, es decir provee de agua potable a gran parte de la región Caribe.
De acuerdo con la página web oficial del parque (http://www.parquesnacionales.gov.co/portal/es/ecoturismo/region-caribe/parque-nacional-natural-sierra-nevada-de-santa-marta-2/):
Por la variedad de ecosistemas, pisos térmicos junto al mar, su belleza singular y su riqueza cultura constituye un territorio único. En la Sierra Nevada se han declarado figuras de protección como de la Línea Negra (Decreto 1500 de 2018), la Zona de Reserva Forestal (Ley 2ª de 1959), la Reserva de la Biosfera y Patrimonio de la Humanidad por la UNESCO (1979), la Zona de Protección y Desarrollo de los Recursos Naturales Renovables y del Medio Ambiente (Resolución 0504 del 2 de abril de 2018 emitida por el Ministerio de Ambiente y Desarrollo Sostenible). Dentro de este macizo montañoso se encuentra un área de particular belleza paisajística y con gran significado cultural, considerado sitio sagrado para los cuatro pueblos indígenas de la Sierra, se trata del Parque Arqueologico “Ciudad Perdida”>.”
2. Métodos.
Los métodos aplicados se encuentran en la siguiente figura:
Métodos
3. Datos.
La imagen Landsat 8 fue descargada del sitio https://earthexplorer.usgs.gov/ y corresponde a PATH 8 y ROW 53 de fecha de toma del 30/12/2018. El límite del PNN Sierra Nevada de Santa Marta fue obtenido de la publicación realizada por Parques Nacionales. Mientras que las muestras de entrenamiento fueron tomadas por el autor, de acuerdo a criterios de fotointerpretación sobre la imagen empleada.
Resultados y Discusión
1. Descripción, cargue y calibración radiométrica de los datos
Para realizar el análisis fue utilizada una imagen Landsat 8 tomada el 30/12/2018, para visualizar es necesario cargarla. Por esa razón se procede a cargar las librerías raster, rgdal, landsat8 junto a Random Forest que fue utilizada para la clasificación supervisada:
## Loading required package: sp## rgdal: version: 1.4-6, (SVN revision 841)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Program Files/R/R-3.6.1/library/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: C:/Program Files/R/R-3.6.1/library/rgdal/proj
## Linking to sp version: 1.3-1## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.Para el cargue de los datos, se establece el espacio de trabajo y se cargan las bandas espectrales junto al límite del Parque Nacional Natural (PNN) de la Sierra Nevada de Santa Marta:
setwd('R:/Datos')
for(i in 1:7){
nam1<-paste('b',i,sep="")
assign(nam1, raster(paste('LC08_L1TP_008053_20181230_20190130_01_T1_B',i,'.tif',sep="")))
}
names(b1) <- paste("Aerosol costero")
names(b2) <- paste("Azul")
names(b3) <- paste("Verde")
names(b4) <- paste("Rojo")
names(b5) <- paste("NIR")
names(b6) <- paste("SWIR1")
names(b7) <- paste("SWIR2")
zona <- readOGR(dsn = ".", layer = "Zona_Estudio")## OGR data source with driver: ESRI Shapefile
## Source: "R:\DATOS", layer: "Zona_Estudio"
## with 1 features
## It has 6 fieldsLa bandas del visible (R,G,B) corresponde a los raster b2, b3 y b4. El infrarrojo cercano corresponde a b5 mientras que b6 y b7 corresponde a las bandas infrarrojo de banda corta. La banda espectral b1 es la de Aerosol costero. Para el ejercicio fueron utilizadas las siete bandas espectrales.
Es importante conocer los sistemas de referencia de la información cargada, puesto que R no realiza adecuadadamente las operaciones de overlay si no se encuentran en el mismo:
crs(b2)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0crs(zona)## CRS arguments: +proj=longlat +ellps=GRS80 +no_defsObservar que la imagen se encuentra en un sistema de referencia UTM Zona 18 mientras que la zona de estudio se encuentra en un origen geográfico basado en el elipsoide GRS80 (MAGNA).
Otras operaciones que se realizan es obtener el número de pixeles de la imagen con la función ncell, la dimensión de la imagen con la función dim y la resolución espacial con con la función resy por último el número de bandas que posee con nlayers. Además se puede comparar dos imágenes con con la función compareRaster. Esto se realiza con el proposito de conocer las características de la imagen inicial.
ncell(b2)## [1] 58378481dim(b2)## [1] 7721 7561 1res(b2)## [1] 30 30nlayers(b2)## [1] 1compareRaster(b2,b4)## [1] TRUELa calibración radiométrica tiene como objetivo obtener los valores de reflectancia de cada banda espectral de la imagen. Para lo anterior se debió consultar la información disponible en el metadato de la imagen y ubicar los valores de ángulo de elevación solar que corresponde a 48.06474205° y las constantes para obtener la reflectancia TOA que son: 2.0000E-05 y -0.100000.
x1<-2.0000E-05
x2<- -0.100000
ang<- 48.06474205
for(i in 1:7){
nam1<-paste('b',i,sep="")
nam2<-paste('b',i,".dn",sep="")
nam3<-paste('b',i,".ref",sep="")
assign(nam2, as(get(nam1), 'SpatialGridDataFrame'))
assign(nam3, as(reflconvS(get(nam2),x1 ,x2 , ang),'RasterLayer'))
}2. Corte de la Imagen
Se realizo la unión de la bandas espectrales y se elimino de la memoria la información secundaria de las bandas que ya no es necesaria:
remove(b1,b2,b3,b4,b5,b6,b7,b1.dn,b2.dn,b3.dn,b4.dn,b5.dn,b6.dn,b7.dn)
imagenlandsat8<- stack(b1.ref, b2.ref, b3.ref, b4.ref, b5.ref, b6.ref, b7.ref)Puesto que la zona de interés solo corresponde a la región de la sierra nevada, se procede a realizar el corte de la imagen. Pero para ello es necesario reproyectarlo al sistema de referencia de la imagen y así proceder a su corte. También se exporta la imagen:
remove(b1.ref, b2.ref, b3.ref, b4.ref, b5.ref, b6.ref, b7.ref)
zona.utm<- spTransform(zona, crs(imagenlandsat8))
imagencortada <- crop(imagenlandsat8,zona.utm)
## [1] 6748896## [1] 1992 3388 7Se puede visualizar que gran parte del PNN se encuentra dentro de la imagen ademásde que el bosque origianal esta delimitado por tonos verdes mientras los herbazales de tono morado debido a la respuesta que poseen en esas bandas espectrales, lo cual va con el objetivo del informe presentado.
3. Análisis unibanda y multibanda.
Primero se grafica la información de cada banda espectral:
par(mfrow = c(2,2))
plot(imagencortada$Azul, col = gray(0:100 / 100),main = "Azul", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
plot(imagencortada$Verde, col = gray(0:100 / 100),main = "Verde", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
plot(imagencortada$Rojo, col = gray(0:100 / 100),main = "Rojo", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
plot(imagencortada$NIR, col = gray(0:100 / 100),main = "NIR", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
par(mfrow = c(2,2))
plot(imagencortada$SWIR1, col = gray(0:100 / 100),main = "SWIR1", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
plot(imagencortada$SWIR2, col = gray(0:100 / 100),main = "SWIR2", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
plot(imagencortada$Aerosol.costero, col = gray(0:100 / 100),main = "Aerosol Costero", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
axis(1)
axis(2)
A partir de las respuestas espectrales se observa que gran parte de la zona del PNN posee vegetación, puesto que la respuesta es muy alta en el NIR junto a SWIR1 y SWIR2. Mientras que las coberturas de nieves perpetuas junto a la de lagunas se observa con más claridad en las bandas del visible (Azul, Verde y Rojo) junto a la de Aerosolo costero.
Para observar los histogramas junto a la correlación entre las bandas fue realizado:
pairs(imagencortada[[1:7]], main = "Relación entre las bandas de la imagen cortada")
Es claro, además de que va de acuerdo a la revisión de literatura, que las bandas del visible entre ellas poseen una correlación lineal muy alta (por encima de 0.95) además de que presentan la misma característica con la banda de Aerosol costero. Mientras que la correlacion entre las bandas del visible y las del infrarrojo es cercana a 0.8.
4. Indice de Vegetación
Por la presencia de bosques naturales y herbazales, es adecuado aplicarle a la imagen algún índice para la detección de vegetación. Por lo anterior, fue aplicado el NDVI, construyéndose para ello una función:
vi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}
ndvi <- vi(imagencortada, 5, 4)
names(ndvi) <- paste("NDVI")
ndvi1<-mask(ndvi,zona.utm)
plot(zona.utm, main = "Sierra Nevada-NDVI", ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plot(ndvi1, col = rev(terrain.colors(10)), add=TRUE)
grid()
axis(1)
axis(2)
Una forma de detectar la presencia de vegetación, es realizando la clasificación de la imagen utilizando como criterio de decisión que los pixeles que están por encima de 0.4 corresponden a esta cobertura. Esto se muestra a continuación junto al histograma del NDVI obtenido:
hist(ndvi,
main = "Distribución de los valores de NDVI",
xlab = "NDVI",
ylab= "Frecuencia",
col = "green",
xlim = c(-1, 1),
breaks = 100,
xaxt = 'n')
axis(side=1, at = seq(-1,1, 0.05), labels = seq(-1,1, 0.05))
vegetacion <- reclassify(ndvi, c(-Inf, 0.4, NA, 0.4, 1, 1, 1, Inf, NA))
vegetacion<-mask(vegetacion,zona.utm)
plot(zona.utm, main='Vegetación', ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plot(vegetacion, add= TRUE)
grid(equilogs = TRUE)
axis(1)
axis(2)
El criterio anterior provoco que se delimitará solamente la cobertura de bosque original, es decir que el NDVI diferencia sin ningún problema esta cobertura de las demás de vegetación que se desean analizar.
5. Análisis de Componentes Principales (PCA).
Una de las técnicas para bajar la dimensional y eliminar el efecto de multicolinealidad de un conjunto de datos es el empleo del análisis de componentes principales. Debido a que su uso da como resultado un espacio ortogonal denso que explica el espacio original. Para lo anterior fue utilizada la imagen y reducida su dimensional en solo dos bandas. Fue realizado la obtención de los componentes principales por una muestra de 10000 pixeles de la imagen y a partir de esto extrapolarlo a toda la imagen:
set.seed(1)
sr <- sampleRandom(imagencortada, 10000)
(pca <- prcomp(sr, scale = TRUE))## Standard deviations (1, .., p=7):
## [1] 2.55740044 0.52118800 0.40970604 0.11328448 0.08088729 0.02561662
## [7] 0.01321586
##
## Rotation (n x k) = (7 x 7):
## PC1 PC2 PC3 PC4 PC5
## Aerosol.costero 0.3822825 0.23628380 -0.3891244 0.5126647 -0.1559784
## Azul 0.3851508 0.14235799 -0.3753223 0.2047806 -0.0135011
## Verde 0.3885916 0.05591591 -0.2424377 -0.2919060 0.1535355
## Rojo 0.3865063 -0.19153945 -0.1905573 -0.7127234 -0.1893167
## NIR 0.3534427 0.68069677 0.5808208 -0.1200139 0.2103500
## SWIR1 0.3762137 -0.33881101 0.4870410 0.1983000 -0.6514720
## SWIR2 0.3723841 -0.55312678 0.1951039 0.2198405 0.6688883
## PC6 PC7
## Aerosol.costero 0.22585568 -0.555380225
## Azul 0.05400206 0.803425549
## Verde -0.80403794 -0.178422802
## Rojo 0.47434726 -0.093768839
## NIR 0.12521176 0.006993108
## SWIR1 -0.19808697 0.059027628
## SWIR2 0.14026481 -0.043585974screeplot(pca)
pci <- predict(imagencortada, pca, index = 1:2)
pci1<-mask(pci,zona.utm)
pc1 <- reclassify(pci1[[1]], c(-Inf,0,1,0,Inf,NA))
pc2 <- reclassify(pci1[[2]], c(-Inf,0,1,0,Inf,NA))par(mfrow = c(1,2))
plot(zona.utm, main='Primer Componente Principal', ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plot(pc1, legend = FALSE, add= TRUE)
grid(equilogs = TRUE)
axis(1)
axis(2)
plot(zona.utm, main='Segundo Componente Principal', ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plot(pc2, legend = FALSE, add= TRUE)
grid(equilogs = TRUE)
axis(1)
axis(2)
Referente a la descomposición se observa que como se sabe el primer componente posee la mayor parte de la variabilidad siendo de más del 90%. Para analizar los componentes fue reclasificados los obtenidos, para interpretar a que hacen referencia los valores altos. Por lo tanto el componente 2 describe los bosques naturales presentes en la zona de estudio junto a las nieves perpetuas y los cuerpos de agua de la parte alta del PNN.
6. Clasificación no supervisada por el algoritmo k-means
Se realizo y analizo la clasificación no supervisada solamente con las siete bandas en reflectancia TOA. El procedimiento que se realiza con el algoritmo k-means es el siguiente: * Se obtienen los valores, luego inicializar el generador de números aleatorios para luego aplicar el algoritmo. * Se establece el número de clusters, que en este caso son doce puesto que las coberturas de interés son 6. * Se aplica el algoritmo para la generación de los grupos. * Se debe almacenar esa información en un raster temporal, para que sea tenida en cuenta en las siguientes generación.
El procedimiento anterior se encuentra en el siguiente bloque de código:
imagencor_df <- values(imagencortada)
idx <- complete.cases(imagencor_df)
# Iniciar los datasets ráster que contendrán todas las soluciones de clúster
# 2 grupos / 8 clusteres
imagencortadaKM <- raster(imagencortada[[1]])
for(nClust in 2:12){
cat("-> Clustering data for nClust =",nClust,"......")
# Realizar clustering K-means
km <- kmeans(imagencor_df[idx,], centers = nClust, iter.max = 50)
# Crear un vector entero temporal para mantener los números del clúster
kmClust <- vector(mode = "integer", length = ncell(imagencortada))
# Generar un vector de agrupamiento temporal para K-means
kmClust[!idx] <- NA
kmClust[idx] <- km$cluster
# Crear un ráster temporal para mantener la nueva solución de clúster
# K-means
tmpimagencortadaKM <- raster(imagencortada[[1]])
# Establecer valores ráster con el vector deL clúster
# K-means
values(tmpimagencortadaKM) <- kmClust
# Unir los rásteres temporales en los finales
if(nClust==2){
imagencortadaKM <- tmpimagencortadaKM
}else{
imagencortadaKM <- stack(imagencortadaKM, tmpimagencortadaKM)
}
cat("Done \n\n")
}## -> Clustering data for nClust = 2 ......Done
##
## -> Clustering data for nClust = 3 ......Done
##
## -> Clustering data for nClust = 4 ......Done
##
## -> Clustering data for nClust = 5 ......Done
##
## -> Clustering data for nClust = 6 ......Done
##
## -> Clustering data for nClust = 7 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 337444800)## Done
##
## -> Clustering data for nClust = 8 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 337444800)## Done
##
## -> Clustering data for nClust = 9 ......Done
##
## -> Clustering data for nClust = 10 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 337444800)## Done
##
## -> Clustering data for nClust = 11 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 337444800)## Done
##
## -> Clustering data for nClust = 12 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 337444800)## DonePuesto que este ejercicio es exploratorio, dentro de los resultados se identifica la que contenga los doce cluster:
imagencortadaKM<-mask(imagencortadaKM,zona.utm)
plot(imagencortadaKM$layer,, main='Clasificación No Supervisada por K-Means', ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
grid()
plot(zona.utm, main = "Imagen Original", col="red", lty = 2, ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plotRGB(imagencortada1, r = 6, g = 5, b = 4, axes = TRUE, stretch = "lin", add=TRUE)
grid()
axis(1)
axis(2)
Se observa que se puede clasificar hasta más de seis clases de vegetación en la zona de estudio, por lo que si se realizará el procedimiento de asignación de coberturas probablemente se puedan obtener las de interés. Sin embargo existen cluster en que se confunde los cuerpos de agua con nieve, sombras o con vegetación, por lo que se observa que se debería generar más para evitar este tipo de confusiones. Es decir esta clasificación necesita de bastante tiempo de edición.
7. Clasificación supervisada por Random Forest.
Para su aplicación se debe poseer una muestra de entrenamiento, la cual fue realizada por fotointerpretación desde la imagen. Para su uso de debe realizar la transformación de las coordenadas:
entrenamiento <- readOGR("R:/Datos/Entrenamiento.shp")## OGR data source with driver: ESRI Shapefile
## Source: "R:\DATOS\Entrenamiento.shp", layer: "Entrenamiento"
## with 1187 features
## It has 3 fieldscrs(entrenamiento)## CRS arguments:
## +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666
## +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defsentrenamiento.utm<- spTransform(entrenamiento, crs(zona.utm))Los códigos de las coberturas corresponden a:
| Código | Cobertura |
|---|---|
| 1 | Nieve Perpetúa |
| 2 | Boque Original |
| 3 | Laguna |
| 4 | Herbazal |
| 5 | Nube |
| 6 | Vegetación Secundaria |
Se genero una nueva imagen que contiene las bandas espectrales en reflectancia TOA, el indice NDVI junto a los componentes principales generados previamente. Con la imagen generada se convierte a formato raster los puntos de entrenamiento, identificando el atributo que indica la cobertura (este debe ser de tipo entero):
imagenclas<-stack(imagencortada,ndvi,pci)
remove(ndvi,pc1,pc2,pca,pci,sr,vegetacion)
entren.raster <- rasterize(entrenamiento.utm, imagenclas, field='ID_COBERT')
plot(zona.utm, main='Puntos de Entrenamiento', ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plot(entrenamiento.utm, add=TRUE)
grid(equilogs = TRUE)
axis(1)
axis(2)
La recomendación es que cada cobertura posea al menos una muestra de 200 datos, por lo que los datos de entrenamiento para el ejercicio presentado tiene estas características:
tab <- table(values(entren.raster))
print(tab)##
## 1 2 3 4 5 6
## 124 232 267 235 123 171El siguiente paso es clasificar los datos de entrenamiento como una variable categorica:
entren.df <- na.omit(values(stack(entren.raster, imagenclas)))
entren.df[,"layer"] <- as.factor(as.character(entren.df[,"layer"]))Para evaluar la calidad de la clasificación es necesario establecer los siguientes criterios:
- Para evaluar el rendimiento del clasificador RF (validación cruzada de holdout (HOCV)), fue utilizado un conjunto de datos de entrenamiento y otro de prueba. Para el ejercicio presentado el 50% correspondio a los datos de entrenamineto.
- Se debe establecer el número de repeticiones para la validacion cruzada, que para el ejercicio es 8.
- Los resultados de la calidad de la clasificación fueron almacenados en una matriz de evaluación conteniendo la información de los indicadores de Exactitud (Accuracy),exhaustividad (kappa) y la precisión de la clasificación (PSS).
Por eso se tomo el método propuesto por Gonçalves (https://www.r-exercises.com/2018/03/07/advanced-techniques-with-raster-data-part-2-supervised-classification/):
nEvalRounds <- 8
pTrain <- 0.5
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 de los resultados del algoritmo Random Forest se utiliza la propuesta de Gonçalves:
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))
}Para obtener la clasificación supervisada se decide aplciar el algoritmo Random Forest con 100 á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(imagenclas, 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 indica los valores que se obtuvieron para los indicadores de evaluación para cada una de las iteraciones realzidas:
knitr::kable(evalMatrix, digits = 3)| Accuracy | Kappa | PSS | |
|---|---|---|---|
| R_1 | 0.957 | 0.947 | 0.946 |
| R_2 | 0.964 | 0.956 | 0.955 |
| R_3 | 0.950 | 0.939 | 0.938 |
| R_4 | 0.955 | 0.945 | 0.944 |
| R_5 | 0.951 | 0.941 | 0.941 |
| R_6 | 0.967 | 0.960 | 0.959 |
| R_7 | 0.955 | 0.945 | 0.944 |
| R_8 | 0.953 | 0.943 | 0.943 |
Se observa que las ocho realizadas fueron adecuadas puesto que los tres indicadores estuvieron por encima de 0.75. De forma adicional fue calculado el promedio de las iteraciones para cada una de las medidas y la desviación estandar para conocer que tanta separabilidad hay entre ellas:
round(apply(evalMatrix,2,FUN = function(x,...) c(mean(x,...), sd(x,...))), 3)## Accuracy Kappa PSS
## [1,] 0.956 0.947 0.946
## [2,] 0.006 0.007 0.007Lo anterior quiere decir que las iteraciones en conjunto poseen una buena cohesión, puesto que su desviación esta por debajo del 1%. Por otro lado fue escogida la Matriz de confusión de la iteración que tuviera mejor indicador de exhaustividad:
evalMatrix[which.max(evalMatrix[,"Kappa"]), , drop=FALSE]## Accuracy Kappa PSS
## R_6 0.9670139 0.9597779 0.9591751Se visualiza cuál fue la escogida en la salida anterior. A continuación se muestra la matriz de confusión de la iteración seleccionada:
cm <- as.data.frame(confMats[,,which.max(evalMatrix[,"Kappa"])])
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_1 | 66 | 0 | 0 | 0 | 0 | 0 |
| c_2 | 0 | 118 | 1 | 0 | 0 | 2 |
| c_3 | 0 | 0 | 125 | 0 | 0 | 0 |
| c_4 | 0 | 1 | 0 | 114 | 0 | 4 |
| c_5 | 0 | 0 | 0 | 0 | 61 | 1 |
| c_6 | 0 | 2 | 0 | 8 | 0 | 73 |
Es claro que sopresivamente las coberturas que no presentan problemas de omisión o confusión son las que corresponden a la nieve perpetúa (1) y nubes (2), podría pensarse que estás dos por su alta reflectancia se pudieron haber confundido. Otra cobertura que no presenta estas condiciones es la de laguna (3), sin embargo se ve en una baja proporción se podría confundir con la cobertura de bosque original (2). Mientras que las condiciones de confusión y omisión las presenta las tres coberturas de vegetación, lo cual es un comportamiento característico en las zonas de límites de coberturas. La cobertura que presenta más estos problemas es la de vegetación secundaría, pero no es preocupante este aspecto porque en ella se clasifican la vegetación que no es bosque natural (1) o herbazal (4).
A continuación se realizan las demás operaciones para obtener el raster de clasificación:
rstModalClass <- modal(rstPredClass)
rstModalClassFreq <- modal(rstPredClass, freq=TRUE)
medFreq <- zonal(rstModalClassFreq, entren.raster, fun=median)
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 8El aporte de las bandas empleadas en la clasificación es la siguiente:
varImpPlot(rf)
Se observa que el indice NDVI fue el que más aporto para realizar la separabilidad de las coberturas al igual que el NIR. Sorpresivamente las bandas del Azul también la aportaron. Es de señalar que el segundo componente principal brindo más información para diferenciar las coberturas que el primero que es que posee la mayor parte de la variabilidad. Mientrás que las dos bandas que no proporcionaron tanta información para la clasificación fueron las dos del SWIR.
cols <- c("azure1","chartreuse4", "darkblue", "darkseagreen3", "gray", "greenyellow")
rstModalClass<-mask(rstModalClass,zona.utm)
plot(rstModalClass,col=cols,main = "Clasificación supervisada por Random Forest", ylab ="Norte", xlab= "Este")
grid()
plot(zona.utm, main = "Imagen Original", col="red", lty = 2, ylab ="Norte", xlab= "Este",xlim=c(ext@xmin, ext@xmax), ylim=c(ext@ymin, ext@ymax))
plotRGB(imagencortada1, r = 6, g = 5, b = 4, axes = TRUE, stretch = "lin", add=TRUE)
grid()
axis(1)
axis(2)
En terminos generales la clasificación fue adecuada, no obstante se encuentra problemas en la clasificación en las lagunas, debido a que las confunde con las sombras que presenta la imagen. Las demás coberturas no presentan estos problemas por lo que se observa que no necesitaría de tanto de tiempo de edición la clasificación obtenida.
Conclusiones.
La aplicación de técnicas de procesamiento digital de imágenes en software libre es una herramienta poderosa para el monitoreo y análisis de un gran volumen de información (tanto imagenes como información geoespacial complementaria) con el propósito de monitorear y tomar decisiones acertadas sobre el territorio. El software libre también permite la replicación de trabajos de otros autores, que permiten que el análisis que se realice sea mucho más robusto y este de acuerdo a la necesidad planteada con una alta probabilidad de que el trabajo realizado sea replicable tiempo después.
En relación a los métodos de clasificación, sin lugar a duda el que mejor resultado dio fue el supervisada aplicando Random Forest. Puesto que sus indicadores de calidad fueron muy altos y además visualmente fueron delimitadas las coberturas bajo lo estándares necesarios para el ejercicio realizado. No obstante, para la edición de este producto es necesario realizar la labor en otro software, puesto que en R esta labor es compleja realizarla. De forma novedosa se observo que las banda del Aerosol Costero fue útil para la delimitación de las coberturas, por lo que tal vez sea aplicable su uso en zonas de las mismas condiciones del PNN Sierra Nevada de Santa Marta
Referente al rendimiento, este tipo de software demanda de buena memoria RAM, más si se realizan labores complejas como la de generar la clasificación en un raster de forma irregular (la correspondiente a la forma del PNN). Por lo anterior se decidio realizar el corte de la imagen con la maxima extensión del PNN, para luego mostrar los resultados aplicando la mascara correspondiente.
Como recomendación para otro ejercicio, sería adecuado utilizar otro tipo de índices espectrales de vegetación junto a otros métodos de clasificación no supervisada para realizar la evaluación de estos métodos para la obtención de las coberturas de interés.
Referencias.
Ali, M. Z., Qazi, W., & Aslam, N. (2018). A comparative study of ALOS-2 PALSAR and landsat-8 imagery for land cover classification using maximum likelihood classifier. Egyptian Journal of Remote Sensing and Space Science, 21, S29–S35. https://doi.org/10.1016/j.ejrs.2018.03.003
Goldblatt, R., Stuhlmacher, M. F., Tellman, B., Clinton, N., Hanson, G., Georgescu, M., … Balling, R. C. (2018). Using Landsat and nighttime lights for supervised pixel-based image classification of urban land cover. Remote Sensing of Environment, 205(November 2017), 253–275. https://doi.org/10.1016/j.rse.2017.11.026
Hansen, M. C., & Loveland, T. R. (2012). A review of large area monitoring of land cover change using Landsat data. Remote Sensing of Environment, 122, 66–74. https://doi.org/10.1016/j.rse.2011.08.024
Suleiman, A., Tight, M. R., & Quinn, A. D. (2019). Applying machine learning methods in managing urban concentrations of traffic-related particulate matter (PM10 and PM2.5). Atmospheric Pollution Research, 10(1), 134–144. https://doi.org/10.1016/j.apr.2018.07.001
Younes Cárdenas, N., Joyce, K. E., & Maier, S. W. (2017). Monitoring mangrove forests: Are we taking full advantage of technology? International Journal of Applied Earth Observation and Geoinformation, 63(April), 1–14. https://doi.org/10.1016/j.jag.2017.07.004