weixin_313155671. 栅格数据分区统计实战指南
分区统计是GIS分析中的基础操作，它能帮助我们快速获取特定区域内的栅格数据统计特征。在实际项目中，我经常使用Python的rasterstats模块来处理这类需求，它的zonal_stats()函数既高效又灵活。
1.1 zonal_stats()函数详解
这个函数的核心参数有三个：
- vector_file：矢量边界文件路径（支持shp、geojson等格式）
- raster_file：待统计的栅格文件路径
- stats：需要计算的统计量列表
重要提示：矢量文件和栅格文件必须使用相同的坐标系，否则统计结果将不准确。我建议在操作前先用QGIS或ArcGIS检查两者的CRS是否一致。
默认情况下，函数会返回四个基础统计量：
- min：区域内像元最小值
- max：区域内像元最大值
- mean：区域内像元平均值
- count：区域内有效像元数量
如果需要更丰富的统计指标，可以通过stats参数指定，例如：
python复制stats=['min','max','mean','median','sum','std','count']
1.2 实战案例：上海市分区统计
下面这个案例展示了如何统计上海各行政区的土地类型数据：
python复制import pandas as pd
from rasterstats import zonal_stats
# 文件路径配置
shp_file = r'C:\data\shanghai_proj.shp'
raster_file = r'C:\data\上海.tif'
# 执行分区统计
stats = zonal_stats(shp_file, raster_file,
stats=['min','max','mean','sum','std'])
# 将结果转为DataFrame
df = pd.DataFrame.from_dict(stats)
# 导出CSV
df.to_csv(r'C:\data\shanghai_stats.csv',
header=True,
index_label='fid',
encoding='gbk')
这段代码有几个关键点需要注意：
- 路径中的反斜杠在Python中需要转义，或者使用原始字符串（r前缀）
- 指定encoding='gbk'是为了兼容中文路径和字段名
- index_label参数为输出CSV添加了行ID列
1.3 结果解读与可视化
统计结果通常如下表示例：
| fid | min | max | mean | sum | std |
|---|---|---|---|---|---|
| 0 | 12 | 56 | 34.2 | 684 | 8.7 |
| 1 | 8 | 62 | 29.8 | 596 | 9.3 |
对于这类数据，我推荐使用以下可视化方法：
- 用柱状图比较各区域均值
- 用箱线图展示各区域数据分布
- 用热力图呈现空间分布模式
需要模型API调用？ 免费领10W Token，多模型网关一键接入 Claude、DeepSeek 等主流模型。
2. 遥感影像聚类分析技术
2.1 K-means算法原理
K-means是遥感分类中最常用的无监督算法之一，其核心是通过迭代优化将像元划分到K个簇中。算法步骤包括：
- 初始化：随机选择K个初始聚类中心
- 分配：计算每个像元到各中心的距离，分配到最近的中心
- 更新：重新计算每个簇的均值作为新中心
- 迭代：重复2-3步直到收敛
在遥感中，距离通常指光谱特征空间中的欧氏距离。例如两个像元在NIR波段的值越接近，它们被分到同一类的概率就越高。
2.2 单波段聚类实现
以下代码演示如何对单波段土地类型数据进行聚类：
python复制from sklearn.cluster import KMeans
import numpy as np
from osgeo import gdal
# 读取栅格数据
ds = gdal.Open(r'C:\data\sh_land.tif')
arr = ds.ReadAsArray()
valid_mask = ~np.isnan(arr) # 过滤无效值
# 准备数据
X = arr[valid_mask].reshape(-1, 1)
# 执行聚类
kmeans = KMeans(n_clusters=3, random_state=42)
labels = kmeans.fit_predict(X)
# 重建结果栅格
result = np.full_like(arr, -9999, dtype=np.int32)
result[valid_mask] = labels
# 保存结果
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.CreateCopy(r'C:\data\sh_land_K.tif', ds)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(result)
out_band.SetNoDataValue(-9999)
out_ds = None
2.3 多波段聚类进阶
对于多光谱影像，我们需要同时考虑多个波段的光谱特征。以下是处理Landsat影像的完整示例：
python复制def multi_band_kmeans(band_paths, out_path, n_clusters=3):
# 读取多波段数据
bands = [gdal.Open(bp).ReadAsArray() for bp in band_paths]
stack = np.dstack(bands) # 将波段叠加为三维数组
# 构建特征矩阵
height, width, _ = stack.shape
features = stack.reshape(-1, len(band_paths))
valid_mask = ~np.isnan(features).any(axis=1)
# 执行聚类
kmeans = KMeans(n_clusters=n_clusters, random_state=42)
labels = kmeans.fit_predict(features[valid_mask])
# 输出结果
result = np.full(features.shape[0], -9999, dtype=np.int32)
result[valid_mask] = labels
result = result.reshape(height, width)
# 保存栅格
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(out_path, width, height, 1, gdal.GDT_Int32)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(result)
out_band.SetNoDataValue(-9999)
out_ds = None
# 使用示例
band_paths = [
r'C:\data\LE07_B3.TIF',
r'C:\data\LE07_B4.TIF',
r'C:\data\LE07_B5.TIF'
]
multi_band_kmeans(band_paths, r'C:\data\kmeans_result.tif', 5)
3. 实战经验与技巧
3.1 参数调优建议
- 聚类数量K的选择：
- 使用肘部法则（Elbow Method）确定最佳K值
- 对于土地覆盖分类，通常3-6个类别就能获得不错的效果
- 可通过 silhouette_score评估聚类质量
- 特征标准化：
- 不同波段的数值范围差异较大时，建议进行标准化
- 常用方法：MinMax缩放或Z-score标准化
- 初始中心选择：
- 设置random_state保证结果可复现
- 对于大数据集，可以使用k-means++初始化
3.2 常见问题排查
- 内存不足问题：
- 对于大影像，可以分块处理
- 使用memmap方式处理大型数组
- 聚类结果噪声多：
- 预处理时进行去噪（如中值滤波）
- 后处理时进行聚类结果平滑
- 分类边界不清晰：
- 尝试增加聚类数量
- 检查输入波段是否具有足够区分度
3.3 性能优化技巧
- 使用numba加速距离计算：
python复制from numba import jit
@jit(nopython=True)
def euclidean_distance(a, b):
return np.sqrt(np.sum((a - b)**2))
- 对于超大数据集：
- 使用MiniBatchKMeans替代标准KMeans
- 考虑降维处理（PCA）减少特征维度
- 并行处理：
- 设置n_jobs参数利用多核CPU
- 对于批量处理，可以使用multiprocessing
在实际项目中，我通常会先在小样本上测试算法效果，确认参数后再处理完整数据集。同时建议保存中间结果，避免重复计算。