Texture is one of the most important characteristics dealt with during image interpretation and classification. Texture analysis has been successfully applied to forestry and vegetation studies using a variety of remote sensing data (Asner et al., 2002; Franklin et al., 2000) and radar images (Costa, 2004; Hess et al., 2003). The Gray Level Co-occurrence Matrix (GLCM) is one of the most widely used methods to compute second order texture measures. This method considers the spatial relationships between two pixels at a time which are called the the reference and the neighbor pixel. The neighbor pixel can be to the East (0), North-East (45), North (90), or North-West (135) of the reference pixel. For a more detail, please see here.

Load packages

library (raster)        # raster data
library(rasterVis)    # raster visualisation    
library(sp)               # Spatial data
library(rgdal)        # sptaial data 
library(RStoolbox)    # Image analysis
library(ggplot2)      # advance plotting
library(glcm)         # texture analysis

The data could be available for download from here.

# Define data folder
dataFolder<-"F://Spatial_Data_Processing_and_Analysis_R//Data//DATA_09//"

Load bands

RED<- raster(paste0(dataFolder,"//Landsat8//Landsat8_PROCESSED//_B4_sre.tif"))
NIR<- raster(paste0(dataFolder,"//Landsat8//Landsat8_PROCESSED//_B5_sre.tif"))

Texture Analysis

We will use glcm package to calculate following eight texture statistics of RED and NIR bands derived from grey-level co-occurrence matrices

Red Band Texture

glcm.red <- glcm(RED,
                   window = c(7, 7),
                   shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
                   statistics = c("mean",
                                  "variance",
                                  "homogeneity",
                                  "contrast",
                                  "entropy", 
                                  "dissimilarity",
                                  "second_moment", 
                                  "correlation"))
glcm.red
## class       : RasterStack 
## dimensions  : 504, 541, 272664, 8  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 318195, 334425, 4737735, 4752855  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       :   glcm_mean, glcm_variance, glcm_homogeneity, glcm_contrast, glcm_entropy, glcm_dissimilarity, glcm_second_moment, glcm_correlation 
## min values  :   0.0312500,     0.9255645,        0.1830405,     0.0000000,    0.0000000,          0.0000000,          0.0241566,             -Inf 
## max values  :   0.3770727,   169.8674705,        1.0000000,    71.7908163,    3.7689108,          6.4642857,          1.0000000,              Inf
plot(glcm.red)

     #main= "GLCM Texture of Red Band")

NIR Band texture

glcm.nir <- glcm(NIR,
                   window = c(7, 7),
                   shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)), 
                   statistics = c("mean",
                                  "variance",
                                  "homogeneity",
                                  "contrast",
                                  "entropy", 
                                  "dissimilarity",
                                  "second_moment", 
                                  "correlation"))
plot(glcm.nir)

Write all texture bands

# writeRaster(glcm.red, paste0(dataFolder,  ".//Landsat8//Texture//"),  paste0("RED_", names(glcm.nir)), bylayer=T, overwrite=TRUE, format="GTiff")
# writeRaster(glcm.nir, paste0(dataFolder, ".//Landsat8//Texture//"),  paste0("NIR_", names(glcm.nir)), bylayer=T, overwrite=TRUE, format="GTiff")

Write multiband raster

# writeRaster(glcm.red, filename=paste0(dataFolder,".//Landsat8//Texture//RED_Texture_multi.tif"), bylayer=F, overwrite=TRUE, format="GTiff")
# writeRaster(glcm.nir, filename=paste0(dataFolder,".//Landsat8//Texture//NIR_Texture_multi.tif"), bylayer=F, overwrite=TRUE, format="GTiff")

Principal Componet Analysis of Texture Bands

In order to reduce the redundancy in the texture bands and to determine appropriate texture features, it is neceassry to do principle component analyses (PCA) to the all texture images. We will use rasterPCA() faction of RSToolbox to calculates R-mode PCA for RasterBricks or RasterStacks and returns a RasterBrick with multiple layers of PCA scores.

r<-stack(glcm.red,glcm.nir)  # stack PC raster of red and NIR bands (total raster 8 x 2 = 16)
r
## class       : RasterStack 
## dimensions  : 504, 541, 272664, 16  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 318195, 334425, 4737735, 4752855  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : glcm_mean.1, glcm_variance.1, glcm_homogeneity.1, glcm_contrast.1, glcm_entropy.1, glcm_dissimilarity.1, glcm_second_moment.1, glcm_correlation.1, glcm_mean.2, glcm_variance.2, glcm_homogeneity.2, glcm_contrast.2, glcm_entropy.2, glcm_dissimilarity.2, glcm_second_moment.2, ... 
## min values  :  0.03125000,      0.92556451,         0.18304046,      0.00000000,     0.00000000,           0.00000000,           0.02415660,               -Inf,  0.03140944,      0.92280891,         0.15628695,      0.01020408,     0.04981169,           0.01020408,           0.02186589, ... 
## max values  :   0.3770727,     169.8674705,          1.0000000,      71.7908163,      3.7689108,            6.4642857,            1.0000000,                Inf,   0.8759566,     754.4886370,          0.9948980,     106.0816327,      3.8423098,            7.9693878,            0.9800083, ...

Develop a PCA model

r1<-scale(r)        # scale the data
r1[is.na(r1)] <- 0  # define zero  all miising vaalues
rpc <- rasterPCA(r1, nComp=4)
summary(rpc$model)
## Importance of components:
##                           Comp.1    Comp.2     Comp.3     Comp.4
## Standard deviation     2.2629110 2.2353332 1.14616868 0.85397999
## Proportion of Variance 0.3791344 0.3699498 0.09726471 0.05399501
## Cumulative Proportion  0.3791344 0.7490842 0.84634887 0.90034388
##                            Comp.5     Comp.6     Comp.7      Comp.8
## Standard deviation     0.78446480 0.62713144 0.37113915 0.321962193
## Proportion of Variance 0.04556225 0.02911893 0.01019839 0.007674816
## Cumulative Proportion  0.94590613 0.97502506 0.98522345 0.992898265
##                             Comp.9     Comp.10     Comp.11    Comp.12
## Standard deviation     0.173120352 0.163495250 0.126134487 0.11798550
## Proportion of Variance 0.002218986 0.001979103 0.001177947 0.00103066
## Cumulative Proportion  0.995117251 0.997096354 0.998274302 0.99930496
##                             Comp.13      Comp.14 Comp.15 Comp.16
## Standard deviation     0.0765902470 0.0593417660       0       0
## Proportion of Variance 0.0004343153 0.0002607229       0       0
## Cumulative Proportion  0.9997392771 1.0000000000       1       1

Extact firs 4 PCs from the model

Since, the fist four PCs exlain 90% variabilty of these 16 texture, we will extract these 4 components only.

PC1<-rpc$map$PC1
PC2<-rpc$map$PC2
PC3<-rpc$map$PC3
PC4<-rpc$map$PC4

Write PCA raster

# writeRaster(PC1, filename = paste0(dataFolder,".\\Landsat8\\PCA\\PC1_Texture.tiff"), "GTiff", overwrite=T)
# writeRaster(PC2, filename = paste0(dataFolder,".\\Landsat8\\PCA\\PC2_Texture.tiff"), "GTiff", overwrite=T)
# writeRaster(PC3, filename = paste0(dataFolder,".\\Landsat8\\PCA\\PC3_Texture.tiff"), "GTiff", overwrite=T)
# writeRaster(PC4, filename = paste0(dataFolder,"\\Landsat8\\PCA\\PC4_Texture.tiff"), "GTiff", overwrite=T)

Plot PCs

if(require(gridExtra)){
  plots <- lapply(1:4, function(x) ggR(rpc$map, x, geom_raster = TRUE))
  grid.arrange(plots[[1]],plots[[2]], plots[[3]], plots[[4]],ncol=2)
}

rm(list = ls())