Python scikit-learn KMeans 与 ArcPy 结合：GIS面数据4属性聚类与3步自动化制图
地理信息系统（GIS）分析中，聚类技术常被用于发现空间数据的潜在模式。传统方法多依赖商业软件内置工具，但通过Python开源生态与ArcPy的结合，我们可以构建更灵活、可复现的自动化分析流程。本文将展示如何用scikit-learn的KMeans算法对GIS面数据的四个属性进行聚类，并通过ArcPy实现结果回写与自动化制图。
1. 环境准备与数据加载
在开始前，确保已安装以下Python库：
- scikit-learn (≥1.3.0)
- arcpy (ArcGIS Pro自带)
- pandas (≥1.5.0)
- numpy (≥1.21.0)
import arcpy import pandas as pd from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler加载面状GIS数据（以Shapefile为例）：
input_features = "path/to/your/input.shp" # 替换为实际路径 fields = ["GDP_per_capita", "Pop_density", "Rainfall", "Area"] # 四个待聚类属性使用pandas读取属性表：
arr = arcpy.da.TableToNumPyArray(input_features, fields) df = pd.DataFrame(arr) print(df.head()) # 检查数据前五行2. 数据预处理与KMeans聚类
2.1 数据标准化
由于不同属性量纲差异大，需进行Z-score标准化：
scaler = StandardScaler() scaled_data = scaler.fit_transform(df)2.2 确定最佳聚类数
通过肘部法则（Elbow Method）寻找最佳K值：
inertia = [] for k in range(1, 10): kmeans = KMeans(n_clusters=k, random_state=42) kmeans.fit(scaled_data) inertia.append(kmeans.inertia_) # 可视化（需matplotlib） import matplotlib.pyplot as plt plt.plot(range(1, 10), inertia, marker='o') plt.xlabel('Number of clusters') plt.ylabel('Inertia') plt.show()2.3 执行KMeans聚类
假设通过肘部法则确定K=5：
kmeans = KMeans(n_clusters=5, random_state=42) clusters = kmeans.fit_predict(scaled_data) df['CLUSTER_ID'] = clusters # 添加聚类结果列3. 结果回写与自动化制图
3.1 将聚类结果写回原始数据
使用ArcPy的更新游标：
# 添加CLUSTER_ID字段 arcpy.AddField_management(input_features, "CLUSTER_ID", "SHORT") # 更新字段值 with arcpy.da.UpdateCursor(input_features, ["OID@"] + fields + ["CLUSTER_ID"]) as cursor: for row in cursor: oid = row[0] row[-1] = df.loc[df.index[df["OID"] == oid], "CLUSTER_ID"].values[0] cursor.updateRow(row)3.2 自动化分类渲染
创建颜色映射并应用分类渲染：
# 定义颜色方案（RGB元组） color_ramp = [ (56, 108, 176), # 蓝色 (191, 213, 151), # 浅绿 (223, 194, 125), # 黄色 (217, 95, 14), # 橙色 (136, 46, 114) # 紫色 ] # 转换为ArcGIS颜色对象 colors = [arcpy.Color(*rgb) for rgb in color_ramp] # 创建唯一值渲染 lyr = arcpy.mapping.Layer(input_features) sym = lyr.symbology sym.renderer = "UniqueValueRenderer" sym.renderer.fields = ["CLUSTER_ID"] for cluster_id, color in zip(range(5), colors): sym.renderer.addValue(cluster_id, str(cluster_id)) sym.renderer.symbols[cluster_id].color = color lyr.symbology = sym3.3 导出专题地图
设置布局元素并导出：
mxd = arcpy.mapping.MapDocument("CURRENT") df = arcpy.mapping.ListDataFrames(mxd)[0] arcpy.mapping.AddLayer(df, lyr) # 添加图例 legend = arcpy.mapping.ListLayoutElements(mxd, "LEGEND_ELEMENT")[0] legend.autoAdd = True # 导出地图 output_map = "path/to/output_map.pdf" arcpy.mapping.ExportToPDF(mxd, output_map)4. 完整流程封装与优化
将上述步骤封装为可复用的函数：
def cluster_gis_features(input_fc, fields, output_fc, n_clusters=5): """全流程GIS属性聚类函数""" # 数据读取与预处理 arr = arcpy.da.TableToNumPyArray(input_fc, fields) df = pd.DataFrame(arr) scaler = StandardScaler() scaled_data = scaler.fit_transform(df) # 聚类分析 kmeans = KMeans(n_clusters=n_clusters, random_state=42) df['CLUSTER_ID'] = kmeans.fit_predict(scaled_data) # 结果回写 if arcpy.Exists(output_fc): arcpy.Delete_management(output_fc) arcpy.CopyFeatures_management(input_fc, output_fc) arcpy.AddField_management(output_fc, "CLUSTER_ID", "SHORT") with arcpy.da.UpdateCursor(output_fc, ["OID@", "CLUSTER_ID"]) as cursor: for row in cursor: oid = row[0] row[1] = df.loc[df["OID"] == oid, "CLUSTER_ID"].values[0] cursor.updateRow(row) return output_fc实际调用示例：
output_features = "path/to/output_clustered.shp" cluster_gis_features(input_features, fields, output_features, n_clusters=5)5. 高级技巧与注意事项
5.1 处理大型数据集
当数据量较大时（>10万要素），可采用以下优化：
# 使用MiniBatchKMeans加速 from sklearn.cluster import MiniBatchKMeans kmeans = MiniBatchKMeans(n_clusters=5, batch_size=1000) # 分块处理 chunk_size = 5000 for i in range(0, len(df), chunk_size): chunk = scaled_data[i:i + chunk_size] kmeans.partial_fit(chunk)5.2 属性权重调整
若某些属性更重要，可在标准化前设置权重：
weights = [0.4, 0.3, 0.2, 0.1] # 对应四个属性的权重 weighted_data = scaled_data * weights5.3 结果验证
计算轮廓系数评估聚类质量：
from sklearn.metrics import silhouette_score score = silhouette_score(scaled_data, clusters) print(f"轮廓系数: {score:.3f}")提示：轮廓系数范围在[-1,1]，越接近1表示聚类效果越好
6. 实际应用案例
假设我们有一个城市区块数据集，包含以下属性：
- 人口密度（人/平方公里）
- 平均房价（万元/平米）
- 绿化率（%）
- 交通站点数量
通过本方法可识别出城市发展类型：
| 聚类ID | 典型特征 | 城市区域类型 |
|---|---|---|
| 0 | 高密度、高房价、高绿化 | 高端住宅区 |
| 1 | 中密度、中房价、中等交通 | 成熟社区 |
| 2 | 低密度、低房价、低交通 | 待开发区域 |
| 3 | 高密度、低房价、低绿化 | 老旧城区 |
| 4 | 中密度、高房价、高交通 | 商业中心区 |
这种分析可辅助城市规划决策，如：
- 识别需要基础设施升级的区域
- 发现潜在的土地利用优化空间
- 指导公共服务资源分配
7. 扩展应用方向
本方法可进一步扩展至：
- 时序聚类：分析多个时间点的属性变化模式
- 空间约束聚类：结合地理邻接关系改进结果
- 多维可视化：用PCA降维后绘制聚类分布
# 示例：PCA降维可视化 from sklearn.decomposition import PCA pca = PCA(n_components=2) coords = pca.fit_transform(scaled_data) plt.scatter(coords[:,0], coords[:,1], c=clusters, cmap='viridis') plt.colorbar(label='Cluster ID') plt.title('2D PCA Projection of Clusters') plt.show()通过Python与ArcPy的结合，我们实现了从数据预处理到最终地图输出的全流程自动化。这种方法不仅提高了分析效率，还保证了结果的可重复性，特别适合需要定期更新的分析任务。
