Invariant GLCM textural metrics for remote sensing image analysis

Introduction

Gray Level Co-Occurrence Matrix (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.

A previous notebook explored how useful are GLCM textural metrics computed in one direccion.

In this notebook, textural metrics for the IR band are calculated in all 4 directions (0°, 45°, 90° and 135°) and then combined to one rotation-invariant texture.

Load libraries & read data

Let’s load the glcm and raster libraries

Loading required package: magrittr

Attaching package: ‘magrittr’

The following object is masked from ‘package:raster’:

    extract


Attaching package: ‘imager’

The following object is masked from ‘package:magrittr’:

    add

The following object is masked from ‘package:raster’:

    bbox

The following object is masked from ‘package:sp’:

    bbox

The following objects are masked from ‘package:stats’:

    convolve, spectrum

The following object is masked from ‘package:graphics’:

    frame

The following object is masked from ‘package:base’:

    save.image

Reading the RGB432 composition as a raster stack

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

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

The ground truth

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

Rotation invariant textural metrics for infrared band

Now, let’s do the computation using the first band (NIR) of the infrared color composite. The key for the rotation invariant is the shift parameter:

rr IRglcm1 <- glcm(red[[1]], window = c(9,9), shift = list(c(0,1), c(1,1), c(1,0), c(1,-1)), statistics = c(, , , , , , _moment, ) )

rr plot(IRglcm1, main=invariant Texture for IR band)

Writing the textural stack

Just in case, the cloud crashes:

rr if (require(rgdal)) { rf <- writeRaster(IRglcm1, filename=file.path(./IRglcm1.tif), format=, overwrite=TRUE) }

Loading required package: rgdal
rgdal: version: 1.4-6, (SVN revision 841)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
 Path to GDAL shared files: /usr/share/gdal/2.2
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.3-1 

Comparison between GLCM Main vs GLCM Variance

A comparison plot for GLCM Main vs GLCM Variance:

rr graphics::par(mfrow = c(1, 3)) plot(IRglcm1[[1]], main=Main Texture) plot(IRglcm1[[2]], main=Variance Texture) r plot(truth, main=Truth)

Principal Components

rr (IRglcm1 <- stack(./IRglcm1.tif))

class      : RasterStack 
dimensions : 937, 992, 929504, 8  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 992, 0, 937  (xmin, xmax, ymin, ymax)
crs        : NA 
names      :   IRglcm1.1,   IRglcm1.2,   IRglcm1.3,   IRglcm1.4,   IRglcm1.5,   IRglcm1.6,   IRglcm1.7,   IRglcm1.8 
min values :  0.25559414, 75.71328307,  0.50627324,  0.00000000,  0.00000000,  0.00000000,  0.05426002, -0.57119603 
max values :    1.000000,  961.038870,    1.000000,  261.438272,    9.472222,    3.041933,    1.000000,    1.435670 

rr ### install.packages(2) library(ggplot2)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

rr (IRglcm2 <- stack(IRglcm1[[1]], IRglcm1[[2]]))

class      : RasterStack 
dimensions : 937, 992, 929504, 2  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 992, 0, 937  (xmin, xmax, ymin, ymax)
crs        : NA 
names      :  IRglcm1.1,  IRglcm1.2 
min values :  0.2555941, 75.7132831 
max values :     1.0000,   961.0389 

rr ## Run PCA set.seed(25) ### rasterPCA(img, nSamples = NULL, nComp = nlayers(img), spca = FALSE, ### maskCheck = TRUE, …) rpc <- rasterPCA(IRglcm2, nSamples= 10000, nComp=3)

rr rpc

$call
rasterPCA(img = IRglcm2, nSamples = 10000, nComp = 3)

$model
Call:
princomp(x = trainData, cor = spca, scores = FALSE)

Standard deviations:
      Comp.1       Comp.2 
185.67904439   0.01802312 

 2  variables and  10000 observations.

$map
class      : RasterBrick 
dimensions : 937, 992, 929504, 2  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 992, 0, 937  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      :          PC1,          PC2 
min values : -484.5749256,   -0.1120525 
max values : 400.75098361,   0.03006738 


attr(,\class\)
[1] \rasterPCA\ \RStoolbox\

rr ## Model parameters: summary(rpc$model)

Importance of components:
                        Comp.1       Comp.2
Standard deviation     185.679 1.802312e-02
Proportion of Variance   1.000 9.421810e-09
Cumulative Proportion    1.000 1.000000e+00

rr loadings(rpc$model)

Loadings:
          Comp.1 Comp.2
IRglcm1.1         1.000
IRglcm1.2  1.000       

               Comp.1 Comp.2
SS loadings       1.0    1.0
Proportion Var    0.5    0.5
Cumulative Var    0.5    1.0

rr rpc$map

class      : RasterBrick 
dimensions : 937, 992, 929504, 2  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 992, 0, 937  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      :          PC1,          PC2 
min values : -484.5749256,   -0.1120525 
max values : 400.75098361,   0.03006738 

rr graphics::par(mfrow = c(1, 3)) plot(rpc\(map[[1]], main=\Texture PC1\) plot(rpc\)map[[2]], main=PC2) r plot(truth, main=Truth)

Main Texture vs Principal Component 1

rr graphics::par(mfrow = c(1, 3)) plot(IRglcm1[[1]], main=Main Texture) plot(rpc$map[[1]], main=PC1) r plot(truth, main=Truth)

Conclusion

Principal component 1 seems to summarize the big share of information carried by textural metrics.