Fuzzy clustering with rasters

Jeremy Gelb

02/04/2021

Source: vignettes/web_vignettes/rasters.Rmd

Introduction

In this vignette, we want to show how to use geocmeans to do fuzzy clustering on a raster dataset. The idea is simple, each pixel of a raster is an observation, and each band of the raster is a variable. Based on this, the same methodology can be applied, but rasters can not be handled like classical tabular and vector dataset. The example given use a medium size raster so some part of the vignette can be quite long to run (in particular when select_parameters is used).

We give a practical example here with a raster of 5 bands:

  • band 1: blue (wavelength: 0.45-0.51)
  • band 2: green (0.53-0.59)
  • band 3: red (0.64-0.67)
  • band 4: near infrared (0.85-0.88)
  • band 5: shortwave infrared (1.57-1.65)
  • band 6: shortwave infrared (2.11-2.29)

The image has a resolution of 30x30 metres and was acquired by Landsat 8 on June 14, 2021. It is provided with geocmeans as a Large RasterBrick Arcachon.

Loading the dataset

Let us start with a simple visualization of the dataset

library(geocmeans)
library(ggpubr)
library(future)
library(tmap)
library(viridis)
library(RColorBrewer)
library(terra)

Arcachon <- terra::rast(system.file("extdata/Littoral4_2154.tif", package = "geocmeans"))
names(Arcachon) <- c("blue", "green", "red", "infrared", "SWIR1", "SWIR2")

# show the pseudo-color image
terra::plotRGB(Arcachon, r = 3, g = 2, b = 1, stretch = "hist")

We can distinguish urban areas, crops, deep water, coast and sand areas. We can expect approximately five clusters.

Fuzzy C-means

As a first exploration and to determine the appropriate number of clusters, we propose the applcication of a classical FCM.

The dataset must be structured as a list of bands. The names given to the list elements will be used as variable names.

# converting the RasterBrick to a simple list of SpatRaster
dataset <- lapply(names(Arcachon), function(n){
  aband <- Arcachon[[n]]
  return(aband)
})

# giving a name to each band
names(dataset) <- names(Arcachon)
# finding an appropriate k and m values (using a multicore plan)
future::plan(future::multisession(workers = 2))
FCMvalues <- select_parameters.mc(algo = "FCM", data = dataset, 
                               k = 5:10, m = seq(1.1,2,0.1), spconsist = FALSE, 
                               indices = c("XieBeni.index", "Explained.inertia",
                                           "Negentropy.index", "Silhouette.index"),
                               verbose = TRUE)
# plotting the silhouette index values
ggplot(FCMvalues) + 
  geom_raster(aes(x = m, y = k, fill = Silhouette.index)) + 
  geom_text(aes(x = m, y = k, label = round(Silhouette.index,2)), size = 2)+
  scale_fill_viridis() +
  coord_fixed(ratio=0.125) 

# plotting the explained inertia
ggplot(FCMvalues) + 
  geom_raster(aes(x = m, y = k, fill = Explained.inertia)) + 
  geom_text(aes(x = m, y = k, label = round(Explained.inertia,2)), size = 2)+
  scale_fill_viridis() +
  coord_fixed(ratio=0.125)

Considering the values above, we select k = 7 and m = 1.5, it seems to provide a good compromise between the silhouette index and the explained inertia.

FCM_result <- CMeans(dataset, k = 7, m = 1.5, standardize = TRUE,
                     verbose = FALSE, seed = 789, tol = 0.001, init = "kpp")

maps1 <- mapClusters(object = FCM_result, undecided = 0.45)

# plotting membership values for group 2
maps1$ProbaMaps[[2]] + theme(legend.position = "bottom")

# plotting membership values for group 5
maps1$ProbaMaps[[5]] + theme(legend.position = "bottom")

# plotting the most likely categories
maps1$ClusterPlot + theme(legend.position = "bottom") + scale_fill_brewer(palette = "Set2")

We can already determine some groups from these maps:

  1. deep water
  2. sand
  3. wet sand
  4. high density urban areas
  5. low density urban areas
  6. vegetation
  7. shallow water

We could try to obtain a more crispy partition by using the generalized version of FCM. This implies selecting a parameter called beta.

GFCMvalues <- select_parameters.mc(algo = "GFCM", data = dataset,
                             k = 7, m = seq(1.1,2,0.1), beta = seq(0.1,0.9,0.1),
                             spconsist = FALSE, verbose = TRUE, init = "kpp",
                             indices = c("XieBeni.index", "Explained.inertia",
                                         "Negentropy.index", "Silhouette.index"))
# plotting the explained inertia
ggplot(GFCMvalues) + 
  geom_raster(aes(x = m, y = beta, fill = Explained.inertia)) + 
  geom_text(aes(x = m, y = beta, label = round(Explained.inertia,2)), size = 2)+
  scale_fill_viridis() +
  coord_fixed(ratio=1)

# plotting the silhouette index
ggplot(GFCMvalues) + 
  geom_raster(aes(x = m, y = beta, fill = Silhouette.index)) + 
  geom_text(aes(x = m, y = beta, label = round(Silhouette.index,2)), size = 2)+
  scale_fill_viridis() +
  coord_fixed(ratio=1)

It seems that beta = 0.5 gives the best results with m = 1.5, we gain 2 more percentage points of inertia while keeping the same Silhouette index.

GFCM_result <- GCMeans(dataset, k = 7, m = 1.5, beta = 0.5, standardize = TRUE,
                       verbose = FALSE, seed = 789, tol = 0.001)

# reorganizing the groups for an easier comparison
GFCM_result <- groups_matching(FCM_result, GFCM_result)

maps2 <- mapClusters(object = GFCM_result, undecided = 0.45)

# plotting membership values for group 2
maps2$ProbaMaps[[2]] + theme(legend.position = "bottom")

# plotting membership values for group 5
maps2$ProbaMaps[[5]] + theme(legend.position = "bottom")

# plotting the most likely categories
maps2$ClusterPlot + theme(legend.position = "bottom") +
  scale_fill_brewer(palette = "Set2")

We can see that with the same threshold, we obtain far less undecided pixels and only few variations among groups in comparison with the previous results.

Spatial Generalized Fuzzy C-means

The image we obtain comprise a lot of noise due to its rather small resolution. To obtain a more spatially homogeneous clustering, we will use here the SGFCM method.

We must determine the right value of alpha and the size of the window. When working with rasters, a square window is used to define the neighbours around a pixel because using a classical neighbour matrix would not be efficient.

w1 <- matrix(1, nrow = 3, ncol = 3)
w2 <- matrix(1, nrow = 5, ncol = 5)
w3 <- matrix(1, nrow = 7, ncol = 7)
future::plan(future::multisession(workers = 2))
SGFCMvalues <- select_parameters.mc(algo = "SGFCM", data = dataset, k = 7, m = 1.5,
                             beta = 0.5, alpha = seq(0.5,2,0.1),
                             window = list(w1,w2,w3),
                             spconsist = TRUE, nrep = 5, 
                             verbose = TRUE, chunk_size = 4,
                             seed = 456, init = "kpp",
                             indices = c("XieBeni.index", "Explained.inertia",
                                         "Negentropy.index", "Silhouette.index"))
dict <- data.frame(
  w = c(1,2,3),
  window = c("3x3","5x5","7x7")
)

SGFCMvalues$window <- dict$window[match(SGFCMvalues$window,dict$w)]

# showing the silhouette index
ggplot(SGFCMvalues) + 
  geom_raster(aes(x = alpha, y = window, fill = Silhouette.index)) + 
  geom_text(aes(x = alpha, y = window, label = round(Silhouette.index,2)), size = 1.5)+
  scale_fill_viridis() +
  coord_fixed(ratio=0.125)

# showing the explained inertia
ggplot(SGFCMvalues) + 
  geom_raster(aes(x = alpha, y = window, fill = Explained.inertia)) + 
  geom_text(aes(x = alpha, y = window, label = round(Explained.inertia,2)), size = 1.5)+
  scale_fill_viridis() +
  coord_fixed(ratio=0.125)

# showing the spatial inconsistency
ggplot(SGFCMvalues) + 
  geom_raster(aes(x = alpha, y = window, fill = spConsistency)) + 
  geom_text(aes(x = alpha, y = window, label = round(spConsistency,2)), size = 1.5)+
  scale_fill_viridis() +
  coord_fixed(ratio=0.125)

Its seems that alpha = 0.9 and a small window of 3x3 gives good results, but we have an important loss in explained inertia. Note that when larger windows need to be used, it is possible to define distance based circular windows with the function circular_window.

SGFCM_result <- SGFCMeans(dataset, k = 7, m = 1.5, standardize = TRUE,
                       lag_method = "mean",
                       window = w1, alpha = 0.9, beta = 0.5,
                       seed = 789, tol = 0.001, verbose = FALSE, init = "kpp")

# reorganizing the groups for an easier comparison
SGFCM_result <- groups_matching(FCM_result, SGFCM_result)

maps3 <- mapClusters(object = SGFCM_result, undecided = 0.2)

# plotting membership values for group 2
maps3$ProbaMaps[[2]] + theme(legend.position = "bottom")

# plotting membership values for group 5
maps3$ProbaMaps[[5]] + theme(legend.position = "bottom")

# plotting the most likely categories
maps3$ClusterPlot + theme(legend.position = "bottom") + 
  scale_fill_brewer(palette = "Set2")

To map the results of the SGFCM, it is necessary to lower the threshold of the parameter undecided. In comparison with the previous maps, we obtain a much smoother result which can be easier to use for mapping purposes or for delimiting areas. As suggested by the drop in the explained inertia, we lost a lot of details. This is not surprising if we consider that the pixels have a size of 30 metres, and the neighbours can be very different from the center of the window. Moreover, group 4 no longer describes wet sand, but interstitial areas between low and high density urban areas.

cluster_results <- list(FCM_result, GFCM_result, SGFCM_result)

indices <- sapply(cluster_results, function(clust){
  c(calcexplainedInertia(clust$Data, clust$Belongings),
    calcSilhouetteIdx(clust$Data, clust$Belongings),
    spConsistency(clust, window = w1, nrep = 5)$Mean)
})

colnames(indices) <- c("FCM", "GFCM", "SGFCM")
rownames(indices) <- c("explained inertia", "silhouette index", "spatial inconsistency")

knitr::kable(indices, digits = 3)
FCM GFCM SGFCM
explained inertia 0.894 0.913 0.894
silhouette index 0.736 0.719 0.704
spatial inconsistency 0.416 0.423 0.287

SGFCM proposes a solution with a spatial inconsistency index divided by three. However, in that context, the spatial algorithm might not be the best choice considering the hard drop of both explained inertia and Silhouette index. Here FCM and FGCM produces very similar results.

Spatial diagnostic

geocmeans provides multiple tools to assess the spatial consistency of a classification. For a raster dataset, the function spatialDiag will calculate the Moran I values for each column in the membership matrix, the spatial inconsistency index and the local ELSA statistic.

diagGFCM <- spatialDiag(GFCM_result, window = matrix(1, nrow = 3, ncol = 3), nrep = 5)

The global Moran I indicates the level of spatial autocorrelation for each category.

knitr::kable(diagGFCM$MoranValues, digits = 3, row.names = FALSE)
Cluster MoranI
Cluster_1 0.959
Cluster_2 0.827
Cluster_3 0.712
Cluster_4 0.654
Cluster_5 0.653
Cluster_6 0.878
Cluster_7 0.888

Here groups 4 and 5 have the lowest levels of spatial autocorrelation, meaning that these groups tend to be less surrounded by pixels from the same groups and are more fragmented. Both groups are urban areas (high and low density) and are spatially close; this is why they have lower Moran I values. We could explore the local Moran I values for both with the function calc_local_moran_raster. This function provides the same results as MoranLocal from raster but is significantly faster.

library(classInt)

# calculating the local Moran I values
loc_moran4 <- calc_local_moran_raster(GFCM_result$rasters$group4,w1)
loc_moran5 <- calc_local_moran_raster(GFCM_result$rasters$group5,w1)

# mapping the values
cols <- rev(RColorBrewer::brewer.pal(n = 8, "Spectral"))
vals <- terra::values(loc_moran4, mat = FALSE)
limits <- classIntervals(vals,  n = 8, style = "kmeans") 
plot(loc_moran4, col = cols, breaks = limits$brks)

vals <- terra::values(loc_moran5, mat = FALSE)
limits <- classIntervals(vals,  n = 8, style = "kmeans") 
plot(loc_moran5, col = cols, breaks = limits$brks)

These maps can be used to investigate the limits of the two categories and determine if the classification is meaningful. Group 7 is by far the most spatially dispersed, which is consistent with its nature: low-density urban areas.

The local Moran I values are useful but are limited to one dimension of the fuzzy classification. An interesting alternative is the ELSA statistic (Naimi et al. 2019). It measures local differences in final group attributions with an entropy based index and combines it with the Euclidean distance between the groups. In other words, it is a measure of spatial autocorrelation for a categorical variable which takes into account the level of dissimilarity between categories. It can be calculated with the following formula:

\[\begin{array}{l} E_{i}=E_{a i} \times E_{c i} \\ E_{a i}=\frac{\sum_{j} w_{i j} d_{i j}}{\max \{d\} \sum_{j} w_{i j}}, j \neq i \\ E_{c i}=-\frac{\sum_{k=1}^{m_{w}} p_{k} \log _{2}\left(p_{k}\right)}{\log _{2} m_{i}}, j \neq i \\ m_{i}=\left\{\begin{array}{l} m \text { if } \sum_{j} w_{i j}>m \\ \sum_{j} w_{i j}, \text { otherwise } \end{array}\right. \\ d_{i j}=\left|c_{i}-c_{j}\right| \end{array}\]

  • \(d_{ij}\) is the dissimilarity (e.g. Euclidean distance) between categories i and j. In a classification context, it can be calculated as the distance between the centres of the clusters (categories).
  • \(p_k\) is the probability of seeing group k within the neighbouring binary window \(w_{ij}\).
cols <- RColorBrewer::brewer.pal(n = 7, "Greys")
vals <- terra::values(diagGFCM$Elsa, mat = FALSE)
limits <- classIntervals(vals[!is.na(vals)],  n = 7, style = "kmeans") 
plot(diagGFCM$Elsa, col = cols, breaks = limits$brks)

This map clearly draws the contours of each cluster. The darker pixels are located at the frontier of clusters or in areas with high heterogeneity.

The original version of ELSA uses a hard partition and thus does not exploit all the information we might get from a fuzzy classification. We propose here an adjusted version: the fuzzy ELSA. As an input, we do not use a categorical vector but a matrix of membership values.

\[\begin{array}{l} \xi_{ij} = (|x_i-x_j|\otimes|x_i-x_j|)\\ Ea_i = \frac{\sum_j(w_{ij}\frac{grandsum(\xi_{ij}\odot d)}{2})}{max \{d\}*\sum_{j} w_{i j}} \\ E_{c i}=-\frac{\sum_{k=1}^{m_{w}} p_{k} \log _{2}\left(p_{k}\right)}{\log _{2} m_{i}}, j \neq i \\ m_{i}=\left\{\begin{array}{l} m \text { if } \sum_{j} w_{i j}>m \\ \sum_{j} w_{i j}, \text { otherwise } \end{array}\right. \\ d_{i j}=\left|c_{i}-c_{j}\right| \end{array}\]

with:

  • x the membership matrix, and \(x_i\) the membership values of observation i to all clusters,
  • \(p_k\) the sum of membership values to each cluster arround observation i including i,
  • \(\otimes\) the outer product between two vectors,
  • \(w\) a binary spatial weight matrix indicating neighbours for each observation,
  • d the dissimilarity matrix between clusters; its diagonal must be filled with zeros,
  • \(\odot\) the element wise product between two matrices,

The fuzzy ELSA can be seen as a generalization of the categorical ELSA and provides the same results if a binary matrix is given. This indicator has not been presented in a paper yet (we are working on it), so it is currently not calculated by spatialDiag. It can be obtained with the function calcFuzzyELSA.

fuzzy_elsa_rast <- calcFuzzyELSA(GFCM_result,window = matrix(1,nrow = 3, ncol = 3))


cols <- RColorBrewer::brewer.pal(n = 7, "Greys")
vals <- terra::values(fuzzy_elsa_rast, mat = FALSE)
limits <- classIntervals(vals[!is.na(vals)],  n = 7, style = "kmeans") 
plot(fuzzy_elsa_rast, col = cols, breaks = limits$brks)

In this case, we can see that the fuzzy ELSA draws more important limits along the coasts, and that the differences between deep and shallow water are weaker. We also have strong artefacts at the South-East (black dots), which are in reality long buildings in a industrial and commercial area.

References

Naimi, Babak, Nicholas AS Hamm, Thomas A Groen, Andrew K Skidmore, Albertus G Toxopeus, and Sara Alibakhshi. 2019. “ELSA: Entropy-Based Local Indicator of Spatial Association.” Spatial Statistics 29: 66–88. https://doi.org/10.1016/j.spasta.2018.10.001.