GLCM textural metrics for remote sensing image analysis

Introduction

Gray Level Co-Occurrence Matrix (GLCM) (Haralick et al. 1973) texture is a powerful image feature for image analysis. The glcm package provides a easy-to-use function to calculate such textural features for RasterLayer objects in R.

This notebook explores how useful could be textural metrics for identification of crop stress from Landsat 5 images.

Texture metrics background

The following table shows commom texture metrics:

Common texture metrics according to Hall-Beyer (2017)

In the above table, Pi,j is the probability of values i and j occurring in adjacent pixels in the original image within the window defining the neighbourhood. i and j are the labels of the columns and rows (respectively) of the GLCM. Because of the construction of the GLCM, i refers to the digital number (DN) value of a target pixel, and j is the DN value of its immediate neighbour (rook’s case). In the GLCM Cor equation, μ is the mean and σ the standard deviation, both as defined by the equations for GLCM Mean and GLCM Var.

Load libraries & read data

Let’s load the glcm and raster libraries

##install.packages(c("glcm", "rgdal", "raster"))
library(glcm)
library(raster)

Now, let’s add two Landsat 5 color compositions over a USA farm from July 13, 2005. These images have been downloaded from this link: NASA Landsat page

green <- stack("./greenmap.JPG")
plot(green)

In the above figure, four of the fields are shown in R-based default color. The bottom left is a bare field. The top left is a spring wheat field. The top right two fields are sugar beets, and the bottom right is a confectionary sunflower field.

What object is green?

rr green

class      : RasterStack 
dimensions : 942, 988, 930696, 3  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 988, 0, 942  (xmin, xmax, ymin, ymax)
crs        : NA 
names      : greenmap.1, greenmap.2, greenmap.3 
min values :          0,          0,          0 
max values :        255,        255,        255 

rr plotRGB(green, 3, 2, 1, stretch=‘lin’, axes=TRUE, main=321)

Let’s try another library for reading the same color composition:

##install.packages("imager")
library(imager)

ngreen <- load.image("./greenmap.JPG")
plot(ngreen, main="CIMG Natural Color Composition")

rr ngreen

Image. Width: 988 pix Height: 942 pix Depth: 1 Colour channels: 3 

Converting the RGB image into a grayscale

(ggreen <- grayscale(ngreen, method = "XYZ", drop = TRUE))

Image. Width: 988 pix Height: 942 pix Depth: 1 Colour channels: 1 

source("conversions.R")

rgreen <- as.raster.cimg(ggreen)

class(rgreen)

[1] "raster"

## plotting RGB converted to grayscale
plot(rgreen)

As there is no direct conversion from Cimg format into RasterLayer format, we have to take a long path: (i) save the Cimg image as jpg; and (ii) read the jpg as raster layer:

imager::save.image(ggreen,"./RGBgrayscale.jpeg")

(rgreen <- raster("./RGBgrayscale.jpeg"))

class      : RasterLayer 
dimensions : 942, 988, 930696  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 988, 0, 942  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /cloud/project/RGBgrayscale.jpeg 
names      : RGBgrayscale 
values     : 0, 255  (min, max)

plot(rgreen, col = gray.colors(10, start = 0.3, end = 0.9, gamma = 2.2, alpha = NULL),
     axes=TRUE, main= "RGB321 as grayscale")

False color composition

Now, let’s load an infrared color composition RGB 4-3-2 (the near infrared, red, and green bands), a composite combination used for vegetation where darker patches show stressed plants.

nred <- load.image("./redmap.JPG")
plot(nred, main="CIMG False Color Composition")

Converting the RGB432 composition into a grayscale image

(gred <- grayscale(nred, method = "XYZ", drop = TRUE))

Image. Width: 992 pix Height: 937 pix Depth: 1 Colour channels: 1 

imager::save.image(gred,"./NIRgrayscale.jpeg")

(rred <- raster("./NIRgrayscale.jpeg"))

class      : RasterLayer 
dimensions : 937, 992, 929504  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 992, 0, 937  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /cloud/project/NIRgrayscale.jpeg 
names      : NIRgrayscale 
values     : 0, 255  (min, max)

Reading the RGB432 composition as a raster stack

Let’s read the same color composition using the raster package:

red <- stack("./redmap.JPG")
plot(red)

The ground truth

truth <- load.image("./Colored.JPG")
plot(truth, main="The ground truth")

The colored image shows just the near infrared band processed so that the yellow, green and teal colors show stressed plants and the blue and purple show healthier plants. Focusing on the sugar beet field, the magenta color shows areas with excess nitrogen and the lighter green zones are sugar beets under stress.

Textural metrics for green band

Now, we can compute textural metrics using the true color composition:

Greenglcm <- glcm(green[[2]], 
              window = c(9,9), 
              shift = c(1,1), 
              statistics = c("mean", "variance", "homogeneity", "contrast", 
                             "dissimilarity", "entropy", "second_moment", "correlation")
              )

Greenglcm

class      : RasterStack 
dimensions : 942, 988, 930696, 8  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 988, 0, 942  (xmin, xmax, ymin, ymax)
crs        : NA 
names      :   glcm_mean, glcm_variance, glcm_homogeneity, glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_second_moment, glcm_correlation 
min values :  0.18750000,   33.75601034,       0.33384539,    0.00000000,         0.00000000,   0.00000000,         0.05380277,             -Inf 
max values :    1.000000,    961.119602,         1.000000,    125.074074,           7.185185,     3.054695,           1.000000,              Inf 

graphics::par(mfrow = c(3, 3))
plot(Greenglcm)

Textural metrics for RGB321 converted to grayscale image

RGBglcm <- glcm(rgreen, 
              window = c(9,9), 
              shift = c(1,1), 
              statistics = c("mean", "variance", "homogeneity", "contrast", 
                             "dissimilarity", "entropy", "second_moment", "correlation")
              )

RGBglcm

class      : RasterStack 
dimensions : 942, 988, 930696, 8  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 988, 0, 942  (xmin, xmax, ymin, ymax)
crs        : NA 
names      :   glcm_mean, glcm_variance, glcm_homogeneity, glcm_contrast, glcm_dissimilarity, glcm_entropy, glcm_second_moment, glcm_correlation 
min values :  0.18036265,   31.35315691,       0.21610246,    0.00000000,         0.00000000,   0.00000000,         0.02301478,             -Inf 
max values :    1.000000,    961.119602,         1.000000,    117.111111,           7.407407,     3.882937,           1.000000,              Inf 

graphics::par(mfrow = c(3, 3))
plot(RGBglcm)

Textural metrics for infrared band

Now, let’s do the computation using the first band (NIR) of the infrared color composite:

IRglcm <- glcm(red[[1]], 
              window = c(9,9), 
              shift = c(1,1), 
              statistics = c("mean", "variance", "homogeneity", "contrast", 
                             "dissimilarity", "entropy", "second_moment", "correlation")
              )

plot(IRglcm, main="Texture for IR band")

Textural metrics for RGB432 converted to grayscale band

GIRglcm <- glcm(rred, 
              window = c(9,9), 
              shift = c(1,1), 
              statistics = c("mean", "variance", "homogeneity", "contrast", 
                             "dissimilarity", "entropy", "second_moment", "correlation")
              )

plot(GIRglcm, main= "Texture for grayscaled RGB432")

Comparison

A comparison plot for Green GLCM Main vs RGB321 GLCM Main:

graphics::par(mfrow = c(1, 3))
plot(Greenglcm[[1]], main="Main Texture - Green Band")
plot(RGBglcm[[1]], main="Main Texture - RGB321")

plot(truth, main="Ground Truth")

A comparison plot for IR GLCM Main vs RGB432 GLCM Main:

graphics::par(mfrow = c(1, 3))
plot(IRglcm[[1]], main="Main Texture - IR Band")
plot(GIRglcm[[1]], main="Main Texture - RGB432")

plot(truth, main="Ground Truth")

A different visualization

###install.packages("RColorBrewer")
library(RColorBrewer)

my.palette <- brewer.pal(n = 7, name = "Spectral")
graphics::par(mfrow = c(1, 3))
plot(IRglcm[[1]], col = my.palette, main="Main Texture - IR Band")
plot(GIRglcm[[1]], col = my.palette, main="Main Texture - RGB432")

plot(truth, main="Ground Truth")

Conclusion

Multispectral textural metrics computed using grayscaled RGB natural color composites seem to improve discrimination of vegetation stress. However, when using RGB false color composites, multispectral textural metrics don’t seem to provide much information than the one provided by the single IR band.