Python处理栅格数据

Question

想要将以下遥感影像首先进行kmeans聚类，然后对聚类结果进行Delaunay三角网通过约束条件进行二次划分。最后再把结果保存为tif格式影像。

dataset = gdal.Open("F:\Graduate smoothly\MoreBuFen.tif")
rows = dataset.RasterYSize
cols = dataset.RasterXSize
# bands = dataset.RasterCount
# data = dataset.ReadAsArray()
# data_2d = np.transpose(data, (1, 2, 0)).reshape(bands,rows * cols)

band = dataset.GetRasterBand(1)  # 获取栅格数据集的第一个波段
grid_data = band.ReadAsArray()  # 将波段数据读取为一个NumPy数组
data = grid_data.flatten()  # 将二维数组转换为一维数组
data_1d = data.reshape(-1, 1)  # 将数组重新调整形状为一列

kmeans = KMeans(n_clusters=5, n_init='auto').fit(data_1d)
labels = kmeans.labels_.reshape(rows,cols)  # 将聚类标签转换为与影像数据相同的二维数组形式
clusters_partitions = []  # 存储每个簇的空间分区结果
for i in range(kmeans.n_clusters):  # 对每个簇进行处理
    cluster_points = np.argwhere(labels == i)  # 获取第i个簇的所有像素点的坐标
    A = 3
    tri = Delaunay(cluster_points)
    edge_lengths = np.sqrt(
        np.sum((cluster_points[tri.simplices[:, 0]] - cluster_points[tri.simplices[:, 1]]) ** 2,
               axis=1))
    triangles = tri.simplices[edge_lengths <= np.mean(edge_lengths) + A * np.std(edge_lengths)]
    clusters_partitions.append(triangles)  # 存储空间分区结果
partition_result = np.zeros((rows, cols)) # 创建一个与原始影像数据相同大小的数组，用于存储空间分区结果
# # 将每个簇的空间分区结果赋值给partition_result数组
for i, partition in enumerate(clusters_partitions):
     partition_result[np.unravel_index(partition.ravel(), (rows, cols))] = i + 1
# 创建一个与原始影像数据相同的GeoTIFF文件
driver = gdal.GetDriverByName('GTiff')
output_file = driver.Create('F:\Graduate smoothly\Result.tif', cols, rows, 1, gdal.GDT_Float32)
# 将partition_result数组写入GeoTIFF文件
output_file.GetRasterBand(1).WriteArray(partition_result)
# 设置GeoTIFF文件的地理信息
output_file.SetGeoTransform(dataset.GetGeoTransform())
output_file.SetProjection(dataset.GetProjection())
# 保存GeoTIFF文件
output_file.FlushCache()

我按照以上代码得到如下图像，坐标轴代表着行列号，但是是不是数据已经没有了空间关系？应该怎么更改呢？

Answer 1

你修改下代码试着把簇空间分区写入GEOTIFF

import gdal
import numpy as np
from sklearn.cluster import KMeans
from scipy.spatial import Delaunay

dataset = gdal.Open("F:\Graduate smoothly\MoreBuFen.tif")
rows = dataset.RasterYSize
cols = dataset.RasterXSize

geo_transform = dataset.GetGeoTransform()
projection = dataset.GetProjection()
band = dataset.GetRasterBand(1)
grid_data = band.ReadAsArray()

data = grid_data.flatten()
data_1d = data.reshape(-1, 1)

kmeans = KMeans(n_clusters=5, n_init='auto').fit(data_1d)
labels = kmeans.labels_.reshape(rows, cols)

clusters_partitions = []  
for i in range(kmeans.n_clusters):  
    cluster_points = np.argwhere(labels == i)  
    A = 3
    tri = Delaunay(cluster_points)
    edge_lengths = np.sqrt(
        np.sum((cluster_points[tri.simplices[:, 0]] - cluster_points[tri.simplices[:, 1]]) ** 2,
               axis=1))
    triangles = tri.simplices[edge_lengths <= np.mean(edge_lengths) + A * np.std(edge_lengths)]
    clusters_partitions.append(triangles)

partition_result = np.zeros((rows, cols)) 
for i, partition in enumerate(clusters_partitions):
    partition_result[np.unravel_index(partition.ravel(), (rows, cols))] = i + 1

driver = gdal.GetDriverByName('GTiff')
output_file = driver.Create('F:\Graduate smoothly\Result.tif', cols, rows, 1, gdal.GDT_Float32)
output_file.GetRasterBand(1).WriteArray(partition_result)
output_file.SetGeoTransform(geo_transform)
output_file.SetProjection(projection)
output_file.FlushCache()