The sci-kit image package for Python allows to calculate the Grey Level Co-Occurence Matrix (GCLM) and its texture measures as described on this page. However, this description, just as in the previous k-means example, leaves out the process of creating the texture image itself.

On this page you find a very good tutorial explaining GCLM and texture. A texture image is created by moving a window over you raster replacing the centre pixel covered by this window with the calculated texture measure of the moving window, as you see in this image taken from Myrka Hall-Beyers excellent tutorial:

A 7x7 window for any pixel would be defined by

glcm_window = sarraster[i-3: i+4, j-3 : j+4]

And the GCLM for a one-pixel offset [1] to the right [0] for this window as documented here would then be

glcm = greycomatrix(glcm_window, [1], [0], symmetric = True, normed = True )

Now you have to calculate this for every single pixel. The only way for making this specific calculation I can come up with is the slow process of looping through each pixel. Since some software implementations (like the SNAP SAR Toolbox) make this much faster, there must be some other way, so maybe a reader of this post has an idea. (I asked this question also at http://stackoverflow.com/questions/35551249/implementing-glcm-texture-feature-with-scikit-image-and-python ).

The following is only practical for small raster images, named sarraster in this case:

for i in range(sarraster.shape[0] ): print i, for j in range(sarraster.shape[1] ): #windows needs to fit completely in image if i <3 or j <3: continue if i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 4): continue #Define size of moving window glcm_window = sarraster[i-3: i+4, j-3 : j+4] #Calculate GLCM and textures glcm = greycomatrix(glcm_window, [1], [0], symmetric = True, normed = True ) # Calculate texture and write into raster where moving window is centered contrastraster[i,j] = greycoprops(glcm, 'contrast') dissimilarityraster[i,j] = greycoprops(glcm, 'dissimilarity') homogeneityraster[i,j] = greycoprops(glcm, 'homogeneity') energyraster[i,j] = greycoprops(glcm, 'energy') correlationraster[i,j] = greycoprops(glcm, 'correlation') ASMraster[i,j] = greycoprops(glcm, 'ASM') glcm = None glcm_window = None

The following is the original image with various texture images. This is just a first test example, so no clear patterns are seen in this image:

The complete script for creating the image above. For a 1100 x 1400 raster this script takes about an hour while in the SNAP SAR Toolbox such calculations take < 1 minute, so I wonder about alternative implementations:

# -*- coding: utf-8 -*- import matplotlib.pyplot as plt import gdal, gdalconst import numpy as np from skimage.feature import greycomatrix, greycoprops #Read SAR image into Numpy Array filename = "//mnt//glaciology//SARimage.jpg" sarfile = gdal.Open(filename, gdalconst.GA_ReadOnly) sarraster = sarfile.ReadAsArray() #Create rasters to receive texture and define filenames contrastraster = np.copy(sarraster) contrastraster[:] = 0 dissimilarityraster = np.copy(sarraster) dissimilarityraster[:] = 0 homogeneityraster = np.copy(sarraster) homogeneityraster[:] = 0 energyraster = np.copy(sarraster) energyraster[:] = 0 correlationraster = np.copy(sarraster) correlationraster[:] = 0 ASMraster = np.copy(sarraster) ASMraster[:] = 0 # Create figure to receive results fig = plt.figure() fig.suptitle('GLCM Textures') # In first subplot add original SAR image ax = plt.subplot(241) plt.axis('off') ax.set_title('Original Image') plt.imshow(sarraster, cmap = 'gray') for i in range(sarraster.shape[0] ): print i, for j in range(sarraster.shape[1] ): # windows needs to fit completely in image if i > (contrastraster.shape[0] - 4) or j > (contrastraster.shape[0] - 4): continue # Define size of moving window glcm_window = sarraster[i-3: i+4, j-3 : j+4] # Calculate GLCM and textures glcm = greycomatrix(glcm_window, [1], [0], symmetric = True, normed = True ) # Calculate texture and write into raster where moving window is centered contrastraster[i,j] = greycoprops(glcm, 'contrast') dissimilarityraster[i,j] = greycoprops(glcm, 'dissimilarity') homogeneityraster[i,j] = greycoprops(glcm, 'homogeneity') energyraster[i,j] = greycoprops(glcm, 'energy') correlationraster[i,j] = greycoprops(glcm, 'correlation') ASMraster[i,j] = greycoprops(glcm, 'ASM') glcm = None glcm_window = None texturelist = {1: 'contrast', 2: 'dissimilarity', 3: ' homogeneity', 4: 'energy', 5: 'correlation', 6: 'ASM'} for key in texturelist: ax = plt.subplot(2,3,key) plt.axis('off') ax.set_title(texturelist[key]) plt.imshow(eval(texturelist[key] + "raster"), cmap = 'gray') plt.show()