from osgeo import gdal
import cv2 as cv
import numpy as np#定义图像打开方式
def image_open(image):data_array = gdal.Open(image)if data_array == "None":print("数据无法读取")return data_array#定义膨胀特征提取
def dilate(data, result_data):kernel = cv.getStructuringElement(cv.MORPH_RECT, (5, 5))for i in range(0, data.shape[0]):dilate_data = cv.dilate(data[i], kernel)result_data[i] = dilate_datareturn result_data#定义腐蚀特征提取
def erode(data, result_data):kernel = cv.getStructuringElement(cv.MORPH_RECT, (5, 5))for i in range(data.shape[0]):erode_data = cv.erode(data[i], kernel)result_data[i] = erode_datareturn result_data#定义OTSU算法的膨胀和腐蚀优化函数
def dilate_threshold(data):band = data.GetRasterBand(2).ReadAsArray(0, 0, width, height) #提取绿波段band_array = band.astype(np.uint8)threshold, binary = cv.threshold(band_array, 0, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)kernel = cv.getStructuringElement(cv.MORPH_RECT, (5, 5))result = cv.dilate(binary, kernel)return resultdef erode_threshold(data):band = data.GetRasterBand(2).ReadAsArray(0, 0, width, height) #提取绿波段band_array = band.astype(np.uint8)threshold, binary = cv.threshold(band_array, 0, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)kernel = cv.getStructuringElement(cv.MORPH_RECT, (5, 5))result = cv.erode(binary, kernel)return result#定义图像输出函数(面向多波段)
def datasave_MB(Filename, data):output1 = gdal.GetDriverByName("GTiff")output2 = output1.Create(Filename, width, height,data.shape[0], gdal.GDT_Float32)output2.SetProjection(Projection)output2.SetGeoTransform(Transform)for i in range(0, data.shape[0]):output2.GetRasterBand(i+1).WriteArray(data[i])del output2#定义图像输出函数(面向单波段)
def datasave_SB(Filename, data):output1 = gdal.GetDriverByName("GTiff")output2 = output1.Create(Filename, width, height,1, gdal.GDT_Float32)output2.SetProjection(Projection)output2.SetGeoTransform(Transform)output2.GetRasterBand(1).WriteArray(data)del output2#定义相关参数和矩阵
Filepath = r"E:yynctryedata20180911_yync(DA).tif"
data = image_open(Filepath)
width = data.RasterXSize
height = data.RasterYSize
Projection = data.GetProjection()
Transform = data.GetGeoTransform()
data1 = data.ReadAsArray(0, 0, width, height)
result = np.zeros(data1.shape)#定义膨胀和腐蚀的特征提取函数
Filename1 = r"E:yynctryedata膨胀处理.tif"
Filename2 = r"E:yynctryedata腐蚀处理.tif"
Filename3 = r"E:yynctryedataotsu膨胀处理.tif"
Filename4 = r"E:yynctryedataotsu腐蚀处理.tif"#膨胀特征结果输出
dilate_data = dilate(data1, result)
datasave_MB(Filename1, dilate_data)#OTSU算法的膨胀优化结果输出
datasave_SB(Filename3, dilate_threshold(data))#腐蚀特征结果输出
dilate_data = erode(data1, result)
datasave_MB(Filename2, dilate_data)#OTSU算法的腐蚀优化结果输出
datasave_SB(Filename4, erode_threshold(data))

cv.getStructuringElement()