Search code examples
rclassificationgisk-meansraster

Creating and setting clusters values to a new raster- k-mean classification


I want to create a simple kmeans unsupervised classification. I have a problem with creating clusters and setting the cluster values to a new raster. I was inspired on this site.

landsat5 <- stack('5.tif')
names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'thermal', 'SWIR2')
ndvi <- (landsat5[['NIR']] - landsat5[['red']]) / (landsat5[['NIR']] + landsat5[['red']])
nr <- getValues(ndvi)
set.seed(99)
kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")
knr <- setValues(ndvi, kmncluster$cluster)
knr <- raster(ndvi)
values(knr) <- kmncluster$cluster
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
             "#c3ff5b", "#ff7373", "#00ff00", "#808080")
par(mfrow = c(1,2))
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(knr, main = 'Unsupervised classification', col = mycolor )

I get this error :

Error in 'setValues(ndvi, kmncluster$cluster)': length(values) is not equal to ncell(x), or to 1.

Any ideas why it shows the error?


Solution

  • As s_t points out, the problem is created by na.omit(nr) because that removes cases such that the number of cases is no longer equal to the number of raster cells.

    Here is a minimal, reproducible, self-contained example

    library(raster)
    b <- brick(system.file("external/rlogo.grd", package="raster"))
    vi <- (b$red - b$green) / (b$red + b$green)
    nr <- getValues(vi)
    

    There are NA values in nr (where red+green == 0) and these need to be removed to be able to use kmeans. But instead of using na.omit you can do

    i <- !is.na(nr)
    set.seed(99)
    kmncluster <- kmeans(nr[i], centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")
    

    And now you can replace the values in nr with their cluster memberships and put these values into a RasterLayer

    nr[i] <- kmncluster$cluster
    knr <- setValues(vi, nr)
    plot(knr)