| import numpy as np |
| from sklearn import cluster |
| from osgeo import gdal, gdal_array |
| import matplotlib.pyplot as plt |
| |
| gdal.UseExceptions() |
| gdal.AllRegister() |
| img_ds = gdal.Open('./raster/LC8_2020.tif', gdal.GA_ReadOnly) |
| |
| band = img_ds.GetRasterBand(2) |
| |
| img = band.ReadAsArray() |
| k_means = cluster.KMeans(n_clusters=8) |
| k_means.fit(X) |
| |
| X_cluster = k_means.labels_ |
| X_cluster = X_cluster.reshape(img.shape) |
| X_cluster.shape |
| plt.figure(figsize=(20,20)) |
| plt.imshow(X_cluster, cmap="hsv") |
| plt.show() |

加入多个波段
| img_ds = gdal.Open('./raster/LC8_2020.tif', gdal.GA_ReadOnly) |
| driver = gdal.GetDriverByName("GTiff") |
| |
| img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount), |
| gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType)) |
| [cols, rows] = img[:, :, 0].shape |
| trans = img_ds.GetGeoTransform() |
| proj = img_ds.GetProjection() |
| |
| 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]) |
| X = img[:, :, :3].reshape(new_shape) |
| (104.769919, |
| 0.0002827044278131621, |
| 0.0, |
| 24.076882, |
| 0.0, |
| -0.0002827880434782616) |
| k_means = cluster.KMeans(n_clusters=8) |
| k_means.fit(X) |
| |
| X_cluster = k_means.labels_ |
| X_cluster = X_cluster.reshape(img[:, :, 0].shape) |
| plt.figure(figsize=(20,20)) |
| plt.imshow(X_cluster, cmap="hsv") |
| |
| plt.show() |

| MB_KMeans = cluster.MiniBatchKMeans(n_clusters=8) |
| MB_KMeans.fit(X) |
| X_cluster = MB_KMeans.labels_ |
| X_cluster = X_cluster.reshape(img[:, :, 0].shape) |
| plt.figure(figsize=(20,20)) |
| plt.imshow(X_cluster, cmap="hsv") |
| |
| plt.show() |

保存至硬盘
| out_img = driver.Create("E:/Desktop/test.tif", rows, cols, 1, gdal.GDT_Byte) |
| out_img.SetGeoTransform(trans) |
| out_img.SetProjection(proj) |
| out_img.GetRasterBand(1).WriteArray(X_cluster) |
| out_img.FlushCache() |
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 推荐几款开源且免费的 .NET MAUI 组件库
· 实操Deepseek接入个人知识库
· 易语言 —— 开山篇
· Trae初体验