Nienadzorowana klasyfikacja z kilometrami w R

Mam szereg zdjęć satelitarnych (5 pasm) i chcę je sklasyfikować według kmeanów w R. Mój skrypt działa dobrze (przeglądaj moje obrazy, konwertuj obrazy na data.frame, grupuj je i konwertuj z powrotem na 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)
}

Mój problem brzmi: nie mogę porównać wyników klasyfikacji, ponieważ przypisania klastra różnią się między obrazami. Na przykład „woda” znajduje się w pierwszym klastrze obrazów nr 1, w kolejnych 2 i w trzeciej 10, co uniemożliwia porównanie wyników dotyczących wody między datami.

Jak mogę naprawić przypisanie klastra?

Czy mogę określić stały punkt początkowy dla wszystkich obrazów (mając nadzieję, że woda jest zawsze wykrywana jako pierwsza, a zatem klasyfikowana jako 1)?

A jeśli tak, to w jaki sposób?

Myślę, że nie możesz ... Musisz najpierw oznaczyć wszystkie klasy, aby je porównać. Kmean klasyfikuje bez nadzoru, więc bez wcześniejszych informacji i dlatego nie może zdefiniować żadnego rodzaju klas.

Jeśli masz warstwę odniesienia, możesz wykonać etykietowanie większością głosów. Oto znacznie bardziej wydajny kod do głosowania większością głosów niż przy użyciu funkcji pakietu „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]

gdzie refjest twój plik referencyjny klasy rastrowej, class_fileto wynik twojego kmeans.

zr daje ci w pierwszej kolumnie numer „strefy”, aw drugiej kolumnie etykietę dla klasy.

Aby zaimplementować klastrowanie na stosie obrazów, nie robi się tego pasmo po pasmie, a raczej na całym stosie obrazów jednocześnie. W przeciwnym razie, jak zauważył @nmatton, statystyka nie ma większego sensu.

Nie zgadzam się jednak z tym, że nie jest to możliwe, po prostu wymaga dużej ilości pamięci. W przypadku rzeczywistych danych satelitarnych będzie to ogromny problem, a być może niemożliwe w przypadku danych o wysokiej rozdzielczości, ale możesz przetwarzać w pamięci, zmuszając swoje rastry do jednego obiektu, który można przekazać do funkcji klastrowania. Będziesz musiał śledzić wartości NA między rastrami, ponieważ zostaną one usunięte podczas grupowania i będziesz musiał znać pozycje w rastrze, abyś mógł przypisać wartości skupienia do właściwych komórek.

Tutaj możemy przejść przez jedno podejście. Dodajmy wymagane biblioteki i przykładowe dane (logo RGB R daje nam 3 pasma do pracy).

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

Po pierwsze, możemy zmusić nasz wielopasmowy obiekt stosu rastrowego do data.frame za pomocą getValues. Zauważ, że dodaję wartość NA w wierszu 1, kolumnie 3, abym mógł zilustrować sposób postępowania bez danych.

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

Tutaj możemy przejść do biznesu i utworzyć indeks komórek wartości innych niż NA, który zostanie użyty do przypisania wyników klastra.

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

Teraz tworzymy obiekt skupienia z 3 pasmowych wartości RGB o k = 4. Używam metody clara K-Medoids, ponieważ jest dobra w przypadku dużych danych i lepsza w przypadku nieparzystych dystrybucji. Jest bardzo podobny do K-średnich.

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

Dla uproszczenia możemy utworzyć pusty raster, wyciągając jeden z pasm rastrowych z naszego oryginalnego obiektu stosu rastrowego i przypisując mu wartości NA.

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

Na koniec, korzystając z indeksu, przypisujemy wartości klastra do odpowiedniej komórki w pustym rastrze i wykreślamy wyniki.

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

W przypadku dużych rastrów możesz zajrzeć do pakietu bigmemory, który zapisuje macierze na dysku i działa na blokach oraz dostępna jest funkcja k-średnich. Należy również pamiętać, że nie do tego dokładnie służy R i że oprogramowanie do przetwarzania obrazu lub GIS może być bardziej odpowiednie. Wiem, że SAGA i przybornik Orfeo to wolne oprogramowanie, w którym klastry k-średnich są dostępne dla stosów obrazów. Istnieje nawet biblioteka RSAGA, która pozwala na wywoływanie oprogramowania z R.