Classement non supervisé avec kmeans en R

J'ai une série temporelle d'images satellite (5 bandes) et je souhaite les classer par kmeans dans R. Mon script fonctionne correctement (parcourir mes images, convertir les images en data.frame, les regrouper et les reconvertir en un raster):

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

Mon problème est le suivant: je ne peux pas comparer les résultats de la classification les uns aux autres car les affectations de cluster diffèrent d'une image à l'autre. Par exemple, "eau" est dans le premier groupe d'images numéro 1, dans les 2 suivants et dans le troisième 10, ce qui rend impossible la comparaison des résultats de l'eau entre les dates.

Comment puis-je corriger l'affectation du cluster?

Puis-je spécifier un point de départ fixe pour toutes les images (en espérant que l'eau est toujours détectée en premier et donc classée comme 1)?

Et si oui, comment?

Je pense que vous ne pouvez pas ... Vous devez d'abord étiqueter chaque classe pour les comparer. Kmean classe sans supervision donc sans aucune information préalable et ne peut donc définir aucun type de classes.

Si vous avez une couche de référence, vous pouvez faire un étiquetage par un vote majoritaire. Voici un code bien plus efficace pour voter à la majorité que d'utiliser la fonction de package «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]

où refest votre fichier de référence de classe raster, class_fileest votre résultat kmeans.

zr vous donne dans le premier col le numéro de «zone» et dans le deuxième col, l'étiquette de la classe.

Pour implémenter le clustering sur une pile d'images, vous ne le faites pas bande par bande mais plutôt sur la pile d'images entière simultanément. Sinon, comme l'a souligné @nmatton, la statistique n'a pas beaucoup de sens.

Cependant, je ne suis pas d'accord pour dire que ce n'est pas possible, juste une mémoire intensive. Sur les données satellite réelles, ce sera un énorme problème, et peut-être impossible sur les données haute résolution, mais vous pouvez traiter en mémoire en contraignant votre raster (s) en un seul objet qui peut être transmis à une fonction de clustering. Vous devrez suivre les valeurs NA sur les rasters car elles seront supprimées lors du clustering et vous devrez connaître les positions dans le raster pour pouvoir affecter les valeurs du cluster aux cellules correctes.

Nous pouvons passer par une approche ici. Ajoutons les bibliothèques requises et quelques exemples de données (le logo RGB R pour nous donner 3 bandes avec lesquelles travailler).

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

Tout d'abord, nous pouvons contraindre notre objet de pile raster multibande à un data.frame en utilisant getValues. Notez que j'ajoute une valeur NA à la ligne 1, colonne 3 afin que je puisse illustrer comment traiter sans données.

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

Ici, nous pouvons nous mettre au travail et créer un index de cellule des valeurs non-NA qui seront utilisées pour affecter les résultats du cluster.

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

Maintenant, nous créons un objet cluster à partir des valeurs RVB 3 bandes avec k = 4. J'utilise la méthode clara K-Medoids car elle est bonne avec de grandes données et est meilleure avec des distributions impaires. Il est très similaire à K-Means.

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

Par souci de simplicité, nous pouvons créer un raster vide en tirant l'une des bandes raster de notre objet de pile raster d'origine et en lui affectant des valeurs NA.

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

Enfin, à l'aide de l'index, nous attribuons les valeurs de cluster à la cellule appropriée dans le raster vide et traçons les résultats.

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

Pour les rasters énormes, vous voudrez peut-être examiner le package bigmemory qui écrit les matrices sur le disque et les opère sur les blocs et il existe une fonction k-means. Gardez également à l'esprit que ce n'est pas exactement pour cela que R a été conçu et qu'un logiciel de traitement d'image ou SIG peut être plus approprié. Je sais que SAGA et la boîte à outils Orfeo sont tous deux des logiciels gratuits qui ont un clustering k-means disponible pour les piles d'images. Il existe même une bibliothèque RSAGA qui permet d'appeler le logiciel depuis R.