predict for kmeans

[aloboa]

Applying kmeans to all pixels of a large multispectral image is not realistic.
A better way would be to apply kmeans to a random sample and then have
a method for predict() (or an specific predict()) to assign each pixel to its most
similar class.
Could this be implemented?

[kadyb]

I think you should use the kmeans() function which has a predict method, e.g. clue::cl_predict() or ClusterR::predict_KMeans(). You can also write prediction method for the stats::kmeans() function: https://stackoverflow.com/questions/53352409/creation-prediction-function-for-kmean-in-r/53352914#53352914

[aloboa]

Is there a terra::kmeans() function? The one I see is stats::means().

[kadyb]

Is there a terra::kmeans() function?

No, but you will find plenty of other packages for clustering on CRAN Task View.

Here is example with stats::kmeans() (instead of defined method, you can use {clue} or {ClusteR} package):

library("terra")

predict.kmeans = function(x, newdata) {
  apply(newdata, 1, function(r) which.min(colSums((t(x$centers) - r)^2)))
}
  
f = system.file("ex/logo.tif", package = "terra")
r = rast(f)
smp = spatSample(r, size = 30)
mdl = kmeans(smp, centers = 3)
pr = predict(r, mdl)
pr

[aloboa]

Perfect, thanks.
My request is adding such a predict.kmeans() function as an standard function in terra.

[kadyb]

Personally, I think it would be better to add predict method upstream to R, then all packages could benefit, e.g. {raster}, {stars}, etc. The question remains why add it to {terra} only, if the function is already available in several packages / we can write our own function.

[rhijmans]

I have added a k_means method (could not use kmeans the way I would like to because it has no ellipses). I think it is worthwhile having as it is so commonly used, and because dealing with NAs is a bit cumbersome.

[aloboa]

Cool, thanks.

[kadyb]

@rhijmans, maybe you could consider using loop with vector allocation instead of apply() -- it should run a little faster and use less memory.

predict_1 = function(x, newdata) {
  apply(newdata, 1, function(r) which.min(colSums((t(x$centers) - r)^2)))
}

predict_2 = function(x, newdata) {
  vec = integer(nrow(newdata))
  newdata = as.matrix(newdata)
  x = x$centers
  for (i in seq_len(nrow(newdata))) {
    vec[i] = which.min(colSums((t(x) - newdata[i, ])^2))
  }
  return(vec)
}

n = 1e6
df = data.frame(x = rnorm(n), y = rnorm(n), z = rnorm(n))
mdl = kmeans(df, centers = 5)

system.time(predict_1(mdl, df))
#>    user  system elapsed
#>   11.09    0.33   11.42
system.time(predict_2(mdl, df))
#>    user  system elapsed
#>    7.91    0.11    8.01

[Nowosad]

@rhijmans the k_means() function does not work correctly:

library(terra)
#> terra 1.7.67
f <- system.file("ex/logo.tif", package = "terra")
r <- rast(f)
km <- k_means(r, centers=5)
#> Warning: [setValues] values were recycled
km
#> class       : SpatRaster 
#> dimensions  : 77, 101, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#> coord. ref. : Cartesian (Meter) 
#> source(s)   : memory
#> name        : lyr1 
#> min value   :    1 
#> max value   :    5
plot(km)

[rhijmans]

Thank you. Fixed.

[Nowosad]

It now works for rasters without NAs, but when NAs do exist:

library(terra)
#> terra 1.7.67
f <- system.file("ex/logo.tif", package = "terra")
r <- rast(f)
r[1] <- NA
km <- k_means(r, centers=5)
#> Error in do_one(nmeth): NA/NaN/Inf in foreign function call (arg 1)

[rhijmans]

This is embarrassing, but thank you for your patience. Looks like it works again.

[Nowosad]

Robert -- thank you a lot for your work!