1

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.

2 Answers 2

2

****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.

2
  • So in this case I would get a raster with new values,recalculated with meanshift, but everything else is the same, right?
    – Ka_Papa
    Commented May 15, 2018 at 7:49
  • @GeorgeNostradamos, please see the modified answer. The components object can be assigned to the raster. Commented May 15, 2018 at 15:54
1

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))

enter image description here

2
  • Thank you for your reply. You helped me cleanup my R code and actually make it better but I feel that I wasn't clear enough. You see, what I want to do is to understand the last part of any object-oriented classification, where every calculation has been done and it comes right down to making those classes geometric objects. That is what wanted to do. No matter if we are talking about watershed or region-growing, I want to know how to create those geometric objects in a shapefile. Is it just rasterToPolygon command or is there something else?
    – Ka_Papa
    Commented May 16, 2018 at 15:07
  • @GeorgeNostradamos Yes, but I would not use this function to obtain objects. Is too slow and results aren't what you need. Look at the plot, there are only two classes. If you convert this to vector the result will be only two multi-part polygons. First, create objects with another software (eCognition, GRASS, OTB, RSGISLib) and add them to R to apply the statistical part
    – aldo_tapia
    Commented May 16, 2018 at 15:15

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.