title: “CLASIFICACION NO SUPERVISADA” output: html_notebook — Natalia Martínez Z Agosto de 2018
ota: Para la elaboración de este cuaderno se utilizó el tutorial Advanced Techniques With Raster DAta: Part1 - Unsupervised Classification realizado por Joao Goncalves publicado el 28 de Febrero de 2018.
Introduccion
La clasificación de imagenes es una secuencia de pasos en la que se busca extraer informacion temática de una imágen raster y el resultado sirve para crear mapas temáticos. La clasificacion de imagenenes surge de la necesidad de saber que sucede con los pixeles planteado por Blaschke y Strobl en 2001, y se basaron en el primer algoritmo de extraccion y clasificacion de objetos homogeneos desarrollado por Kettig y Landgrebe en 1976. Sin embrgo las imagenes de alta resolución requieren una clasificacion mas exacta que la que se puede obtener por pixeles, es por esto que los métodos basados en objetos permiten hacer un mejor uso de la informacion espacial tal como lo demuestra Salehi et al 2012 quienes realizaron clasificaciones basadas en pixeles y en objetos, encontrando mejores resultados en clasificacion basada en objetos por la utilizacion de mayor cantidad de nformacion como tamaño, forma, ubicación, condiciones de contorno, relaciones topologicas entre otros aspectos que se tienen en cuenta (Blaschke, et al 2015).
Existen dos tipos de clasificación, la supervisada en la cual se utilizan muestras de entrenamiento que definen el numero de clases a identificar, y la clasificacion no supervisada en la que se obtienen clases espectrales donde el algoritmo utiliza sus propios dispositivos para encontrar y presentar la estructura de los datos (Brownlee, 2016).
Este cuaderno tiene como objetivo realizar un ejercicio de clasificación no supervisada utilizando dos algoritmos diferentes Kmeans y clara. Para tal fin se utilizó una imagen Landsat 8 capturada el 14 de Agosto de 2017, de las cuales se tomaron las bandas 1, 2, 3, 4, 5, 6 y 7. Se realizó un corte de la imagen original tomando un area de 4.599,13Km2, la cual inclyue los municipios de Acacias, Castilla la Nueva, San Carlos de Guaroa, Restrepo y villavicencio en el departamento del Meta-Colombia.
Clasificación no supervisada
(Humboldt State University, 2015)
La clasificación no supervisada (CNS) también conocida como “clustering”, utiliza las propiedades y momentos de la distribución estadística de los pixeles en un espacio característico de bandas espectrales para hacer una diferenciación entre grupos relativamente similares.
La CNS es una forma eficaz para las imagenes tomadas por sensores remotos en un espacio con caracteristicas espectrales, para extraer información util sobre la cobertura terrestre. Esta clasificación utiliza clusters conocidos como tecnicas de reducción de datos donde se comprime la informacion diversa de los pixeles en grupos que tienen valores similares.
A diferencia de la clasificación supervisada la CNS no necesita que el usuario ingrese datos de entrenamiento. La CNS requiere unicamente definir el número de grupos (cluster) y las bandas que se van a utiliza, sin necesidad de otras entradas.
Los algoritmos intentan proporcionar la mejor solución a los valores de los grupos de pixeles, donde las distancias dentro de cada grupos se minimiza y la separacionentre grupos se maximiza.
Algoritmo de clasificación Kmeans
Es un algoritmo muy conocido y utilizado por su eficiencia y robustez, su nombre hace referencia al numero K de clases o grupos a buscar
Es una tecnica no jerarquica de particion donde:
- Selecciona K objetos al azar del conljunto total y se asignan como patron o centroides de las K clases que se van a buscar (previamente definidas).
- Calcula todas las distancias euclideanas de todos los objetos restantes a los K centroides, y se asigna la pertenencia a cada objeto al cluster que tenga mas cercano.
- Se recalcula el centroide de cada cluster como la medida de todos los objetos que lo compnen, buscando minimizar el valor de una funcion de coste, que es la sumatoria de todas las distancias euclídeas de los objetos de cada clase al centroide de su resspectiva clase.
- El paso 2 y 3 se repiten sucesivamente hasta que los centros de todos los grups permanezcan constantes o hasta que se cumpla alguna otra condición de parada.
Algoritmo de clasificación CLARA (Clustering Large Applications)
Este algoritmo divide la base de datos original en muestras de tamaño “s”, aplicando el algoritmo PAM (extension del algoritmo K-means, donde cada cluster esta representado por un medioide y no un centroide) sobre cada una de ellas, seleccionando la mejor clasificación de las resultantes. Este algoritmo está indicado para bases de datos con gran cantidad de objetos, y su principal característica es que minimiza la carga computacional en detrimiento de una agrupacion óptima y precisa.
Librerias que se utilizaron para realizar el proceso
library(raster)## Loading required package: splibrary(rgdal)## rgdal: version: 1.3-3, (SVN revision 759)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.3-1library(sp)
library(ggplot2)
library(rgeos)## rgeos version: 0.3-28, (SVN revision 572)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.3-1
## Polygon checking: TRUElibrary(landsat8)
library(cluster)
library(clusterCrit)Seleccion del directorio de trabajo
setwd("/home/maria.orellana/Natalia_Landsat_8")Bandas a utilizar
banda1 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B1.TIF")
banda2 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B2.TIF")
banda3 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B3.TIF")
banda4 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B4.TIF")
banda5 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B5.TIF")
banda6 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B6.TIF")
banda7 <- raster("/home/maria.orellana/Natalia_Landsat_8/LC08_L1TP_007057_20170814_20170825_01_T1_B7.TIF")Correccioón radiométrica
Radiancia
Para realizar la corrección radiométrica de la imágen se utilizó el paquete landsat8.
banda.dn <- as(banda1, ‘SpatialGridDataFrame’) band1.rad <- radconv(banda1.dn, multipbanda, aditivobanda)
banda1.dn <- as(banda1, 'SpatialGridDataFrame')
banda1.rad <- radconv(banda1.dn, 1.2238E-02, -61.19172)
banda2.dn <- as(banda2, 'SpatialGridDataFrame')
banda2.rad <- radconv(banda2.dn, 1.2532E-02, -62.66107)
banda3.dn <- as(banda3, 'SpatialGridDataFrame')
banda3.rad <- radconv(banda3.dn, 1.1548E-02, -57.74164)
banda4.dn <- as(banda4, 'SpatialGridDataFrame')
banda4.rad <- radconv(banda4.dn, 9.7382E-03, -48.69100)
banda5.dn <- as(banda5, 'SpatialGridDataFrame')
banda5.rad <- radconv(banda5.dn, 5.9593E-03, -29.79646)
banda6.dn <- as(banda6, 'SpatialGridDataFrame')
banda6.rad <- radconv(banda6.dn, 1.4820E-03, -7.41011)
banda7.dn <- as(banda7, 'SpatialGridDataFrame')
banda7.rad <- radconv(banda7.dn, 4.9952E-04, -2.49760)Reflectancia
La reflectancia
reflconvS(x, Mp, Ap, sunelev)
refb1 <- reflconvS(banda1.dn, 2.0000E-05, -0.100000, 61.24856156)
refb2 <- reflconvS(banda2.dn, 2.0000E-05, -0.100000, 61.24856156)
refb3 <- reflconvS(banda3.dn, 2.0000E-05, -0.100000, 61.24856156)
refb4 <- reflconvS(banda4.dn, 2.0000E-05, -0.100000, 61.24856156)
refb5 <- reflconvS(banda5.dn, 2.0000E-05, -0.100000, 61.24856156)
refb6 <- reflconvS(banda6.dn, 2.0000E-05, -0.100000, 61.24856156)
refb7 <- reflconvS(banda7.dn, 2.0000E-05, -0.100000, 61.24856156)Para poder ralizar la union de las bandas es necesario transformar el formato Spatialgriddataframe a formato raster refb1 <- as(refb1, ‘RasterLayer’) plot(refb1)
refb1 <- as(refb1, 'RasterLayer')
refb2 <- as(refb2, 'RasterLayer')
refb3 <- as(refb3, 'RasterLayer')
refb4 <- as(refb4, 'RasterLayer')
refb5 <- as(refb5, 'RasterLayer')
refb6 <- as(refb6, 'RasterLayer')
refb7 <- as(refb7, 'RasterLayer')Union de bandas
Union de todas las bandas de la imagen y ploteo de la imagen
bandas <- brick(refb1, refb2, refb3, refb4, refb5, refb6, refb7)
plot(bandas)
En este paso se realiza el cambio de los nombres de las bandas ya que tienen el nombre de la imagen original y al plotear no se ve bien en el grafico
names(bandas) <- paste("Banda",1:7,sep="")
plot(bandas)
Cortar imagen con la zona de estudio
Se seleccionar zona para cortar ubicandola en la ruta en la cual se encuentra
zona <- shapefile("/home/maria.orellana/Natalia_Landsat_8/Area_zonest.shp")cambio sistema de referencia zona
Es necesario saber si el sistema de coordenadas de la imagen y de la zona con la cual quiero cortar la imagen es el mismo
crs(bandas)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0crs(zona)## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0Para cambiar el sistema de coordenadas se utiliza la funcion spTrasform que es la que permite reproyectar el formato shp, en este caso particular al mismo sistema de coordenadas en que se encuentra la imagen
zona <- spTransform(zona, CRS = crs(bandas))Corte de la imagen con el zona seleccionada
bandas <- crop(bandas, extent(zona))
bandas <- mask(bandas, zona)Ploteo en color verdadero
Ploteo de la imágen cortada en color verdadero (RGB)
plotRGB(bandas, r=4, g=3, b=2, scale=10000, stretch="lin", main="RGB composite (b4,b3,b2) of Landsat-8")
Extraer los valores del raster en un dataframe
Este procedimiento se realiza para simplificar el flujo de trabajo del procesamiento
bandasDF <- values(bandas)Revisar NAs de los datos
Esto debido a que los algoritmos de agrupamiento no se ejecutan con valores NA y arrojan un error. Una forma de hacerlo es utilizando un índice lógico que hace un sub-conjunto con los datos.
idx <- complete.cases(bandasDF)Inicio del raster dataset el cual va a guardar todas las soluciones del cluster
De dos grupos de clusters hasta 12
bandasKM <- raster(bandas[[1]])
bandasCLARA <- raster(bandas[[1]])Codigo completo:
for(nClust in 2:12){
cat("-> Clustering data for nClust =",nClust,"......")
# Realización del agrupamiento o cluser de K-means
km <- kmeans(bandasDF[idx,], centers = nClust, iter.max = 50)
# Realización del agrupamiento o cluser de CLARA utilizando distancia manhattan
# La distancia manhattan esta definida como la distancia medida entre dos puntos a lo largo de los ejes en angulo recto
cla <- clara(bandasDF[idx, ], k = nClust, metric = "manhattan")
# Crear un vector entero temporal que contenga los numeros de los clusters
kmClust <- vector(mode = "integer", length = ncell(bandas))
claClust <- vector(mode = "integer", length = ncell(bandas))
# Genera el vector de los clusters de Kmeans, el cual realiza un seguimiento de los NA
kmClust[!idx] <- NA
kmClust[idx] <- km$cluster
# Genera el vector temporal de los clusters de Clara y realiza seguimiento de los NA
claClust[!idx] <- NA
claClust[idx] <- cla$clustering
# Crea el raster temporal que contiene los nuevos agrupamientos o clusters solucion
# Para K-means
tmpbandasKM <- raster(bandas[[1]])
# Para CLARA
tmpbandasCLARA <- raster(bandas[[1]])
# Establece los valores del raster con el vector cluster
# Para K-means
values(tmpbandasKM) <- kmClust
# Para CLARA
values(tmpbandasCLARA) <- claClust
# Stack del raster temporal en los finales
if(nClust==2){
bandasKM <- tmpbandasKM
bandasCLARA <- tmpbandasCLARA
}else{
bandasKM <- stack(bandasKM, tmpbandasKM)
bandasCLARA <- stack(bandasCLARA, tmpbandasCLARA)
}
cat(" done!\n\n")
}## -> Clustering data for nClust = 2 ...... done!
##
## -> Clustering data for nClust = 3 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 255508050)## done!
##
## -> Clustering data for nClust = 4 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 255508050)## 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 (= 255508050)## done!
##
## -> Clustering data for nClust = 8 ...... done!
##
## -> Clustering data for nClust = 9 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 255508050)## done!
##
## -> Clustering data for nClust = 10 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 255508050)## done!
##
## -> Clustering data for nClust = 11 ...... done!
##
## -> Clustering data for nClust = 12 ......## Warning: Quick-TRANSfer stage steps exceeded maximum (= 255508050)## done!Escribir el cluster solución de cada uno de los algoritmos
#·writeRaster(bandasKM,"/home/rus/Taller 1 PRA/Kmeans.tif", overwrite=TRUE)
#writeRaster(bandasCLARA,"/home/rus/Taller 1 PRA/clara.tif", overwrite=TRUE)Evaluación de la clasificación no supervisada. Rendimiento de los clusters de agrupamiento
Para evaluar el rendimiento de cada solución de clúster y seleccionar el “mejor” número de clústeres, para particionar los datos de muestra, se utiliza el índice índice de silueta.
Este índice hace referencica a un método de interpretación y validación de la consistencia dentro de grupos de datos el cual se realiza utilizando el paquete clusterCrit que se denomina criterio interno. Este método realiza una representación gráfica que muestra qué tan bien se encuentra cada objeto agrupado dentro de su grupo. El valor de dicha silueta es una medida de cuán similar es un objeto para su propio clúster ya que evalua la coh cohesión dentro del clúster, en comparación con otros clústeres, mostrando la separación entre clustes.
El índice de silueta varía de -1 a +1, donde un valor alto indica que el objeto está bien adaptado a su propio clúster pero no está bien adaptado a los clústeres vecinos. Cuando la mayoría de los objetos tienen un valor alto, la configuración de agrupamiento se considera apropiada. Por el contrario si muchos puntos tienen un valor bajo o negativo, la configuración del clúster puede tener demasiados o muy pocos clústeres.
El paquete clusterCrit en R, permite implementar estos criterios de clúster internos en el intCriteria (entre muchos otros, como los índices Dunn, Ball-Hall, Davies-Bouldin, GDI, Tau).
Un detalle importante es que el cálculo del índice de silueta es un proceso bastante lento para grandes números de observaciones (> 5000), por tanto se utilizará un enfoque de muestreo aleatorio estratificado. Lo cual quiere decir que se toma un subconjunto de celdas de cada clúster y se calcula el índice en función de estas. Con esto se asumie que la muestra es algo robusta y representativa de las células completas. Lo ideal seria que el proceso se repetiera varias veces y luego calcular un valor promedio. Sin embargo, para simplificar el proceso se utiliza una sola muestra de celdas en este notbook, adicionalmente porque la esmimacion arroja errores relativamente bajos.
Funcionamiento del código:
clustPerfSI <- data.frame(nClust = 2:12, SI_KM =NA, SI_CLARA=NA)clustPerfSI <- na.omit(clustPerfSI)
for(i in 1:nlayers(bandasKM)){ # Iteración a traves de cada capa
cat("-> Evaluación de rendimiento de los clusters por cluster =",(2:12)[i],"......")
# Extraccion de muestras de celulas aleatorias estratificadas por clustr
cellIdx_bandasKM <- sampleStratified(bandasKM[[i]], size = 2000)
cellIdx_bandasCLARA <- sampleStratified(bandasCLARA[[i]], size = 2000)
# Obtiene los valores de celda de la muestra aleatoria estratificada del ráster
# objeto data frame object (bandasDF)
bandasDFStRS_KM <- bandasDF[cellIdx_bandasKM[,1], ]
bandasDFStRS_CLARA <- bandasDF[cellIdx_bandasCLARA[,1], ]
# Para asegúrese que todas las columnas sean numéricas ya que la función intCriteria lo exige.
bandasDFStRS_KM[] <- sapply(bandasDFStRS_KM, as.numeric)
bandasDFStRS_CLARA[] <- sapply(bandasDFStRS_CLARA, as.numeric)
# Calculo del índice de siluetas basado en muestras para:
#
# K-means
clCritKM <- intCriteria(traj = bandasDFStRS_KM,
part = as.integer(cellIdx_bandasKM[,2]),
crit = "Silhouette")
# CLARA
clCritCLARA <- intCriteria(traj = bandasDFStRS_CLARA,
part = as.integer(cellIdx_bandasCLARA[,2]),
crit = "Silhouette")
# Escribir los valores del índice de silueta en el data frame clustPerfSI.
# Resultados
clustPerfSI[i, "SI_KM"] <- clCritKM[[1]][1]
clustPerfSI[i, "SI_CLARA"] <- clCritCLARA[[1]][1]
cat(" done!\n\n")
}## -> Evaluación de rendimiento de los clusters por cluster = 2 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 3 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 4 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 5 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 6 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 7 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 8 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 9 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 10 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 11 ...... done!
##
## -> Evaluación de rendimiento de los clusters por cluster = 12 ...... done!write.csv(clustPerfSI, file = "/home/maria.orellana/Natalia_Landsat_8/clustPerfSI.csv", row.names = FALSE)Ploteo de la tabla con los resultados del índice de silueta para poder comparar cada solución de los cluster
knitr::kable(clustPerfSI, digits = 3, align = "c",
col.names = c("#clusters","Promedio (k-means)","Promedio (Clara)"))| #clusters | Promedio (k-means) | Promedio (Clara) |
|---|---|---|
| 2 | 0.314 | 0.299 |
| 3 | 0.371 | 0.174 |
| 4 | 0.356 | 0.241 |
| 5 | 0.346 | 0.303 |
| 6 | 0.331 | 0.274 |
| 7 | 0.326 | 0.271 |
| 8 | 0.317 | 0.246 |
| 9 | 0.316 | 0.286 |
| 10 | 0.303 | 0.277 |
| 11 | 0.319 | 0.271 |
| 12 | 0.309 | 0.266 |
La tabla anterior muestra la distancia de cada uno de los elementos del grupo a su centroide, donde se pueden diferenciar entre valores altos o bajos indicando si existe una mayor o menor cohesion del grupo respectivamente.
Ploteo para comparar los dos algoritmos:
Indice de silueta
En ete indice se muestra la cohesion y separacion de un punto x del conjunto de datos. La cohesion hace refrencia a la distancia promedio de x a los demas puntos de un mismo cluster. La separación hace referencia a la distancia promedio de x a todos los puntos del cluster mas cercano. Los valores del indice de silueta pueden variar entre -1 y 1, siendo -1 un mal agrupamiento, 0 indiferente y 1 buen agrupamiento.
plot(clustPerfSI[,1], clustPerfSI[,2],
xlim = c(1,13), ylim = range(clustPerfSI[,2:3]), type = "n",
ylab="Promedio Indice de silueta", xlab="# de Clusters",
main="Indice de silueta por # de clusters")
## Ploteo del promedio de los valores del índice de silueta vs el # de clusters para K-means
lines(clustPerfSI[,1], clustPerfSI[,2], col="red")
## Ploteo del promedio de los valores del índice de silueta vs el # de clusters para CLARA
lines(clustPerfSI[,1], clustPerfSI[,3], col="blue")
## Lineas de la grilla
abline(v = 1:13, lty=2, col="light grey")
abline(h = seq(0.30,0.44,0.02), lty=2, col="light grey")
## Leyenda
legend("topright", legend=c("K-means","CLARA"), col=c("red","blue"), lty=1, lwd=1)
Teniendo en cuenta los valores entre los que varia el indice de silueta, en la tabla y gráfica anterior los valores mas altos son los que presentan un mejor agrupamiento.
Se evidencian resultados diferentes en cuanto al rendimiento de agrupamiento de los dos algoritmos utilizados. El algoritmo K-means muestra un mejor rendimiento respecto a CLARA lo cual puede atribuirse a que este ulimo trabaja en un subconjunto de datos, por lo tanto, es mas difícil encontrar los mejores centros de clúster. Para el algoritmo K-means, la partición de los datos en 4 grupos / clusters parece ser la mejor opción.
Diagrama de las mejores soluciones según el índice de silueta:
plot(bandasKM[[4]])
Conclusiones
- Por medio de este ejercicio se puede evidenciar que el usuario no tiene ninguna manipulacion de información, los unicos datos que se proporcionan al algoritmo es el numero de clusters.
- Por su parte el algoritmo realiza todas las operaciones necesarias arrojando un resultado final de clasificación.
- Para este ejercicio en particular el algoritmo Kmeans fue el que presentó mejores rendimientos.
Discusión
Este ejercicio permite desarrollar habilidades en el manejo de R permitiendo obtener resultados que cumplen el objetivo propuesto. A pesar de que los algoritmos utilizados realizan todo el proceso para la clasificacion de la imagen, el paso a paso permite entender el desarrollo que se está llevando a cabo en cada uno de estos. Al igual que los algoritmos utilizados asi como otros algoritmos, que para el caso particular presentan algunas diferencias en cuanto a rendimiento, se puede considerar que realizan un buen agrupamiento de datos.
Bibliografía
Blaschke, T., Kelly, M y Merschdorf, H. 2015. Object based Image Analysis: Evolution, History, State of the Art and Future Vison. Brownlee, J. 2016. Supervised and Unsupervised Machine Learning Algorithms. www.machinelearingmastery.com Humboldt State University, 2015. GSP 216. Introduction to remote sensing. Humboldt State Geospatil online.