Introduction

This notebook explores image segmentation for land cover classification. It aims to help Geomatica Basica students at Universidad Nacional de Colombia to become familiar with geospatial functionalities provided by R base and contributed libraries.

Thus, our objective here is to develop a land classification of a given region based on an satellite image. We will try to do this based on a RGB natural color photograph, that is, an image with three bands representing the visible spectrum. In a natural color image, radiation in the green wavelengths is displayed as green, radiation in the red wavelengths is displayed as red, and radiation in the blue wavelengths is displayed as blue. Usually, these images don’t carry as much information as color infrared images. Thus, it is a real challenge to attempt this task.

The segmentation section in this notebook, both contents and code, follows closely the Richard Plant’s post at r-bloggers. Thanks!

Preliminaries

Le’ts install and load the required libraries:

# run from the console
#install.packages("rasterVis")
#install.packages('OpenImageR')
#install.packages('RImagePalette')
## Loading required package: sp
## Loading required package: lattice
## Loading required package: latticeExtra

1. Reading the image to analyze

The image has been downloaded from the NASA Earth Observatory website. It illustrates the diverse agricultural landscape in the western part of Minas Gerais state in Brazil.

(crops <- stack("../pngs/crops.jpg"))
## class      : RasterStack 
## dimensions : 960, 1440, 1382400, 3  (nrow, ncol, ncell, nlayers)
## resolution : 1, 1  (x, y)
## extent     : 0, 1440, 0, 960  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      : crops.1, crops.2, crops.3 
## min values :       0,       0,       0 
## max values :     255,     255,     255

The fields in this image are located southwest of the city of Perdizes. A mix of regularly-gridded polygonal fields and circular center-pivot fields marks the human use of the region. Small streams of the Araguari River extend like fingers throughout the landscape.

plotRGB(crops, 1, 2, 3)

The photography shows a variety of crops: sunflowers, wheat, potatoes, coffee, rice, soybeans, and corn. Fallow fields —not in active agricultural use— display the violet, reddish, and light tan soils common to this part of Brazil. Darker soils are often rich in iron and aluminum oxides, and are typical of highly weathered soil that forms in hot, humid climates.

True Color Vegetation Indices

A Vegetation Index (VI) is a spectral transformation of two or more bands designed to enhance the contribution of vegetation properties and allow reliable spatial and temporal inter-comparisons of terrestrial photosynthetic activity and canopy structural variations.

There are many Vegetation Indices (VIs), with many being functionally equivalent. Many of the indices make use of the inverse relationship between red and near-infrared reflectance associated with healthy green vegetation. But, what about the potential use of indices combining spectral bands from the visible bands?

Jian et al. (2018) proposed TCVI, a new RGB natural color vegetation index, and compared its performance to the other RGB-based indices.

library(knitr)
include_graphics("../pngs/VIs.png")

From the above table, we will test here two indices: the ExG, and the NGRDI.

This is function to calculate ExG:

ExG <- function(img, i, j, k) {
  bi <- img[[i]]
  bj <- img[[j]]
  bk <- img[[k]]
  vi <- (2*bj-bi-bk) / (bj + bk +bi)
  return(vi)
}

This is function to calculate NGRDI:

NGRDI <- function(img, i, j, k) {
  bi <- img[[i]]
  bj <- img[[j]]
  bk <- img[[k]]
  vi <- (bj-bi) / (bj + bi)
  return(vi)
}

Thanks to Matthias Forkel we can use a nice palette for vegetation indices:

colores <- function(n) {
    .fun <- colorRampPalette(c("chocolate4", "orange", "yellow", "grey", "green", "green3", "darkgreen"))
    col <- .fun(n)
    return(col)
}

Let’s calculate the ExG index:

# For the crops image red = 1, green=2, blue = 1.
(egvi <- ExG(crops, 1, 2,3))
## class      : RasterLayer 
## dimensions : 960, 1440, 1382400  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 0, 1440, 0, 960  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : -0.3543046, 1.33945  (min, max)

Visualization time:

cortes <- seq(-0.35, 1, by=0.1)
cols <- colores(length(cortes)-1)
##plot(ndvitrend, 2, col=cols, breaks=classbreaks)
plot(egvi, 2, col = cols, breaks=cortes, main = "ExG Vegetation Index")

Let’s calculate the NGRDI index:

# For the crops image red = 1, green=2, blue = 1.
(ngvi <- NGRDI(crops, 1, 2,3))
## class      : RasterLayer 
## dimensions : 960, 1440, 1382400  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 0, 1440, 0, 960  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : -0.3809524, 1  (min, max)

Visualization time:

cortes <- seq(-0.35, 1.0, by=0.1)
cols <- colores(length(cortes)-1)
##plot(ndvitrend, 2, col=cols, breaks=classbreaks)
plot(ngvi, 2, col = cols, breaks=cortes, main = "NGRDI Vegetation Index")

2. Image segmentation

In digital image processing and computer vision, image segmentation is the process of partitioning a digital image into multiple segments (sets of pixels, also known as image objects or superpixels). The goal of segmentation is to simplify and/or change the representation of an image into something that is more meaningful and easier to analyze. Image segmentation is typically used to locate objects and boundaries (lines, curves, etc.) in images. More precisely, image segmentation is the process of assigning a label to every pixel in an image such that pixels with the same label share certain characteristics.

There are a number of image analysis packages available in R to conduct image segmentation. We will use the OpenImageR package developed by Mouselimis (2018). The image objects created by this package have a very simple and transparent data structure that makes them ideal for data analysis and, in particular, for integration with the raster package.

The OpenImageR package uses simple linear iterative clustering. The SLIC algorithm measures total distance between pixels as the sum of two components, dlab, the distance in the CIELAB color space, and dxy, the distance in pixel (x,y) coordinates. Specifically, the distance Ds is computed as

\(D_{s} = d_{lab} + \frac{m}{S}. d_{xy}\)

where \(S\) is the grid interval of the pixels and \(m\) is a compactness parameter that allows one to emphasize or de-emphasize the spatial compactness of the superpixels. The algorithm begins with a set of K cluster centers regularly arranged in the (x,y) grid and performs K-means clustering of these centers until a predefined convergence measure is attained.

Let’s read the image using the OpenImageR package:

RGB.Color <- readImage("../pngs/crops.jpg")

Let’s use now SLIC for image segmentation. For detailed explanation of the code see the package’s vignette.

Region.slic = superpixels(input_image = RGB.Color, method = "slic", superpixel = 80,
   compactness = 30, return_slic_data = TRUE,
   return_labels = TRUE, write_slic = "",
   verbose = FALSE)
## Warning in interface_superpixels(input_image, method, superpixel, compactness, :
## The input data has values between 0.000000 and 1.000000. The image-data will be
## multiplied by the value: 255!
imageShow(Region.slic$slic_data)

The function imageShow() is part of the OpenImageR package. The OpenImageR function superpixels() creates superpixels, but it does not actually create an image segmentation in the sense that we define at the beginning of this section. This is done by functions in the SuperPixelImageSegmentation package whose discussion is beyond this notebook scope. For more information, please review the package documentation.

3. Data analysis

The simple structure of the objects generated by the function superpixels() makes it an easy matter to use these objects for our purposes. In last section we read the image data into the object RGB.Color for analysis. Let’s take a look at the structure of this object.

str(RGB.Color)
##  num [1:960, 1:1440, 1:3] 0.886 0.855 0.843 0.82 0.749 ...

It is a three-dimensional array. The “height” and “width” of the box are the number of rows and columns respectively in the image and at each value of the height and width the “depth” is a vector whose three components are the red, green, and blue color values of that cell in the image, scaled to [0,1] by dividing the 8-bit RGB values, which range from 0 to 255, by 255. This is the source of the warning message that accompanied the output of the function superpixels(), which said that the values will be multiplied by 255. For our purposes, however, it means that we can easily analyze images consisting of transformations of these values, such as the NGRDI.

Let’s get a matrix with NGRDI values:

ngvi.mat <- matrix(ngvi@data@values,
            nrow = ngvi@nrows,
            ncol = ngvi@ncols, byrow = TRUE)

The function imageShow() works with data that are either in the eight bit 0 – 255 range or in the [0,1] range (i.e., the range of x between and including 0 and 1). It does not, however, work with NGRVI values if these values are negative. Therefore, we will scale NGRDI values to [0,1].

#rescale between 0 and 1
normalize <- function(x) {
  min <- min(x)
  max <- max(x)
  return(1* (x - min) / (max - min))
}
ngvi.mat1 <- normalize(ngvi.mat)
imageShow(ngvi.mat1)

ngvi.data <- RGB.Color
#ngvi.data[,,1] <- ngvi.mat1
#ngvi.data[,,2] <- ngvi.mat1
#ngvi.data[,,3] <- ngvi.mat1

In the next section we will describe how to create an image segmentation of the NGRDI data and how to use cluster analysis to create a land classification.

4. Image classification

Here is an application of the function superpixels() for land cover classification of the NGRDI data.

ngvi.80 = superpixels(input_image = ngvi.data,
   method = "slic", superpixel = 80,
   compactness = 30, return_slic_data = TRUE,
   return_labels = TRUE, write_slic = "",
   verbose = FALSE)
## Warning in interface_superpixels(input_image, method, superpixel, compactness, :
## The input data has values between 0.000000 and 1.000000. The image-data will be
## multiplied by the value: 255!
imageShow(ngvi.80$slic_data)

Here the structure of the object ngvi.80 .

str(ngvi.80)
## List of 2
##  $ slic_data: num [1:960, 1:1440, 1:3] 226 218 215 209 191 169 151 136 130 127 ...
##  $ labels   : num [1:960, 1:1440] 0 0 0 0 0 0 0 0 0 0 ...

It is a list with two elements, slic_data, a three dimensional array of pixel color data (not normalized), and labels, a matrix whose elements correspond to the pixels of the image. The second element’s name hints that it may contain values identifying the superpixel to which each pixel belongs. Let’s see if this is true.

sort(unique(as.vector(ngvi.80$labels)))
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [51] 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

There are 75 unique labels. Although the call to superpixels() specified 80 superpixels, the function generated 76. We can see which pixels have the label value 0 by setting the values of all of the other pixels to [255, 255, 255], which will plot as white.

R0 <- ngvi.80
for (i in 1:nrow(R0$label))
    for (j in 1:ncol(R0$label))
       if (R0$label[i,j] != 0)
          R0$slic_data[i,j,] <- c(255,255,255)
imageShow(R0$slic_data)

The operation isolates the superpixel in the upper right corner of the image, together with the corresponding portion of the boundary. We can easily use this approach to figure out which value of ngvi.80$label corresponds to which superpixel. Now let’s deal with the boundary. A little exploration of the ngvi.80 object suggests that pixels on the boundary have all three components equal to zero. Let’s isolate and plot all such pixels by coloring all other pixels white.

Bdry <- ngvi.80
  for (i in 1:nrow(Bdry$label))
   for (j in 1:ncol(Bdry$label))
     if (!(Bdry$slic_data[i,j,1] == 0 &
     Bdry$slic_data[i,j,2] == 0 &
       Bdry$slic_data[i,j,3] == 0))
       Bdry$slic_data[i,j,] <- c(255,255,255)
Bdry.norm <- NormalizeObject(Bdry$slic_data)
imageShow(Bdry$slic_data)

The above figure shows that we have indeed identified the boundary pixels. Note that the function imageShow() displays these pixels as white with a black edge, rather than pure black.

Having done a preliminary analysis, we can organize our segmentation process into two steps.

The first step will be to replace each of the superpixels generated by the OpenImageR function superpixels() with one in which each pixel has the same value, corresponding to a measure of central tendency (e.g, the mean, median , or mode) of the original superpixel.

The second step will be to use the K-means unsupervised clustering procedure to organize the superpixels from step 1 into a set of clusters and give each cluster a value corresponding to a measure of central tendency of the cluster.

I used the function make.segments() to carry out the segmentation. The first argument of make.segments() is the superpixels object, and the second is the functional measurement of central tendency. Although in this case each of the three colors of the object ngvi.80 have the same values, this may not be true for every application, so the function analyzes each color separately.

The function, written by Richard Plant, is as follows:

# Identify a measure of central tendency of each superpixel
make.segments <- function(x, ftn){
# The argument ftn is any functional measure of central tendency
   z <- x
# For each identified superpixel, compute measure of central tendency
   for (k in unique(as.vector(x$labels))){
# Identify members of the superpixel having the given label
      in.super <- matrix(0, nrow(x$label), ncol(x$label))
      for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (x$label[i,j] == k)
               in.super[i,j] <- 1
#Identify the boundary cells as having all values 0
      on.bound <- matrix(0, nrow(x$label), ncol(x$label))
      for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (in.super[i,j] == 1){
               if (x$slic_data[i,j,1] == 0 & x$slic_data[i,j,2] == 0 
                  & x$slic_data[i,j,3] == 0)
                     on.bound[i,j] <- 1
         }
#Identify the superpixel cells not on the boundary
      sup.data <- matrix(0, nrow(x$label), ncol(x$label))
         for (i in 1:nrow(x$label))
         for (j in 1:ncol(x$label))
            if (in.super[i,j] == 1 & on.bound[i,j] == 0)
               sup.data[i,j] <- 1
# Compute the measure of central tendency of the cells in R, G, B
      for (n in 1:3){
# Create a matrix M of the same size as the matrix of superpixel values
         M <- matrix(0, dim(x$slic_data)[1], dim(x$slic_data)[2]) 
         for (i in 1:nrow(x$label))
            for (j in 1:ncol(x$label))
# Assign to M the values in the superpixel
               if (sup.data[i,j] == 1) M[i,j] <- x$slic_data[i,j,n]
          if (length(M[which(M > 0 & M < 255)]) > 0)
# Compute the measure of central tendency
            ftn.n <- round(ftn(M[which(M > 0 & M < 255)]), 0)
         else
            ftn.n <- 0
         for (i in 1:nrow(x$label))
            for (j in 1:ncol(x$label))
               if (in.super[i,j] == 1) z$slic_data[i,j,n] <- ftn.n
           }
      }
   return(z)
   }

Here is the application of that function to the object ngvi.80 with the second argument set to mean.

We can see the result.

ngvi.means <- make.segments(ngvi.80, mean)
imageShow(ngvi.means$slic_data)

The next step is to develop clusters that represent identifiable land cover types. In a real project the procedure would be to collect a set of ground truth data from the site, but that option is not available to us. Instead we will work with the true color rendition of the astronaut photo.

The land cover can be subdivided using K-means into any number of types. Note that I have not interpreted here the relationship between the K-means clusters and the land cover classes existing in Perdizes. However, in a real application, this task is a must.

set.seed(123)
ngvi.clus <-
 kmeans(as.vector(ngvi.means$slic_data[,,1]), 10)
vege.class <- matrix(ngvi.clus$cluster,
    nrow = ngvi@nrows,
   ncol = ngvi@ncols, byrow = FALSE)
class.ras <- raster(vege.class, xmn = 0,
    xmx = 1440, ymn = 0, ymx = 960, crs = crs(crops))

Next we can use the raster function ratify() to assign descriptive factor levels to the clusters.

class.ras <- ratify(class.ras)
rat.class <- levels(class.ras)[[1]]
rat.class$landcover <- c("LC 10", "LC 09", "LC 08", "LC 07", "LC 06", "LC 05", "LC 04", "LC 03", "LC 02", "LC 01")
#Create a palette of 10 colors
ncrops <-  jpeg::readJPEG("../pngs/crops.jpg")
crops.pal <- image_palette(ncrops, n=10)
levels(class.ras) <- rat.class
levelplot(class.ras, margin=FALSE,
   col.regions= crops.pal[c(5,7,9,4,8,10,2,4,1,6)],
    main = "Land Cover Types in Perdizes")

Have a look at the output and evaluate whether they are a close representation of the existing land cover classes or not.

In any case, it may be that you want to try these functionalities with your own images. Good luck!!!

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.28          RImagePalette_0.1.1 OpenImageR_1.1.6   
## [4] rasterVis_0.47      latticeExtra_0.6-29 lattice_0.20-41    
## [7] raster_3.1-5        sp_1.4-2           
## 
## loaded via a namespace (and not attached):
##  [1] zoo_1.8-8          tidyselect_1.1.0   xfun_0.14          purrr_0.3.4       
##  [5] colorspace_1.4-1   vctrs_0.3.1        htmltools_0.4.0    viridisLite_0.3.0 
##  [9] yaml_2.2.1         rlang_0.4.6        hexbin_1.28.1      later_1.1.0.1     
## [13] pillar_1.4.4       glue_1.4.1         RColorBrewer_1.1-2 jpeg_0.1-8.1      
## [17] lifecycle_0.2.0    stringr_1.4.0      munsell_0.5.0      gtable_0.3.0      
## [21] codetools_0.2-16   evaluate_0.14      fastmap_1.0.1      httpuv_1.5.4      
## [25] parallel_3.6.3     Rcpp_1.0.4.6       xtable_1.8-4       promises_1.1.1    
## [29] scales_1.1.1       mime_0.9           ggplot2_3.3.1      png_0.1-7         
## [33] digest_0.6.25      stringi_1.4.6      tiff_0.1-5         dplyr_0.8.5       
## [37] shiny_1.4.0.2      grid_3.6.3         rgdal_1.4-8        tools_3.6.3       
## [41] magrittr_1.5       tibble_3.0.1       crayon_1.3.4       pkgconfig_2.0.3   
## [45] ellipsis_0.3.1     assertthat_0.2.1   rmarkdown_2.1      R6_2.4.1          
## [49] compiler_3.6.3