I am trying to develop codes to create GLCM textures images of a geotiff raster image using python. It is working but the resulted texture images are weird and seem to be uncorrect. Could anyone please help me how to fix it?
Thank you
Following are the codes I used
import numpy as np
from numpy import min, max
import rasterio
from rasterio.plot import show
from rasterio.warp import calculate_default_transform, reproject, Resampling
import matplotlib.pyplot as plt
from skimage.feature import graycomatrix, graycoprops
#-----------------------------------------------------------------------------#
def quantize_raster(r, n_levels, method, min_val=None, max_val=None, filename=None, overwrite=False, wopt=None):
og_class = r.__class__.__name__
if og_class == "RasterLayer":
r = r.to_sparraster() # Convert to SpatRaster
elif not isinstance(r, np.ndarray):
raise ValueError("Input must be a RasterLayer or a numpy array.")
if r.ndim != 2:
raise ValueError("Input array must be two-dimensional.")
if method == "equal range":
if min_val is None:
min_val = np.nanmin(r)
if max_val is None:
max_val = np.nanmax(r)
qrules = np.linspace(min_val, max_val, n_levels)
elif method == "equal prob":
qrules = np.quantile(r, q=np.linspace(0, 1, n_levels), interpolation='higher')
else:
raise ValueError("method must be 'equal range' or 'equal prob'")
rq = np.digitize(r, qrules, right=True).astype('uint8')
if og_class == "RasterLayer":
rq = rasterio.open(rq)
if filename is not None:
with rasterio.open(filename, "w", **rq.profile) as dst:
dst.write(rq.read(1), 1)
return rq
if filename is not None:
rq = rasterio.open(filename, "w", driver='GTiff', width=r.shape[1], height=r.shape[0], count=1, dtype='uint8')
rq.write_band(1, rq)
rq.close()
return rq
return rq
#-----------------------------------------------------------------------------#
# Read the image folder
data_name = "D:\Papers\_2023\Seamount_Flores\data\Amp50mUTM.tif"
# Open the geotiff using rasterio
r = rasterio.open(data_name)
show(r, cmap='gray')
# Convert the geotiff into numpy array
r = r.read(1)
img = quantize_raster(r, n_levels = 32, method = 'equal prob')
#-----------------------------------------------------------------------------#
# Define a function to apply a sliding window to an image and compute GLCM for each window
def compute_glcm_sliding_window(image, window_size, distance, angles):
rows, cols = image.shape
Contrast = np.zeros([rows, cols], dtype=np.uint8)
Dissimilarity = np.zeros([rows, cols], dtype=np.uint8)
Homogeneity = np.zeros([rows, cols], dtype=np.uint8)
ASM = np.zeros([rows, cols], dtype=np.uint8)
Energy = np.zeros([rows, cols], dtype=np.uint8)
Correlation = np.zeros([rows, cols], dtype=np.uint8)
for i in range(rows - window_size+1):
print(i)
for j in range(cols - window_size+1):
window = image[i:i+window_size, j:j+window_size]
glcm = graycomatrix(window, [distance], angles, symmetric=True, normed=True)
Contrast[i,j] = graycoprops(glcm, 'contrast')
Dissimilarity[i,j] = graycoprops(glcm, 'dissimilarity')
Homogeneity[i,j] = graycoprops(glcm, 'homogeneity')
ASM[i,j] = graycoprops(glcm, 'ASM')
Energy[i,j] = graycoprops(glcm, 'energy')
Correlation[i,j] = graycoprops(glcm, 'correlation')
return Contrast, Dissimilarity, Homogeneity, ASM, Energy, Correlation
# Define the size of the sliding window
window_size = 3
# Define the offset/distance between pixels for computing the co-occurrence matrix
distance = 1
# Define the angle/direction for computing the co-occurrence matrix
#angles = [0, np.pi/4, np.pi/2, 3*np.pi/4]
angles = [0]
# Compute GLCM for each sliding window in the input image
Contrast, Dissimilarity, Homogeneity, ASM, Energy, Correlation = compute_glcm_sliding_window(img, window_size, distance, angles)
--------------------------------------------------#
from osgeo import gdal, osr, gdal_array
# Open the GeoTIFF file
dataset = gdal.Open("D:\Papers\_2023\Seamount_Flores\data\Amp50mUTM.tif")
# Get the geographic transform parameters
geo_transform = dataset.GetGeoTransform()
"""
The GetGeoTransform() method returns a tuple containing six values in the following order:
1. The X-coordinate of the upper-left corner of the image.
2. The pixel width, which is the size of a pixel in the X direction in geographic coordinates.
3. The rotation, which is usually zero for North-up images.
4. The Y-coordinate of the upper-left corner of the image.
5. The rotation, which is usually zero for North-up images.
6. The pixel height, which is the size of a pixel in the Y direction in geographic coordinates.
"""
# Get the size of the image
width = dataset.RasterXSize
height = dataset.RasterYSize
# The pixel width, which is the size of a pixel in the X direction in geographic coordinates
Grid_x, Grid_y = (geo_transform[1], geo_transform[5])
# Calculate the corner coordinates
min_x = geo_transform[0]
max_y = geo_transform[3]
max_x = geo_transform[0] + (x_size * geo_transform[1])
min_y = geo_transform[3] - (y_size * abs(geo_transform[5]))
upper_left = [geo_transform[0],geo_transform[3]]
upper_right = [(geo_transform[0]+(x_size * geo_transform[1])),geo_transform[3]]
lower_left = [geo_transform[0],geo_transform[3] - (y_size * abs(geo_transform[5]))]
lower_right = [(geo_transform[0]+(x_size * geo_transform[1])), geo_transform[3] - (y_size * abs(geo_transform[5]))]
# Get the projection information
projection = dataset.GetProjection()
# Get the EPSG code from the projection information
#geotransform = (min_x, Grid, 0.0, max_y, 0.0, -Grid)
geotransform = (geo_transform[0], geo_transform[1], 0.0, geo_transform[3], 0.0, geo_transform[5])
srs = osr.SpatialReference(wkt=projection)
srs.ImportFromEPSG(32751)
dst_options = ['COMPRESS=LZW']
# Define the file path and the numpy array
#file_path = "output.tif"
file_path = "D:\Papers\_2023\Seamount_Flores\data\Output_ASM_OK.tif"
# Defining the output Geotiff image
dst_nbands=1
dst_format = 'Gtiff'
#dst_datatype = gdal.GDT_Byte
dst_datatype = gdal.GDT_Float64
dst_options = ['COMPRESS=LZW']
dst_file = 'Seamount_flores.tif'
# retrieves a geotiff driver then uses the driver to create a new raster dataset
driver = gdal.GetDriverByName(dst_format)
dst_ds = driver.Create(dst_file, ydim, xdim, dst_nbands, dst_datatype, dst_options)
# Define spatial and projection reference for the output raster
dst_ds.SetGeoTransform(geotransform) # specify coords
srs = osr.SpatialReference() # establish encoding
srs.ImportFromEPSG(32751)
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
# write to a raster band with a numpy array that matches the raster’s dimensions.
dst_ds.GetRasterBand(1).WriteArray(zi)
# write to disk
dst_ds.FlushCache()
dst_ds = None
Could anyone please help me how to fix it?
