Python codes of GLCM for texture feature extraction

Any idea how to have access to Gray Level Co-occurence matrix (GLCM) python codes for SAR texture feature extraction? I would like to run the texture analysis on SAR Terrain correction data in order to produce “entropy”, but through the python.

Your help is very appreciated.

Good morning SAR2016,

I once wrote a python code for it, but I could not really test it because of “out-of-memory” on my computer. But it should work: (8.7 KB)

Best regards,


thank you andreas for this example! I generally understood the python API but often struggled with the right syntax.

1 Like

Morning, Thanks a lot for sharing it.
My SAR data is 1.4 Gb which has caused us to use a 32 gb ram .

I am going to test it on my data (SAR Terrain Correction geotiff) and will keep you posted about it.

PS- This is a good tutorial for those who might be interested in learning more about how GLCM works

by the way, for the memory issue, this probably helps.

Yes, I changed it already to the maxium available. But the problem is I work with Sentinel-2 (~6 GB). So with my 8GB Desktop I have no chance. Not even with my private 16GB Laptop. :slight_smile:

Can you elaborate the procedure to find the “java max mem” and how to increase it? i am using only snap desktop app and getting java heap space error, not familiar with snappy and commands. so, kindly help me to increase the java heap space.

my laptop is of 12gb ram 64 bit, i am working on sentinel 1 data. If i increase java heap space, will it be helpful?

Thanks in advance

Look at the following links to know how to increase SNAP memory…

Increase snappy memory(python): in python when we face with the problem of
java.lang.OutOfMemoryError: Java heap space or Data Buffer

good link to increase memory when we use GPT:(comment 6)

GLCMs( grey level co-occurrence matrics )s features are good for analyzing images with spatial variations without fixed objectiveness like seismic data. They are obtained by summing up all co-occurrences of grey scale values at a specifed offset (distance and angle in 2d case) over an image, with following aggregations. They further detailed ‘dissimilarity’, ‘contrast’, ‘homogeneity’, ‘energy’ and ‘correlation’ by ways of aggregation. One can google “glcm + seismic” to find applications of GLCMs.

scikit-image, GLCM features API
GLCM Texture Features example code

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import tqdm
from import imread

trainids = pd.read_csv(’…/input/train.csv’)[‘id’].tolist()
def read_image(imgid):
fn = ‘…/input/train/images/{}.png’.format(imgid)
return imread(fn)[…,0].astype(np.float32) / 255

def read_mask(imgid):
fn = ‘…/input/train/masks/{}.png’.format(imgid)
return imread(fn).astype(np.uint8)
from skimage.feature import greycomatrix, greycoprops
from multiprocessing import Pool

def glcm_props(patch):
lf = []
props = [‘dissimilarity’, ‘contrast’, ‘homogeneity’, ‘energy’, ‘correlation’]

# left nearest neighbor
glcm = greycomatrix(patch, [1], [0], 256, symmetric=True, normed=True)
for f in props:
    lf.append( greycoprops(glcm, f)[0,0] )

# upper nearest neighbor
glcm = greycomatrix(patch, [1], [np.pi/2], 256, symmetric=True, normed=True)
for f in props:
    lf.append( greycoprops(glcm, f)[0,0] )
return lf

def patch_gen(img, PAD=4):
img1 = (img * 255).astype(np.uint8)

W = 101
imgx = np.zeros((101+PAD*2, 101+PAD*2), dtype=img1.dtype)
imgx[PAD:W+PAD,PAD:W+PAD] = img1
imgx[:PAD,  PAD:W+PAD] = img1[PAD:0:-1,:]
imgx[-PAD:, PAD:W+PAD] = img1[W-1:-PAD-1:-1,:]
imgx[:, :PAD ] = imgx[:, PAD*2:PAD:-1]
imgx[:, -PAD:] = imgx[:, W+PAD-1:-PAD*2-1:-1]

xx, yy = np.meshgrid(np.arange(0, W), np.arange(0, W))
xx, yy = xx.flatten() + PAD, yy.flatten() + PAD

for x, y in zip(xx, yy):
    patch = imgx[y-PAD:y+PAD+1, x-PAD:x+PAD+1]
    yield patch

def glcm_feature(img, verbose=False):

W, NF, PAD = 101, 10, 4

if img.sum() == 0:
    return np.zeros((W,W,NF), dtype=np.float32)

l = []
with Pool(3) as pool:
    for p in tqdm.tqdm(pool.imap(glcm_props, patch_gen(img, PAD)), total=W*W, disable=not verbose):
fimg = np.array(l, dtype=np.float32).reshape(101, 101, -1)
return fimg

def visualize_glcm(imgid):
img = read_image(imgid)
mask = read_mask(imgid)

fimg = glcm_feature(img, verbose=1)

_, (ax0, ax1) = plt.subplots(1, 2, figsize=(6,3))

amin = np.amin(fimg, axis=(0,1))
amax = np.amax(fimg, axis=(0,1))
fimg = (fimg - amin) / (amax - amin)

fimg[...,4] = np.power(fimg[...,4], 3)
fimg[...,9] = np.power(fimg[...,9], 3)

_, axs = plt.subplots(2, 5, figsize=(15,6))
axs = axs.flatten()

for k in range(fimg.shape[-1]):


If you are interested to learn basic to advance then join CETPA INFOTECH and get offer to implement your skill on live projects.


I also had the same issue with python code and here my solution -

import gdal, osr
import numpy as np
from scipy.interpolate import RectBivariateSpline
from numpy.lib.stride_tricks import as_strided as ast
import dask.array as da
from joblib import Parallel, delayed, cpu_count
import os
from skimage.feature import greycomatrix, greycoprops

def im_resize(im,Nx,Ny):
    resize array by bivariate spline interpolation
    ny, nx = np.shape(im)
    xx = np.linspace(0,nx,Nx)
    yy = np.linspace(0,ny,Ny)

        im = da.from_array(im, chunks=1000)   #dask implementation

    newKernel = RectBivariateSpline(np.r_[:ny],np.r_[:nx],im)
    return newKernel(yy,xx)

def p_me(Z, win):
    loop to calculate greycoprops
        glcm = greycomatrix(Z, [5], [0], 256, symmetric=True, normed=True)
        cont = greycoprops(glcm, 'contrast')
        diss = greycoprops(glcm, 'dissimilarity')
        homo = greycoprops(glcm, 'homogeneity')
        eng = greycoprops(glcm, 'energy')
        corr = greycoprops(glcm, 'correlation')
        ASM = greycoprops(glcm, 'ASM')
        return (cont, diss, homo, eng, corr, ASM)
        return (0,0,0,0,0,0)

def read_raster(in_raster):
    ds = gdal.Open(in_raster)
    data = ds.GetRasterBand(1).ReadAsArray()
    data[data<=0] = np.nan
    gt = ds.GetGeoTransform()
    xres = gt[1]
    yres = gt[5]

    # get the edge coordinates and add half the resolution 
    # to go to center coordinates
    xmin = gt[0] + xres * 0.5
    xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
    ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
    ymax = gt[3] - yres * 0.5
    del ds
    # create a grid of xy coordinates in the original projection
    xx, yy = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
    return data, xx, yy, gt

def norm_shape(shap):
   Normalize numpy array shapes so they're always expressed as a tuple,
   even for one-dimensional shapes.
      i = int(shap)
      return (i,)
   except TypeError:
      # shape was not a number

      t = tuple(shap)
      return t
   except TypeError:
      # shape was not iterable

   raise TypeError('shape must be an int, or a tuple of ints')

def sliding_window(a, ws, ss = None, flatten = True):
        a  - an n-dimensional numpy array
        ws - an int (a is 1D) or tuple (a is 2D or greater) representing the size 
             of each dimension of the window
        ss - an int (a is 1D) or tuple (a is 2D or greater) representing the 
             amount to slide the window in each dimension. If not specified, it
             defaults to ws.
        flatten - if True, all slices are flattened, otherwise, there is an 
                  extra dimension for each dimension of the input.

        an array containing each n-dimensional window from a
    if None is ss:
        # ss was not provided. the windows will not overlap in any direction.
        ss = ws
    ws = norm_shape(ws)
    ss = norm_shape(ss)
    # convert ws, ss, and a.shape to numpy arrays
    ws = np.array(ws)
    ss = np.array(ss)
    shap = np.array(a.shape)
    # ensure that ws, ss, and a.shape all have the same number of dimensions
    ls = [len(shap),len(ws),len(ss)]
    if 1 != len(set(ls)):
        raise ValueError(\
        'a.shape, ws and ss must all have the same length. They were %s' % str(ls))

    # ensure that ws is smaller than a in every dimension
    if np.any(ws > shap):
        raise ValueError(\
        'ws cannot be larger than a in any dimension.\
     a.shape was %s and ws was %s' % (str(a.shape),str(ws)))

    # how many slices will there be in each dimension?
    newshape = norm_shape(((shap - ws) // ss) + 1)

    # the shape of the strided array will be the number of slices in each dimension
    # plus the shape of the window (tuple addition)
    newshape += norm_shape(ws)

    # the strides tuple will be the array's strides multiplied by step size, plus
    # the array's strides (tuple addition)
    newstrides = norm_shape(np.array(a.strides) * ss) + a.strides
    a = ast(a,shape = newshape,strides = newstrides)
    if not flatten:
        return a
    # Collapse strided so that it has one more dimension than the window.  I.e.,
    # the new array is a flat list of slices.
    meat = len(ws) if ws.shape else 0
    firstdim = (np.product(newshape[:-meat]),) if ws.shape else ()
    dim = firstdim + (newshape[-meat:])
    # remove any dimensions with size 1
    dim = filter(lambda i : i != 1,dim)

    return a.reshape(dim), newshape

def CreateRaster(xx,yy,std,gt,proj,driverName,outFile):  
    Exports data to GTiff Raster
    std = np.squeeze(std)
    std[np.isinf(std)] = -99
    driver = gdal.GetDriverByName(driverName)
    rows,cols = np.shape(std)
    ds = driver.Create( outFile, cols, rows, 1, gdal.GDT_Float32)      
    if proj is not None:  
    ss_band = ds.GetRasterBand(1)
    del ds

#Stuff to change

if __name__ == '__main__':  
    win_sizes = [7]
    for win_size in win_sizes[:]:   
        in_raster = #Path to input raster
        win = win_size
        meter = str(win/4)

        #Define output file names
        contFile = 
        dissFile = 
        homoFile = 
        energyFile = 
        corrFile =
        ASMFile = 

        merge, xx, yy, gt = read_raster(in_raster)

        merge[np.isnan(merge)] = 0

        Z,ind = sliding_window(merge,(win,win),(win,win))

        Ny, Nx = np.shape(merge)

        w = Parallel(n_jobs = cpu_count(), verbose=0)(delayed(p_me)(Z[k]) for k in xrange(len(Z)))

        cont = [a[0] for a in w]
        diss = [a[1] for a in w]
        homo = [a[2] for a in w]
        eng  = [a[3] for a in w]
        corr = [a[4] for a in w]
        ASM  = [a[5] for a in w]

        #Reshape to match number of windows
        plt_cont = np.reshape(cont , ( ind[0], ind[1] ) )
        plt_diss = np.reshape(diss , ( ind[0], ind[1] ) )
        plt_homo = np.reshape(homo , ( ind[0], ind[1] ) )
        plt_eng = np.reshape(eng , ( ind[0], ind[1] ) )
        plt_corr = np.reshape(corr , ( ind[0], ind[1] ) )
        plt_ASM =  np.reshape(ASM , ( ind[0], ind[1] ) )
        del cont, diss, homo, eng, corr, ASM

        #Resize Images to receive texture and define filenames
        contrast = im_resize(plt_cont,Nx,Ny)
        dissimilarity = im_resize(plt_diss,Nx,Ny)
        homogeneity = im_resize(plt_homo,Nx,Ny)
        energy = im_resize(plt_eng,Nx,Ny)
        correlation = im_resize(plt_corr,Nx,Ny)
        ASM = im_resize(plt_ASM,Nx,Ny)
        del plt_cont, plt_diss, plt_homo, plt_eng, plt_corr, plt_ASM

        del w,Z,ind,Ny,Nx

        driverName= 'GTiff'    
        proj = osr.SpatialReference()

        CreateRaster(xx, yy, contrast, gt, proj,driverName,contFile) 
        CreateRaster(xx, yy, dissimilarity, gt, proj,driverName,dissFile)
        CreateRaster(xx, yy, homogeneity, gt, proj,driverName,homoFile)
        CreateRaster(xx, yy, energy, gt, proj,driverName,energyFile)
        CreateRaster(xx, yy, correlation, gt, proj,driverName,corrFile)
        CreateRaster(xx, yy, ASM, gt, proj,driverName,ASMFile)

        del contrast, merge, xx, yy, gt, meter, dissimilarity, homogeneity, energy, correlation, ASM