0

I would like to perform a k-means clustering on time series data. I tried doing it with rts package.

library(rts)
path <- system.file("external", package="rts")
lst <- list.files(path=path,pattern='.asc$',full.names=TRUE)
raster <- stack(lst)

AOI = extent(c(649349,652733,5729977,5732891))
raster=crop(raster,AOI)

dates <- c("2000-02-01","2000-03-01","2000-04-01","2000-05-01")
dates <- as.Date(dates)
timeseries <- rts(raster,dates)

Since this package looks to be built on raster package I tried to retrieve the values like I would have done with raster : values(timeseries) but it doesn't work. I haven't found a way to browse my data in the package documentation.

I also tried to store values in a list :

values = NULL
for (i in 1:ncell(timeseries[[1]]) ){
  values = c(values, timeseries[i])
}

But it stores values as numeric vectors, so I lose the time dimension (even though out of the loop, timeseries[1] gives an xts object).

Do you know a way to do it? Or have you already done a classification or clustering on raster time series with another package or method?

2
  • 1
    You can get the data out of a raster stack easy enough (r[i,j]) so you need to define the "distance" between two time series and then compute a distance matrix for each pixel pair. You can treat the data as a vector, you don't need to make them time series objects (unless you have a package that does time series kmeans) to feed into stats::kmeans?
    – Spacedman
    Jun 23 at 17:02
  • I didn't believe you when you said there is no need to use the data as a time series object. So I found a package designed for time series clustering (TSclust::pam). It turns out it gives the same clusters stats::kmeans gives. You were right, a simple vector works as good as a time series.
    – Anthony
    Jun 24 at 14:39

1 Answer 1

1

Using sample data from the rts package:

path <- system.file("external", package="rts") # location of files
lst <- list.files(path=path,pattern='.asc$',full.names=TRUE)
r <- stack(lst) # creating a RasterStack object

which has four layers, you can feed the as.matrix of that as a four-column matrix to stats::kmeans, and it will (by default) take the root-mean-square (pythagoras in 4d) distance for clustering. Then put the cluster member values into a new single layer raster:

km = kmeans(as.matrix(r), 3)
cl = raster(r)
cl[] = km$cluster
plot(cl)

enter image description here

Showing the cluster membership from 1 to 3. Pixels that are green (cluster 3) come from stacks of 4 pixels that are generally closer together in root-mean-square different than they are to pixels in the white (cluster 1) or yellow (cluster 2) areas according to the kmeans clustering process.

Bonus points: hierarchical clustering. You can use the same idea to feed hclust, but give it a distance matrix, and then cut the tree wherever (3 here to compare with hclust):

hc = hclust(dist(as.matrix(r)))
hcr = raster(r)
hcr[] = cutree(hc,3)
plot(hcr)

enter image description here

Similar but different. How different? You can do a cross-tabulation of cluster assignments for the two methods:

> table(hcr[], cl[])
   
      1   2   3
  1   0 107   0
  2 521  85  94
  3   0   0 483

Endless fun.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.