library(cluster)
library(spgrass6)

x <- readRAST6('e10_sd')

n_classes <- 5
x.clara <- clara(data.frame(coordinates(x),  slot(x, 'data')), k=n_classes, stand=TRUE)

x$clust <- factor(x.clara$clustering)
spplot(x, zcol='clust', col.regions=terrain.colors(n_classes))

# pdf(file='elev_sd_clust.pdf', height=5, width=5)
boxplot(x$e10_sd ~ factor(x$clust, levels=order(tapply(x$e10_sd, x$clust, median)) ), ylab='Standard Deviation (meters)', xlab='Class', main='Standard Deviation of Elevation on 50x50 Grid', cex=0.5 )
# dev.off()

x$clust <- x.clara$clustering
writeRAST6(x, zcol='clust', vname='t_clust')