无监督分类：用均值漂移方法对遥感影像进行聚类

python的sklearn模块中有封装好的MeanShift方法，下面试一下看看

# -*- coding: utf-8 -*-
import cv2
import gdal
import numpy as np
from sklearn.cluster import MeanShift, estimate_bandwidth

def read_img(filename):
    dataset=gdal.Open(filename)

    im_width = dataset.RasterXSize
    im_height = dataset.RasterYSize

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0,0,im_width,im_height)

    del dataset 
    return im_proj,im_geotrans,im_width, im_height,im_data


def write_img(filename, im_proj, im_geotrans, im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape 

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)

    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

if __name__ == '__main__':
    img_path = 'E:/sb/sb_test1.tif'
    im_proj, im_geotrans, im_width, im_height, im_data = read_img(img_path)
    im_data = im_data[0:3,...] #遥感影像是4波段的这里取三个波段，不然报错，方法还不支持四个波段
    im_data = im_data.transpose((2,1,0))
    im_temp = im_data.reshape((-1,3))
    im_temp = np.float32(im_temp)
    bandwidth=estimate_bandwidth(im_temp, quantile=0.2, n_samples=500) #自动计算带宽
    ms = MeanShift(bandwidth=bandwidth, bin_seeding=True, cluster_all=True)
    ms.fit_predict(im_temp)
    labels=ms.labels_
    cluster_centers = ms.cluster_centers_
    seg = labels.reshape((im_data.shape[0], im_data.shape[1]))
    seg = seg.transpose(1,0)
    seg_path = 'E:/sb/test/sb_test1_seg.tif'
    write_img(seg_path, im_proj, im_geotrans, seg)

原图：

结果：

数据是10米分辨率的哨兵数据，很明显算法在低分辨率的影像上对水体还有林地这样面积比较大边界不是棱角分明的地物类型是很有效的，下一步的关键就是合并同类。另外，效果还可以通过调整参数quantile=0.2, n_samples=500来调整。

注意这里是自动聚类，本次一共有66个类（参数不同类别不同），不管你视觉上看不看得出哪里是水哪里是林地，你实际是无法确定那个就是水的，类别（属性）并没有标明，只有类别上的区分。