基于Python的K-means简单分类

  对于K-means的分类实现，我用的是 jupyter notebook ，这样更方便，可视性更强。
  用Python对卫星数据进行非监督分类，需要 GDAL，Numpy 和 Sklearn，如果查看数据，还需要 Matplotlib ：

import numpy as np
from sklearn import cluster
from osgeo import gdal, gdal_array
import matplotlib.pyplot as plt

# 让GDAL抛出Python异常，并注册所有驱动程序
gdal.UseExceptions()
gdal.AllRegister()

1.单一频段上分类

# 栅格文件读取，使用spot遥感影像
# 获取第一个波段数据
# 将img_ds转化为一个numpy数组，数组的形状为(6000, 6000)
# 可以使用print(img.shape)检查
img_ds = gdal.Open('E:/BaiduNetdiskDownload/spot_pan.tif', gdal.GA_ReadOnly)
band = img_ds.GetRasterBand(1)
img = band.ReadAsArray()

# 将数据展平为行（未知长度），并将列的值保持为1
X = img.reshape((-1,1))

# 对数据运行k-means分类器
# 选择6个分类集群
# 拟合到定义的数据X中
# 给拟合的结果分配一个新变量X_Cluster
# 重新调整原始图像的尺寸
# 此分类过程耗时比较长
k_means = cluster.KMeans(n_clusters=6)
k_means.fit(X)

X_cluster = k_means.labels_
X_cluster = X_cluster.reshape(img.shape)

# 可视化数据 
plt.figure（figsize =（10,10））
plt.imshow（X_cluster，cmap = “ hsv”）
plt.show（）

  分类结果 ：
  4组聚类（效果更好）

  6组聚类：

2. 所有频段上分类

# 栅格文件读取
img_ds = gdal.Open('E:/BaiduNetdiskDownload/spot_pan.tif', gdal.GA_ReadOnly)

# 将多波段图像加载到numpy中(最快的方法)
img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount),
               gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType))

for b in range(img.shape[2]):
    img[:, :, b] = img_ds.GetRasterBand(b + 1).ReadAsArray()

new_shape =(img.shape [0] * img.shape [1],img.shape [2])

# spot有4个波段，将列重塑保持为4
X = img[:, :, :4].reshape(new_shape)

k_means = cluster.KMeans(n_clusters=4)
k_means.fit(X)
X_cluster = k_means.labels_
X_cluster = X_cluster.reshape(img[:, :, 0].shape)

plt.figure(figsize=(10,10))
plt.imshow(X_cluster, cmap="hsv")
plt.show()

  结果如下图所示：

  可见效果，要比单一波段好一些。这段代码的真正好处之一是如何更改分类器。因此，如果在数据加载后使用Mini-Batch K-Means聚类算法，则只需更改一行就可以了。
  K-Means 算法是常用的聚类算法，但其算法本身存在一定的问题，例如在大数据量下的计算时间过长。为此，Mini Batch K-Means，这个基于K-Means的变种聚类算法应运而生。一般当样本量大于1万做聚类时，就需要考虑选用Mini Batch K-Means算法。

MB_KMeans = cluster.MiniBatchKMeans(n_clusters=4)
MB_KMeans.fit(X)

X_cluster = MB_KMeans.labels_
X_cluster = X_cluster.reshape(img[:, :, 0].shape)

  结果如下图所示：

  这是K-Means的更快实现方法，但可能会有更多的噪音。

3.保存分类

  最后需要将分类结果进行保存：

ds = gdal.Open("E:/Download/spot_pan.tif")
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
[cols, rows] = arr.shape

format = "GTiff"
driver = gdal.GetDriverByName(format)


outDataRaster = driver.Create("E:/Download/k_means.gtif", rows, cols, 1, gdal.GDT_Byte)
outDataRaster.SetGeoTransform(ds.GetGeoTransform()) 
outDataRaster.SetProjection(ds.GetProjection()) 


outDataRaster.GetRasterBand(1).WriteArray(X_cluster)

outDataRaster.FlushCache()
del outDataRaster