Clasificación no supervisada con kmeans en R

10

Tengo una serie temporal de imágenes de satélite (5 bandas) y quiero clasificarlas por kmeans en R. Mi script está funcionando bien (recorrer mis imágenes, convertir las imágenes en data.frame, agruparlas y convertirlas nuevamente en un trama):

for (n in files) {
image <- stack(n)    
image <- clip(image,subset)

###classify raster
image.df <- as.data.frame(image)  
cluster.image <- kmeans(na.omit(image.df), 10, iter.max = 10, nstart = 25) ### kmeans, with 10 clusters

#add back NAs using the NAs in band 1 (identic NA positions in all bands), see http://stackoverflow.com/questions/12006366/add-back-nas-after-removing-them/12006502#12006502
image.df.factor <- rep(NA, length(image.df[,1]))
image.df.factor[!is.na(image.df[,1])] <- cluster.image$cluster

#create raster output
clusters <- raster(image)   ## create an empty raster with same extent than "image"  
clusters <- setValues(clusters, image.df.factor) ## fill the empty raster with the class results  
plot(clusters)
}

Mi problema es: no puedo comparar los resultados de clasificación entre sí porque los asignadores de clúster difieren de una imagen a otra. Por ejemplo, "agua" está en el primer grupo de imágenes número 1, en los siguientes 2 y en los terceros 10, lo que hace imposible comparar los resultados del agua entre las fechas.

¿Cómo puedo arreglar la asignación del clúster?

¿Puedo especificar un punto de partida fijo para todas las imágenes (esperando que el agua siempre se detecte primero y, por lo tanto, se clasifique como 1)?

Y si es así, ¿cómo?

Iris
fuente

Respuestas:

6

Creo que no puedes ... Primero debes etiquetar cada clase para compararlas. Kmean clasifica sin supervisión, sin ninguna información previa y, por lo tanto, no puede definir ningún tipo de clase.

Si tiene una capa de referencia, puede hacer un etiquetado por mayoría de votos. Aquí hay un código bastante más eficiente para la votación por mayoría que usar la función de paquete 'raster' zonal:

require (data.table)
fun <- match.fun(modal)
vals <- getValues(ref) 
zones <- round(getValues(class_file), digits = 0) 
rDT <- data.table(vals, z=zones) 
setkey(rDT, z) 
zr<-rDT[, lapply(.SD, modal,na.rm=T), by=z]

donde refestá el archivo de referencia de la clase de ráster, class_filees el resultado de kmeans.

zr le da en la primera columna el número de 'zona' y en la segunda columna, la etiqueta de la clase.

nmatton
fuente
Tenía miedo de que no sea posible. ¡Gracias por el código para la votación por mayoría!
Iris
4

Para implementar la agrupación en una pila de imágenes, no lo hace banda por banda sino más bien en toda la pila de imágenes simultáneamente. De lo contrario, como señaló @nmatton, la estadística no tiene mucho sentido.

Sin embargo, no estoy de acuerdo en que esto no sea posible, solo en memoria intensiva. En datos satelitales reales, este será un gran problema, y ​​tal vez imposible en datos de alta resolución, pero puede procesar en la memoria coaccionando su (s) ráster (es) en un solo objeto que se puede pasar a una función de agrupamiento. Deberá realizar un seguimiento de los valores de NA en los rásteres porque se eliminarán durante el agrupamiento y deberá conocer las posiciones en el ráster para poder asignar los valores del clúster a las celdas correctas.

Podemos pasar por un enfoque aquí. Agreguemos las bibliotecas requeridas y algunos datos de ejemplo (el logotipo RGB R para darnos 3 bandas para trabajar).

library(raster)
library(cluster)
r <- stack(system.file("external/rlogo.grd", package="raster")) 
  plot(r)

Primero, podemos obligar a nuestro objeto de pila de ráster multibanda a un marco de datos usando getValues. Tenga en cuenta que estoy agregando un valor de NA en la fila 1, columna 3 para que pueda ilustrar cómo tratar sin datos.

r.vals <- getValues(r[[1:3]])
  r.vals[1,][3] <- NA

Aquí, podemos ponernos manos a la obra y crear un índice de celda de los valores distintos de NA que se utilizarán para asignar los resultados del clúster.

idx <- 1:ncell(r)
idx <- idx[-unique(which(is.na(r.vals), arr.ind=TRUE)[,1])]  

Ahora, creamos un objeto de clúster a partir de los valores RGB de 3 bandas con k = 4. Estoy usando el método clara K-Medoids porque es bueno con datos grandes y es mejor con distribuciones extrañas. Es muy similar a K-Means.

clus <- cluster::clara(na.omit(scale(r.vals)), k=4)

Por simplicidad, podemos crear un ráster vacío extrayendo una de las bandas de ráster de nuestro objeto de pila ráster original y asignándole valores NA.

r.clust <- r[[1]]
r.clust[] <- NA

Finalmente, usando el índice, asignamos los valores del clúster a la celda apropiada en el ráster vacío y graficamos los resultados.

r.clust[idx] <- clus$clustering
plot(r.clust) 

Para rásters enormes, es posible que desee buscar en el paquete bigmemory que escribe matrices en el disco y opera en bloques y hay una función k-means disponible. Además, tenga en cuenta que esto no es exactamente para lo que R fue diseñado y que un procesamiento de imágenes o software SIG puede ser más apropiado. Sé que SAGA y la caja de herramientas Orfeo son software libre que tiene clústeres de k-means disponibles para pilas de imágenes. Incluso hay una biblioteca RSAGA que permite llamar al software desde R.

Jeffrey Evans
fuente
Si todas las imágenes se apilan y agrupan a la vez, entonces, el resultado es una imagen agrupada, ¿verdad?
Iris
@ Iris, sí, así es como funciona este tipo de agrupación de imágenes y sigue las implementaciones en el software de detección remota. Un ejemplo claro y relevante sería la implementación de isocluster en ArcGIS ( desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/… )
Jeffrey Evans
Entonces esta respuesta no ayuda en absoluto. Mi problema fue que intenté hacer una detección de cambios con el tiempo en función de varias clasificaciones de imágenes sin supervisión, pero pude comparar los diferentes resultados porque las clases se asignaron de manera diferente.
Iris
La clasificación no supervisada no es una forma viable de realizar la detección de cambios. Incluso una ligera variación en una imagen dada podría terminar con píxeles asignados a una clase diferente. Este sería el caso incluso si proporcionara centros de clúster para K-Means. Tengo una función de entropía en el paquete spatialEco que es útil para la detección de cambios. Calcula la entropía dentro de una ventana NxN y luego deriva delta en cada paso de tiempo. La entropía negativa representa la pérdida y la ganancia positiva de los componentes del paisaje dentro de una magnitud dada bajo la entropía máxima.
Jeffrey Evans el
Esa es una vieja pregunta y descarté la idea de usar k-means hace años. Pero es bueno saber el paquete SpacialEco para la próxima vez;)
Iris