Clasificación no supervisada

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

Las técnicas de clasificación de imagénes se basan en agrupar pixeles que poseen regiones del espectro similares y que pueden asociarse con diferentes clases de coberturas del suelo. La clasificación no supervisada (UC) permite identificar entre las diferentes coberturas basados en la clasificación de las imágenes digitales mediante el uso de bandas dentro de una categoría espectral identificada, en pixeles seleccionados, para luego ser agrupados por su similaridad. Esta técnica de agrupación comprime un conjunto grande de datos en una escala de pixel con un mínimo de trabajo por parte del analista.

Existen varios algoritmos para realizar una clasificación no supervisada, estos algoritmos permiten evaluar el número adecuado de agrupaciones, las cuales permitirán obtener la información requerida. K-means, identifica el centroide y su valor para cada grupo dividiendo el espacio entre las celdas, identificando a cada observación como una única agrupación. El algoritmo CLARA, no toma en cuenta todo el conjunto de datos, y elige aleatoriamente una pequeña muetsra.

El siguiente tutorial realiza es un ejemplo de clasificación no supervisada utilizando 2 algoritmos: K-means y CLARA, El ejemplo esta basado en el script realizado por João Gonçalves, y se puede encontrar en el siguiente link https://www.r-exercises.com/2018/02/28/advanced-techniques-with-raster-data-part-1-unsupervised-classification/

Preparación de datos: Lo primero que realizamos es transformar los niveles digitales a valores de reflectancia utilizando la funcion >landsat8<. Despues de cortar las bandas con la zona de estudio visualizamos en diferentes combinaciones de bandas. Visualmente se puede identificar coberturas como: Páramo, vegetación arbustiva, bosque nativo, infraestructura, pastizales, cultivos, vegetación herbacea.

## cargar librerias
library(raster)
library(rgdal)
library(rgdal)
library(ggplot2)
library(sp)
library(rgeos)
library(landsat8)
library(cluster)
#Ubicacion de carpeta
setwd("D:/maestria/percepcion remota avanzada/taller1/datos")
getwd()

[1] "D:/maestria/percepcion remota avanzada/taller1/datos"

band1 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B1.TIF")
band2 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B2.TIF")
band3 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B3.TIF")
band4 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B4.TIF")
band5 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B5.TIF")
band6 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B6.TIF")
band7 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B7.TIF")
paute <- readOGR("paute.shp")

OGR data source with driver: ESRI Shapefile 
Source: "paute.shp", layer: "paute"
with 1 features
It has 6 fields
Integer64 fields read as strings:  DPA_VALOR 

#transformar a valores de  reflectancia 
#valores de reflectancia      reflconvS(x, Mp, Ap, sunelev)
#SUN_ELEVATION = 61.31854934 
#REFLECTANCE_MULT_BAND_1 = 2.0000E-05
#REFLECTANCE_ADD_BAND_1 = -0.100000
band1.dn<- as(band1, 'SpatialGridDataFrame')
band2.dn<- as(band2, 'SpatialGridDataFrame')
band3.dn<- as(band3, 'SpatialGridDataFrame')
band4.dn<- as(band4, 'SpatialGridDataFrame')
band5.dn<- as(band5, 'SpatialGridDataFrame')
band6.dn<- as(band6, 'SpatialGridDataFrame')
band7.dn<- as(band7, 'SpatialGridDataFrame')
b1.ref <- reflconvS(band1.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b2.ref <- reflconvS(band2.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b3.ref <- reflconvS(band3.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b4.ref <- reflconvS(band4.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b5.ref <- reflconvS(band5.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b6.ref <- reflconvS(band6.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b7.ref <- reflconvS(band7.dn,2.0000E-05 ,-0.100000 , 61.31854934)
b1.ref <- as(b1.ref,'RasterLayer')
b2.ref <- as(b2.ref,'RasterLayer')
b3.ref <- as(b3.ref,'RasterLayer')
b4.ref <- as(b4.ref,'RasterLayer')
b5.ref <- as(b5.ref,'RasterLayer')
b6.ref <- as(b6.ref,'RasterLayer')
b7.ref <- as(b7.ref,'RasterLayer')
#unimos las bandas
#cambiamos los nombres a las bandas
bandas <- brick(b1.ref,b2.ref,b3.ref,b4.ref,b5.ref,b6.ref,b7.ref)
names(bandas) <- paste("Banda",1:7,sep="") 
plot(bandas)

#cortar con la zona de estudio
paute17n<-spTransform(paute, CRS=crs(band1))
bandas <- crop(bandas, extent(paute17n))
bandas <- mask(bandas,paute17n)
plotRGB(bandas, r=4, g=3, b=2, stretch="hist",scale=2000) #Color natural

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

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

Realizamos el agrupamiento de los datos de las bandas con los algoritmos K-mean y KLARA, para los 8 clústeres.

bandas_df <- values(bandas)
# Chequear los NA´s en los datos
idx <- complete.cases(bandas_df)
# Iniciar los datasets ráster que contendrán todas las soluciones de clúster
# 2 grupos /  8 clusteres
bandasKM <- raster(bandas[[1]])
bandasCLARA <- raster(bandas[[1]])
for(nClust in 2:8){
  
  cat("-> Clustering data for nClust =",nClust,"......")
  
  # Realizar clustering K-means
  km <- kmeans(bandas_df[idx,], centers = nClust, iter.max = 50)
  
  # Realizar clara clustering usando distancia Manhattan 
  cla <- clara(bandas_df[idx, ], k = nClust, metric = "manhattan")
  
  # Crear un vector entero temporal para mantener los números del clúster
  kmClust <- vector(mode = "integer", length = ncell(bandas))
  claClust <- vector(mode = "integer", length = ncell(bandas))
  
  # Generar un vector de agrupamiento temporal para K-means
  kmClust[!idx] <- NA
  kmClust[idx] <- km$cluster
  
  # Generar un vector de agrupamiento temporal para CLARA
  claClust[!idx] <- NA
  claClust[idx] <- cla$clustering
  
  # Crear un ráster temporal para mantener la nueva solución de clúster
  # K-means
  tmpbandasKM <- raster(bandas[[1]])
  # CLARA
  tmpbandasCLARA <- raster(bandas[[1]])
  
  # Establecer valores ráster con el vector deL clúster
  # K-means
  values(tmpbandasKM) <- kmClust
  # CLARA
  values(tmpbandasCLARA) <- claClust
  
  # Unir los rásteres temporales en los finales
  if(nClust==2){
    bandasKM    <- tmpbandasKM
    bandasCLARA <- tmpbandasCLARA
  }else{
    bandasKM    <- stack(bandasKM, tmpbandasKM)
    bandasCLARA <- stack(bandasCLARA, tmpbandasCLARA)
  }
  
  cat(" hecho =)\n\n")
}

-> Clustering data for nClust = 2 ...... hecho =)

-> Clustering data for nClust = 3 ...... hecho =)

-> Clustering data for nClust = 4 ...... hecho =)

-> Clustering data for nClust = 5 ...... hecho =)

-> Clustering data for nClust = 6 ...... hecho =)

-> Clustering data for nClust = 7 ...... hecho =)

-> Clustering data for nClust = 8 ...... hecho =)

writeRaster(bandasKM, "kmeans.tif", overwrite=TRUE)
writeRaster(bandasCLARA, "Clara.tif", overwrite=TRUE)

Para evaluar el rendimiento de cada solución de clúster se utilizó el índice silhouette, el cual es un método de interpretar y validar las agrupaciones de datos, proporciona una representación gráfica que representa cómo objeto agrupado se encuentra dentro de su grupo, además el valor de la silueta es una medida que muestra que tan similar es un objeto para su propio clúster en comparación con otros.

Este indice varía de -1 a +1, donde un valor alto indica que el objeto está bien adaptado a su propio clúster y no a los vecinos. Cuando la mayoria de pixeles tienen un valor alto, entonces la configuración de agrupamiento se considera apropiada. Si muchos puntos tienen un valor bajo o negativo, es necesario evaluar el numero de clusters utilizados. Para realizar este analisis el prgrama R cuenta con el paquete >clusterCrit<.

En este ejemplo tomamos las siguientes especificaciones

• Se utilizó un muestreo aleatorio estratificado debido a que su calculo es un proceso bastante lento por el numero de observaciones.

• Se asumió que la muestra es robusta y representativa

• Se usó una sola muestra de celdas

clCritCLARA

$silhouette
[1] 0.2998882

En la tabla podemos observar que los valores de indice silhouette para las dos clasificaciones no supervisadas

knitr::kable(clustPerfSI, digits = 3, align = "c", 
             col.names = c("#clusters","Avg. Silhouette (k-means)","Avg. Silhouette (CLARA)"))

#clusters	Avg. Silhouette (k-means)	Avg. Silhouette (CLARA)
2	0.373	0.321
3	0.379	0.243
4	0.332	0.328
5	0.330	0.261
6	0.310	0.269
7	0.301	0.243
8	0.305	0.300

Comparamos los dos metodos de clasificación

plot(clustPerfSI[,1], clustPerfSI[,2], 
     xlim = c(1,13), ylim = range(clustPerfSI[,2:3]), type = "n", 
     ylab="Avg. Silhouette Index", xlab="# of clusters",
     main="Silhouette index by # of clusters")
# Plot Avg Silhouette values across # of clusters for K-means
lines(clustPerfSI[,1], clustPerfSI[,2], col="red")

# Plot Avg Silhouette values across # of clusters for CLARA
lines(clustPerfSI[,1], clustPerfSI[,3], col="blue")
# Grid lines
abline(v = 1:13, lty=2, col="light grey")

abline(h = seq(0.30,0.44,0.02), lty=2, col="light grey")
legend("topright", legend=c("K-means","CLARA"), col=c("red","blue"), lty=1, lwd=1)

Ploteamos las mejores clasificaciones de cada modelo

plot(bandasKM[[2]])
title("K-mean solution clasifiacion")

plot(bandasCLARA[[2]])
title("CLARA solution clasifiacion")

RESULTADOS

La tabla y la gráfica nos muestra que el algoritmo K-means supera ampliamente el rendimiento de clasificación a CLARA. Sin embargo los dos clasificadores no muestran resultados idoneos de clasificación, obteniendo valores menores a 0.38 de indice de silueta. por lo cual es importante determinar cual es el número de cluster ideal para la clasificación de la zona de estudio. Ademas tambien se deberia evaluar el tamaño de la muestra aleatorias estratificadas por clúster, en este caso utilizamos un valor de 2000, sin embargo se podria aumentar para analizar los resulatdos. Entre las razones del bajo rendimiento de Clara, puede deberse a que este algoritmo trabaja con subconjunto de datos y podria ser menos capaz de encontrar mejores centros de los cluster.

DISCUSIONEs

Al ser este ejercicio un analisis preliminar de clasificación, es indispensable evaluar varios parametros en los clasificadores, entre ellos el numero de clusters y el numero de muestras aleatorias, para obtener mejores resultados.