Classificação não supervisionada com kmeans em R

Eu tenho uma série temporal de imagens de satélite (5 bandas) e quero classificá-las por kmeans em R. Meu script está funcionando bem (faça um loop pelas minhas imagens, converta as imagens em data.frame, agrupe-as e converta-as novamente em varredura):

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)
}

Meu problema é: não consigo comparar os resultados da classificação entre si, porque os responsáveis ​​pelo cluster diferem de imagem para imagem. Por exemplo, "água" está no primeiro cluster de imagens número 1, nos próximos 2 e nos terceiros 10, tornando impossível comparar os resultados da água entre as datas.

Como posso corrigir a atribuição de cluster?

Posso especificar um ponto de partida fixo para toda a imagem (esperando que a água seja sempre detectada primeiro e, portanto, classificada como 1)?

E se sim, como?

Eu acho que você não pode ... Você primeiro precisa rotular cada classe para compará-las. O Kmean classifica sem supervisão, sem informações prévias e, portanto, não pode definir nenhum tipo de classe.

Se você tiver uma camada de referência, poderá fazer uma rotulagem por votação majoritária. Aqui está um código muito mais eficiente para a votação majoritária do que usar a função de pacote '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]

onde refestá o arquivo de referência da classe raster, class_fileé o resultado do kmeans.

zr fornece a você na primeira coluna o número da 'zona' e na segunda coluna o rótulo da classe.

Para implementar o clustering em uma pilha de imagens, você não faz isso banda por banda, mas sim na pilha inteira de imagens simultaneamente. Caso contrário, como apontado por @nmatton, a estatística não faz muito sentido.

No entanto, não concordo que isso não seja possível, apenas memória intensiva. Em dados reais de satélite, esse será um grande problema, e talvez impossível em dados de alta resolução, mas você pode processar na memória coagindo suas varreduras em um único objeto que pode ser passado para uma função de agrupamento. Você precisará rastrear os valores de NA entre os rasters, pois eles serão removidos durante o cluster e você precisará conhecer as posições na varredura para poder atribuir os valores do cluster às células corretas.

Podemos passar por uma abordagem aqui. Vamos adicionar as bibliotecas necessárias e alguns dados de exemplo (o logotipo RGB R para fornecer três bandas para trabalhar).

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

Primeiro, podemos coagir nosso objeto de pilha de varredura de várias bandas a um data.frame usando getValues. Observe que estou adicionando um valor de NA na linha 1, coluna 3, para ilustrar como lidar com nenhum dado.

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

Aqui, podemos começar a trabalhar e criar um índice de célula com valores que não sejam de NA que serão usados ​​para atribuir os resultados do cluster.

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

Agora, criamos um objeto de cluster a partir dos valores RGB de 3 bandas com k = 4. Estou usando o método claro K-Medoids porque é bom com grandes dados e melhor com distribuições ímpares. É muito semelhante ao K-Means.

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

Por uma questão de simplicidade, podemos criar uma varredura vazia, puxando uma das faixas de varredura de nosso objeto de pilha de varredura original e atribuindo valores de NA.

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

Finalmente, usando o índice, atribuímos os valores de cluster à célula apropriada na varredura vazia e plotamos os resultados.

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

Para grandes usuários de rasters, você pode procurar no pacote bigmemory, que grava matrizes em disco e opera em blocos, e há uma função k-means disponível. Além disso, lembre-se de que não foi exatamente para isso que o R foi projetado e que um software de processamento de imagem ou GIS pode ser mais apropriado. Eu sei que o SAGA e a caixa de ferramentas Orfeo são softwares livres que possuem agrupamentos k-means disponíveis para pilhas de imagens. Existe até uma biblioteca RSAGA que permite que o software seja chamado de R.