Using meanshift function with raster data in R

Question

I have been playing around with different functions in R and I wanted to use the meanshift's msClustering() function with my Raster data.

Following I have a reproducible example:

a<-matrix(data=sample(1:500,25),nrow = 5, ncol = 5) b<-raster::raster(a) raster::crs(b)<-"+init=epsg:4624" c<-MeanShift::msClustering(a) b$new<-c$components

and I get an error saying

1: In v[] <- value : number of items to replace is not a multiple of replacement length

which is logical given that they don't have the same length or dimensions.

So based on that example and error can anyone explain to me why I get a matrix of different dimension and give me a solution? Because I think my problem is more general and theoretical based on the classification method, rather than technical, so I figured it would be better to ask here.

Answer 1

****Edit (05.15.2018), please see revised code. You can, in fact, assign the components object to the raster. Sorry I did not see this in the first place. ****

You need to pipe the new values into the raster using a bracket index or raster::setValues.

library(raster)
a <- matrix(data=sample(1:500,25),nrow = 5, ncol = 5)
clust <- raster::raster(a)
component <- raster::raster(a)
cl <- MeanShift::msClustering(a/sd(a))
clust[] <- cl$labels
component[] <- cl$components[,2] # mode 2
  plot(stack(clust,component))

The components object is not a vector but, rather a matrix and does not match the dimensions of the raster due to an extra column. The components are a matrix of the cluster representatives values of the mean shift Gaussian mixtures and not the cluster values, which are contained in labels. If you exclude the first column of the components matrix, the matrix dimensions match. Please take note that for large rasters this is not memory safe and could explode in RAM quickly, you are effectively storing the raster twice.

Answer 2

If you want to apply this algorithm in raster objects, you need to create a matrix with layers x pixels shape. Here a basic example for satellite data:

library(raster)
library(MeanShift)
library(RStoolbox)

mtlFile  <- system.file("external/landsat/LT52240631988227CUB02_MTL.txt", 
                        package="RStoolbox")
metaData <- readMeta(mtlFile)
lsat     <- stackMeta(mtlFile)

lsat_ref <- radCor(lsat, metaData = metaData, method = "apref")

lsat_ref <- crop(lsat_ref,extent(lsat,51,100,51,100))

options(mc.cores = 8)

cl <- MeanShift::msClustering(t(as.matrix(lsat_ref[[3:5]])), multi.core = T)

clr <- setValues(lsat_ref[[1]], cl$labels)

names(clr) <- 'MeanShift'

plot(stack(lsat_ref[[3:5]],clr))