Advanced Techniques With Raster Data: Part 1 – Unsupervised Classification

K-means Algorithm

The k-means clustering algorithm attempts to define the centroid of each cluster with its mean value. This means that data is partitioned into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster. In a geometric interpretation, k-means partitions the data space into Voronoi cells (see plot below). K-means provides an “exclusive” solution, meaning that each observation belongs to one (and only one) cluster. The algorithm is generally efficient for dealing with large data-sets, which commonly happens for raster data (ex. satellite or aerial images.)

CLARA Algorithm

The CLARA (Clustering LARge Application) algorithm is based on the Partition Around Medoids (PAM) algorithm, which in turn is an implementation of K-medoids… … so, let’s try to dig in by parts. Also, since these posts are intended to be ‘short-and-sweet’, I will not use mathematical notation or pseudo-code in the descriptions. There are plenty of awesome books and online resources on these subjects that can you can consult for more information.

The k-medoids algorithm is a clustering algorithm similar to k-means. Both are partitional algorithms in the sense that they break-up the data into groups. They both also attempt to minimize the distance between points labeled in a cluster and a point designated as the center of that cluster.

However, in contrast to the k-means algorithm, k-medoids chooses specific data-points as centers (named as medoids) and works with a generalization of the Manhattan Norm to define the distance between data points. A medoid can be defined as the object of a cluster whose average dissimilarity to all the objects in the cluster is minimal. For example, it is the most centrally located point in the cluster.

Generally, k-medoids are more robust to noise and outliers as compared to k-means because it minimizes a sum of pair-wise dissimilarities instead of a sum of squared Euclidean distances.

PAM is the most common realization of the k-medoids clustering algorithm. PAM uses a greedy search, which may not find the optimum solution, however, it is faster than an exhaustive search. PAM has a high computational cost and uses a large amount of memory to compute the dissimilarity object. This leads us to the CLARA algorithm!
CLARA randomly chooses a small portion of the actual data as a representative of the data; then, medoids are chosen from this sample using PAM. If the sample is robustly selected, in a fairly random manner and with enough data points, it should closely represent the original data-set.

CLARA draws multiple samples of the data-set, applies PAM to each sample, finds the medoids, and then returns its best clustering as the output. At first, a sample data-set is drawn from the original data-set; then the PAM algorithm is applied to find the k medoids. Using these k medoids and the whole data-set, the ‘current’ dissimilarity is calculated. If it is smaller than the one you got in the previous iteration, then these k medoids are kept as the best ones. This process of selection is repeated frequently a specified number of times.

In R, the clara function from the cluster package implements this algorithm. It accepts dissimilarities calculated based on "euclidean" or "manhattan" distance. Euclidean distances are root sum-of-squares of differences, while Manhattan distances are the sum of absolute differences.

Workflow

Our approach to clustering the Landsat 8 spectral raster data will employ two stages:

1. Cluster image data with K-means and CLARA for a number of clusters between 2 and 12
2. Assessing each clustering solutions performance through the average Silhouette Index

Let’s start by downloading and uncompressing the Landsat-8 surface reflectance sample data for the Peneda-Geres National Park:

## Create a folder named data-raw inside the working directory to place downloaded data
if(!dir.exists("./data-raw")) dir.create("./data-raw")

## If you run into download problems try changing: method = "wget"
download.file("https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/LT8_PNPG_MultiBand.zip", "./data-raw/LT8_PNPG_MultiBand.zip", method = "auto")

## Uncompress the zip file
unzip("./data-raw/LT8_PNPG_MultiBand.zip", exdir = "./data-raw")

Now, let’s read the multi-band GeoTIFF file into R as a RasterBrick object:

library(raster)

# Load the multi-band GeoTIFF file with brick function
rst <- brick("./data-raw/LC82040312015193LGN00_sr_b_1_7.tif")

# Change band names      
names(rst) <- paste("b",1:7,sep="")  

Plot the data in the RGB display (bands 4,3,2) to see if everything is fine:

plotRGB(rst, r=4, g=3, b=2, scale=10000, stretch="lin", main="RGB composite (b4,b3,b2) of Landsat-8")

Data is OK, which means we can proceed to step #1. For example, clustering raster data with both algorithms and for different numbers of clusters (from 2 to 12).

For simplifying the processing workflow, the ‘internal’ values of the raster object will be entirely loaded up into memory into a data.frame object (in this case with 2.286610^{6} rows – one for each cell; and, 7 columns – one for each layer) You can also see ?values and ?getValues for more info on this. Although this makes things easier and faster, in some cases, depending on the size of the image and its number of layers/bands, it will be unfeasible to push all the data into RAM.

However, if your raster object is too large to fit into a data frame object in memory, you can still use R to perform K-means clustering. Packages such as RSToolbox provide an implementation of k-means that may be more suited for your case. You can check here and the function unsuperClass here.

Also, we will need to be careful regarding NA’s because clustering algorithms will not run with these values (typically throwing an error). One simple way is to use a logical index to sub-set the data.

Let’s see how this works out in actual R code (use comments as guidance):

library(cluster)

# Extract all values from the raster into a data frame
rstDF <- values(rst)

# Check NA's in the data
idx <- complete.cases(rstDF)

# Initiate the raster datasets that will hold all clustering solutions 
# from 2 groups/clusters up to 12
rstKM <- raster(rst[[1]])
rstCLARA <- raster(rst[[1]])


for(nClust in 2:12){
  
  cat("-> Clustering data for nClust =",nClust,"......")
  
  # Perform K-means clustering
  km <- kmeans(rstDF[idx,], centers = nClust, iter.max = 50)
  
  # Perform CLARA's clustering (using manhattan distance)
  cla <- clara(rstDF[idx, ], k = nClust, metric = "manhattan")
  
  # Create a temporary integer vector for holding cluster numbers
  kmClust <- vector(mode = "integer", length = ncell(rst))
  claClust <- vector(mode = "integer", length = ncell(rst))
  
  # Generate the temporary clustering vector for K-means (keeps track of NA's)
  kmClust[!idx] <- NA
  kmClust[idx] <- km$cluster
  
  # Generate the temporary clustering vector for CLARA (keeps track of NA's too ;-)
  claClust[!idx] <- NA
  claClust[idx] <- cla$clustering
  
  # Create a temporary raster for holding the new clustering solution
  # K-means
  tmpRstKM <- raster(rst[[1]])
  # CLARA
  tmpRstCLARA <- raster(rst[[1]])

  # Set raster values with the cluster vector
  # K-means
  values(tmpRstKM) <- kmClust
  # CLARA
  values(tmpRstCLARA) <- claClust
  
  # Stack the temporary rasters onto the final ones
  if(nClust==2){
    rstKM    <- tmpRstKM
    rstCLARA <- tmpRstCLARA
  }else{
    rstKM    <- stack(rstKM, tmpRstKM)
    rstCLARA <- stack(rstCLARA, tmpRstCLARA)
  }
  
  cat(" done!\n\n")
}

# Write the clustering solutions for each algorithm
writeRaster(rstKM,"./data-raw/LT8_PGNP_KMeans_nc2_12-1.tif", overwrite=TRUE)
writeRaster(rstCLARA,"./data-raw/LT8_PGNP_CLARA_nc2_12-1.tif", overwrite=TRUE)