一种非监督分类方法，输入MNDWI或NDWI的TIFF格式文件，得到水体的分类栅格图。要求开放水体面积占比较大，不然效果可能不好。代码简单明了，不做过多说明，转载请注明出处。
try:
from osgeoimport gdal
from osgeoimport ogr
from osgeoimport osr
except ImportError:
import gdal
import ogr
import osr
import numpyas np
import scipy.ioas sio
from skimageimport io
def make_raster(in_ds, fn, data, data_type, nodata=None):
"""Create a one-band GeoTIFF.
in_ds - datasource to copy projection and geotransform from
fn - path to the file to create
data - NumPy array containing data to write
data_type - output data type
nodata - optional NoData value
"""
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
fn, in_ds.RasterXSize, in_ds.RasterYSize, 1, data_type)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)
if nodatais not None:
out_band.SetNoDataValue(nodata)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False)
return out_ds
def Normalized(Onefeature):
L=len(Onefeature.shape);# 判断有无depth
if(L<=2):# 单通道时
minx=np.min(Onefeature); maxX=np.max(Onefeature);
outtempX =255*(Onefeature-minx)/(maxX-minx);
outX = outtempX.astype(np.uint8);# astype:转换数组的数据类型
else:# 多通道时
outX=Onefeature;
for indexin range(Onefeature.shape[2]):#0：行数；1：列数；2：depth数
tempX=Onefeature[:,:,index];#一个通道一个通道的计算
minx = np.min(tempX); maxX = np.max(tempX);
outtempX =255*(tempX-minx)/(maxX-minx);
outX[:,:,index] = outtempX.astype(np.uint8);
return outX;
def OTSU(img_gray):
max_g =0
suitable_th =0
th_begin =0
th_end =256
for thresholdin range(th_begin, th_end):
bin_img = img_gray > threshold
bin_img_inv = img_gray <= threshold
fore_pix = np.sum(bin_img)
back_pix = np.sum(bin_img_inv)
if 0 == fore_pix:
break
if 0 == back_pix:
continue
w0 =float(fore_pix) / img_gray.size
u0 =float(np.sum(img_gray * bin_img)) / fore_pix
w1 =float(back_pix) / img_gray.size
u1 =float(np.sum(img_gray * bin_img_inv)) / back_pix
# intra-class variance
g = w0 * w1 * (u0 - u1) * (u0 - u1)
if g > max_g:
max_g = g
suitable_th = threshold
return suitable_th
def Result2Mat (img_gray,threshold_gray):
Object_class = img_gray > threshold_gray
Back_class = img_gray <= threshold_gray
return Object_class, Back_class
#####----------最大类间方差算法阈值-------------#####
Mutilbands =r'E:\zwcS2MNDWI\GT180408MNDWItiff.tif'
ImageData = io.imread(Mutilbands);
Gray_ImageData = Normalized(ImageData);
suitable_th = OTSU(Gray_ImageData)
min = np.min(ImageData)
max = np.max(ImageData)
over_th = (suitable_th/255)*(max-min) + min
print(over_th);
#####------------保存为.mat文件----------------#####
Object_class, Back_class=Result2Mat(Gray_ImageData,suitable_th)
sio.savemat(r'E:\zwcjieguo\Object_class.mat', {'name': Object_class})
sio.savemat(r'E:\zwcjieguo\Back_class.mat', {'name': Back_class})
#####--------------MAT2tiff-------------------#####
mat = sio.loadmat('E:\zwcjieguo\Object_class.mat')
data=mat['name']
mat_t = np.transpose(data)
in_geofn =r"E:\zwcS2MNDWI\GT180408MNDWItiff.tif"
out_fn =r"E:\zwcjieguo\GT180408Waterbody.tif"
in_geods = gdal.Open(in_geofn)
make_raster(in_geods,out_fn,data,gdal.GDT_Int32,-99)
网友评论