Are there alternatives to this construction of a minimum spanning tree of a raster for the SKATER algorithm?

Question

This could be seen as a continuation of this question where the top accepted answer uses the spedp library in R in order to generate contiguous clusters of a raster image with multivariate data. In the code they provided (also below, with example dataset here, a raster un UTM Zone 35S containing the results of a supervised classification in 11 classes based on LANDSAT 30m imagery), as far as I understand it, a neighbor list is created with the dimensions of the raster, and then weights are calculated based on the euclidean distance and used to generate a Minimum Spanning Tree (MST) which then runs the SKATER algorithm in order to prune this tree based on feature distance. The user, unfortunately, did not provide a working example, nor the results of the code.

This is used effectively in tutorials for this package with polygons, but the problem with applying this to a raster seems to be that the MST for a raster, at least calculated with the methods available in spdep, will always provide a link where the next connected cell is the one below it until it hits an edge, so any possible pruning of the MST will result in your clusters being in a straight line rather than in a more natural shape.

Finally to the question: Is there a way to circumvent this? Or is SKATER and its associated methods simply not suited for clustering raster data?

# 1. Load packages
packs <- list("tidyverse", "raster", "spdep", "parallel")
lapply(packs, require, character.only = T)

# 2. Load raster layers using raster()
layers <- brick("sample.img")

r1 <- layers[[1]]
r2 <- layers[[2]]
r3 <- layers[[3]]
# And so on...

# 3. Define the number of regions (k + 1)
k <- 10

# 4. Merge explanatory variables in data frame
dat <- lapply(list(r1, r2, r3), values) %>%
   do.call(cbind, .) %>%
   as.data.frame(., stringsAsFactors = F) %>%
   magrittr::set_colnames(c("1984", "1993", "1997"))

# 5. Set up parallel framework for faster computation (optional)
# For a non-parallel, single core execution skip steps 5 and 13
ncores <- detectCores() - 1
cl <- parallel::makeCluster(ncores, type = "PSOCK")
set.coresOption(ncores)
set.ClusterOption(cl)

# 6. Standardize variables
sdat <- scale(dat) %>%
   as.data.frame(., stringsAsFactors = F)

# 7. Create neighbor list object
raster_nb <- cell2nb(nrow(r1), ncol(r1), type = "queen")     # you can alternatively set contiguity to "rook"

# 8. Subset cells
# There are various reasons for which you might need to exclude some pixels to avoid errors in subsequent functions
# One example is missing values in your raster layers (e.g. due to water bodies)
complete_pixels <- which(complete.cases(sdat))
raster_nb <- subset.nb(raster_nb, 1:length(raster_nb) %in% complete_pixels)
sdat <- sdat[complete_pixels,]

# 9. Calculate dissimilarity between neighboring cells
lcosts <- nbcosts(raster_nb, sdat)

# 10. Calculate spatial weights based on dissimilarity between neighbors
raster_w <- nb2listw(raster_nb, lcosts, style = "B")

# 11. Obtain minimum spanning tree
raster_mst <- mstree(raster_w)

# 12. Run skater clustering algorithm
skater_clusters <- skater(raster_mst[,1:2], sdat, k)

# 13. Close parallel framework
stopCluster(cl)

# 14. Visualize results
# Initialize an empty raster with the same dimensions
cluster_raster <- raster(layers)
cluster_raster[] <- NA

# Assign cluster IDs from SKATER results
cluster_raster[complete_pixels] <- skater_clusters$groups

# Visualize
plot(cluster_raster, main="Cluster Assignments")

Answer 1

Your code, as posted above, runs unmodified and produces the following plot, here alongside rasters r1 r2 and r3 so we can tell if there's any relation to the input:

Looks reasonable to me - you can see some of the clusters relating to the input, and it looks nothing like your stripey image. Maybe its a problem with that particular input data.

The only change I had to make to your code was to force the number of cores to be an integer with ncores <- as.integer(detectCores() - 1)

Package versions are:

> packageVersion("spdep")
[1] ‘1.2.8’
> packageVersion("raster")
[1] ‘3.6.23’

And you don't need the tidyverse package bundle installed, code runs fine without it.

I notice that your stripy plot is not at the same location as the sample.img you shared, so I can't be sure this is the same input data that produced your stripy plot. Maybe its a problem with that particular input data.