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:
snappy_glcm.py (8.7 KB)
Best regards,
Andreas
thank you andreas for this example! I generally understood the python API but often struggled with the right syntax.
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
http://www.fp.ucalgary.ca/mhallbey/descriptive_stats_group.htm
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. 
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 skimage.io 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):
l.append(p)
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))
ax0.imshow(img)
ax1.imshow(mask)
plt.show()
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]):
axs[k].imshow(fimg[...,k])
plt.show()
visualize_glcm(np.random.choice(trainids))
If you are interested to learn basic to advance then join CETPA INFOTECH and get offer to implement your skill on live projects.
Hi,
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)
try:
im = da.from_array(im, chunks=1000) #dask implementation
except:
pass
newKernel = RectBivariateSpline(np.r_[:ny],np.r_[:nx],im)
return newKernel(yy,xx)
def p_me(Z, win):
'''
loop to calculate greycoprops
'''
try:
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)
except:
return (0,0,0,0,0,0)
def read_raster(in_raster):
in_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.
'''
try:
i = int(shap)
return (i,)
except TypeError:
# shape was not a number
pass
try:
t = tuple(shap)
return t
except TypeError:
# shape was not iterable
pass
raise TypeError('shape must be an int, or a tuple of ints')
def sliding_window(a, ws, ss = None, flatten = True):
'''
Source: http://www.johnvinyard.com/blog/?p=268#more-268
Parameters:
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.
Returns
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:
ds.SetProjection(proj.ExportToWkt())
ds.SetGeoTransform(gt)
ss_band = ds.GetRasterBand(1)
ss_band.WriteArray(std)
ss_band.SetNoDataValue(-99)
ss_band.FlushCache()
ss_band.ComputeStatistics(False)
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)
contrast[merge==0]=np.nan
dissimilarity = im_resize(plt_diss,Nx,Ny)
dissimilarity[merge==0]=np.nan
homogeneity = im_resize(plt_homo,Nx,Ny)
homogeneity[merge==0]=np.nan
energy = im_resize(plt_eng,Nx,Ny)
energy[merge==0]=np.nan
correlation = im_resize(plt_corr,Nx,Ny)
correlation[merge==0]=np.nan
ASM = im_resize(plt_ASM,Nx,Ny)
ASM[merge==0]=np.nan
del plt_cont, plt_diss, plt_homo, plt_eng, plt_corr, plt_ASM
del w,Z,ind,Ny,Nx
driverName= 'GTiff'
epsg_code=26949
proj = osr.SpatialReference()
proj.ImportFromEPSG(epsg_code)
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