I have a time series of satellite images (5 bands) and want to classify them by kmeans in R. My script is working fine (loop through my images, convert the images to data.frame, cluster them, and convert it back to a 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  

My problem is: I can't compare the classification results to each other because the cluster assignents differ from image to image. For example, "water" is in the first image cluster number 1, in the next 2 and in the third 10, making it impossible to compare the water results between the dates.

How can I fix the cluster assignment?

Can I specify a fixed starting point for all image (hoping that water is always detected first and thus classified as 1)?

And if yes, how?

2 Answers 2


I think you can't... You first have to label each classes to compare them. Kmean classify unsupervisedly so without any prior information and so cannot define any kind of classes.

If you have a reference layer, you can make a labelling by a majority voting. Here's a quite more efficient code for majority voting than using the 'raster' package function 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]

where ref is your raster class reference file, class_file is your kmeans result.

zr gives you in first col the 'zone' number and in second col, the label for the class.

  • I was afraid that it's not possible. Thank you for the code for majority voting!
    – Iris
    Commented Dec 1, 2014 at 9:28

To implement clustering on an image stack, you do not do it band-by-band but rather on the entire image stack simultaneously. Otherwise, as pointed out by @nmatton, the statistic does not make much sense.

However, I do not agree that this is not possible, just memory intensive. On real satellite data this will be a huge problem, and perhaps impossible on high resolution data, but you can process in memory by coercing your raster(s) into a single object that can be passed to a clustering function. You will need to track NA values across rasters because they will be removed during clustering and you will need to know the positions in the raster so you can assign the cluster values to the correct cells.

We can step through one approach here. Lets add the required libraries and some example data (the RGB R logo to give us 3 bands to work with).

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

First, We can coerce our multi-band raster stack object to a data.frame using getValues. Note that I am adding an NA value at row 1, column 3 so I can illustrate how to deal with no data.

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

Here, we can get down to business and create a cell index of the non-NA values that will be used to assign the cluster results.

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

Now, we create a cluster object from the 3 band RGB values with k=4. I am using the clara K-Medoids method because it is good with large data and is better with odd distributions. It is very similar to K-Means.

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

For simplicity sake, we can create an empty raster by pulling one of the raster bands from our original raster stack object and assigning it NA values.

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

Finally, using the index, we assign the cluster values to the appropriate cell in the empty raster and plot the results.

r.clust[idx] <- clus$clustering

For huge rasters you may want to look into the bigmemory package which writes matrices to disk and the operates on blocks and there is a k-means function available. Also, keep in mind that this not exactly what R was designed for and that an image processing or GIS software may be more appropriate. I know that SAGA and the Orfeo toolbox are both free software that have k-means clustering available for image stacks. There is even an RSAGA library that allows the software to be called from R.

  • If all images are stacked and clustered at once, then, the result is one clustered image, right?
    – Iris
    Commented Sep 15, 2016 at 8:11
  • @Iris, yes this is how this type of image clustering works and follows the implementations in remote sensing software. A clear and relevant example would be the isocluster implementation in ArcGIS (desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/…) Commented Sep 15, 2016 at 16:31
  • Then this anwer does not help at all. My problem was that I tried to do a change detection over time based on several unsupervied image classifications, but I could compare the different results because the classes were assigned differently.
    – Iris
    Commented Sep 15, 2016 at 16:48
  • Unsupervised classification is not a viable way to perform change detection. Even slight variation in a given image could end up with pixels being assigned into a different class. This would be the case even if you provided cluster centers for K-Means. I have an entropy function in the spatialEco package that is useful for change detection. You calculate the entropy within an NxN window and then derive delta at each time-step. Negative entropy represents loss and positive is gain of landscape components within a given magnitude under the maximum entropy. Commented Sep 15, 2016 at 17:03
  • That is an old question and I discared the idea of using k-means ages ago. But good to know the spatialEco package for next time ;)
    – Iris
    Commented Sep 15, 2016 at 20:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.