Python tif Kmeans K均值聚类

1120463759

import rasterio as rio
from rasterio.plot import show
from sklearn import cluster
from osgeo import gdal, gdal_array
import matplotlib.pyplot as plt
import matplotlib.colors as mc
import numpy as np
# Open the image
elhas_raster = rio.open("C:/Users/Administrator/Desktop/WKL_Eastern/Li_noNAN_V2/stacking/L58.tif")
print(elhas_raster.meta)

# Read, enhance and show the image
elhas_arr = elhas_raster.read()
vmin, vmax = np.nanpercentile(elhas_arr, (5,95))  # 5-95% contrast stretch
fig, ax = plt.subplots(figsize=[20,20], ncols=1,nrows=1)
show(elhas_raster, cmap='gray', vmin=vmin, vmax=vmax, ax=ax)
ax.set_axis_off()
# fig.savefig("elhas.jpg", bbox_inches='tight')
# plt.show()

# print the shape of the original image
elhas_arr.shape
# create an empty array with same dimension and data type
imgxyb = np.empty((elhas_raster.height, elhas_raster.width, elhas_raster.count), elhas_raster.meta['dtype'])
imgxyb.shape

# loop through the raster's bands to fill the empty array
for band in range(imgxyb.shape[2]):
    imgxyb[:,:,band] = elhas_raster.read(band+1)
print(imgxyb.shape)
# convet to 1d array
img1d = imgxyb[:,:,:58].reshape((imgxyb.shape[0]*imgxyb.shape[1],imgxyb.shape[2]))

img1d.shape
# create an object of the classifier and train it
cl = cluster.KMeans(n_clusters=3)
param = cl.fit(img1d)
cl.labels_
# get the labels of the classes and reshape it x-y-bands shape order (one band only)
img_cl = cl.labels_
img_cl = img_cl.reshape(imgxyb[:,:,0].shape)
img_cl.shape
# Create a custom color map to represent our different 4 classes
#cmap = mc.LinearSegmentedColormap.from_list("", ["black","red","green","yellow"])
# Save result
ds = gdal.Open("C:/Users/Administrator/Desktop/WKL_Eastern/Li_noNAN_V2/stacking/L58.tif")
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
[cols, rows] = arr.shape
format = "GTiff"
driver = gdal.GetDriverByName(format)

outDataRaster = driver.Create("C:/Users/Administrator/Desktop/WKL_Eastern/Li_noNAN_V2/stacking/L58_Kmeans2.tif", rows, cols, 1, gdal.GDT_Byte)
outDataRaster.SetGeoTransform(ds.GetGeoTransform())
outDataRaster.SetProjection(ds.GetProjection())

outDataRaster.GetRasterBand(1).WriteArray(img_cl)
outDataRaster.FlushCache()
del outDataRaster
# Show the resulting array and save it as jpg image
# plt.figure(figsize=[20,20])
# plt.imshow(img_cl, cmap=cmap)
# plt.axis('off')
# plt.savefig("C:/Users/Administrator/Desktop/WKL_Eastern/Li_noNAN_V2/stacking/L09_Kmeans1.tif", bbox_inches='tight')
# plt.show()
print("OK")

1120463759

:mg:mg:mg:mg

zoubo604

haohaoxuexi

XSKYKC2022

膜拜大神，进来学习。

lidonglin

进来学习

zhewin

没事来逛逛

遥望地平线12

感谢分享!

dashuigangxx

感谢楼主分享

dashuigangxx

感谢楼主分享

zhouyangyi

学习学习