python 进行空间自相关-聚类-异常值分析

import arcpy
import os
import traceback

# 设置工作环境
arcpy.env.overwriteOutput = True
input_gdb = r”E:\zmx\2025\mxphdproject\tem_data\SpeciesAnalysis.gdb”
output_gdb = r”E:\zmx\2025\mxphdproject\tem_data\SpeciesAnalysis.gdb”

def perform_spatial_analysis(input_feature, count_field, grid_size):
“””执行空间分析三部曲：Moran’s I、热点分析、聚类和异常值分析（保留原始COType字段名）”””
try:
# 1. 空间自相关分析 (Global Moran’s I)
print(f”\n正在对 {grid_size} 网格进行空间自相关分析…”)

moran_result = arcpy.stats.SpatialAutocorrelation(
Input_Feature_Class=input_feature,
Input_Field=count_field,
Generate_Report=”NO_REPORT”,
Conceptualization_of_Spatial_Relationships=”CONTIGUITY_EDGES_CORNERS”,
Distance_Method=”EUCLIDEAN_DISTANCE”,
Standardization=”ROW”
)

print(f” 空间自相关分析完成”)
print(f” Moran’s I 指数: {moran_result[0]}”)
print(f” Z得分: {moran_result[1]}”)
print(f” P值: {moran_result[2]}”)

# 2. 热点分析 (Getis-Ord Gi*)
print(f”\n正在对 {grid_size} 网格进行热点分析…”)
hotspot_result = os.path.join(output_gdb, f”hotspot_{grid_size}”)
arcpy.stats.HotSpots(
Input_Feature_Class=input_feature,
Input_Field=count_field,
Output_Feature_Class=hotspot_result,
Conceptualization_of_Spatial_Relationships=”CONTIGUITY_EDGES_CORNERS”,
Distance_Method=”EUCLIDEAN_DISTANCE”,
Standardization=”NONE”
)
print(f” 热点分析完成，结果保存在: {hotspot_result}”)

# 3. 聚类和异常值分析 (Anselin Local Moran’s I)
print(f”\n正在对 {grid_size} 网格进行聚类和异常值分析…”)
cluster_result = os.path.join(output_gdb, f”cluster_{grid_size}”)

arcpy.stats.ClustersOutliers(
Input_Feature_Class=input_feature,
Input_Field=count_field,
Output_Feature_Class=cluster_result,
Conceptualization_of_Spatial_Relationships=”CONTIGUITY_EDGES_CORNERS”,
Distance_Method=”EUCLIDEAN_DISTANCE”,
Standardization=”ROW”,
Distance_Band_or_Threshold_Distance=None,
Weights_Matrix_File=None,
Apply_False_Discovery_Rate__FDR__Correction=”NO_FDR”,
Number_of_Permutations=99,
number_of_neighbors=8
)

print(f” 聚类和异常值分析完成，结果保存在: {cluster_result}”)

# 字段重命名策略（保留COType原始名称）
for result in [hotspot_result, cluster_result]:
if arcpy.Exists(result):
# 仅修改可能冲突的字段
try:
arcpy.management.AlterField(
in_table=result,
field=”GI_STAR”,
new_field_name=f”GI_STAR_{grid_size}”,
new_field_alias=f”Gi*统计量_{grid_size}”
)
except:
print(f” {result} 中不存在GI_STAR字段，跳过重命名”)

try:
arcpy.management.AlterField(
in_table=result,
field=”LOCALMI”,
new_field_name=f”LOCALMI_{grid_size}”,
new_field_alias=f”局部Moran’I_{grid_size}”
)
except:
print(f” {result} 中不存在LOCALMI字段，跳过重命名”)

# 明确保留COType字段不修改
print(f” 已保留 {result} 中的COType字段原始名称”)

# 自动化符号系统（使用原始COType字段）
apply_symbology(cluster_result, grid_size)

return True

except arcpy.ExecuteError as e:
print(f”执行 {grid_size} 网格空间分析时出错: {str(e)}”)
print(arcpy.GetMessages(2))
return False
except Exception as e:
print(f”执行 {grid_size} 网格空间分析时出现意外错误: {str(e)}”)
print(traceback.format_exc())
return False

def apply_symbology(cluster_layer, grid_size):
“””为聚类结果应用符号系统（使用原始COType字段）”””
try:
aprx = arcpy.mp.ArcGISProject(“CURRENT”)
map = aprx.listMaps()[0]

# 检查图层是否已存在
existing_layers = [lyr.name for lyr in map.listLayers()]
layer_name = os.path.basename(cluster_layer)

if layer_name not in existing_layers:
map.addDataFromPath(cluster_layer)

# 获取图层对象
cluster_lyr = map.listLayers(layer_name)[0]

# 设置符号系统
sym = cluster_lyr.symbology
if hasattr(sym, ‘renderer’):
sym.updateRenderer(‘UniqueValueRenderer’)
# 关键点：使用原始COType字段
sym.renderer.fields = [“COType”]

# 定义分类颜色（与COType原始值匹配）
color_mapping = {
“HH”: {“color”: [255, 0, 0], “label”: “高-高聚类”},
“LL”: {“color”: [0, 0, 255], “label”: “低-低聚类”},
“HL”: {“color”: [255, 170, 255], “label”: “高-低异常”},
“LH”: {“color”: [115, 178, 255], “label”: “低-高异常”}
}

# 清除默认值并添加自定义分类
sym.renderer.removeAllValues()
for val, props in color_mapping.items():
sym.renderer.addValue({
“value”: val,
“label”: props[“label”],
“symbol”: {
“type”: “CIMSymbolReference”,
“symbol”: {
“type”: “CIMFillSymbol”,
“symbolLayers”: [{
“type”: “CIMSolidFill”,
“color”: {
“type”: “CIMRGBColor”,
“values”: props[“color”] + [100] # RGB + 透明度
}
}]
}
}
})

cluster_lyr.symbology = sym

aprx.save()
print(f” 已基于COType字段应用符号系统”)
except Exception as e:
print(f” 符号系统自动化失败: {str(e)}”)
print(traceback.format_exc())

# 定义要分析的网格图层和对应字段
grid_layers = {
“50km”: “grid_50km_Project_SpeciesCount”,
“100km”: “grid_100km_Project_SpeciesCount”,
“150km”: “grid_150km_Project_SpeciesCount”,
“200km”: “grid_200km_Project_SpeciesCount”
}

# 执行分析
for grid_size, layer_name in grid_layers.items():
input_layer = os.path.join(input_gdb, layer_name)

# 检查图层是否存在
if not arcpy.Exists(input_layer):
print(f”警告: 图层 {layer_name} 不存在，跳过分析”)
continue

# 检查字段是否存在
count_field = f”SP_COUNT_{grid_size}”
fields = [f.name for f in arcpy.ListFields(input_layer)]
if count_field not in fields:
print(f”警告: 字段 {count_field} 在图层 {layer_name} 中不存在，跳过分析”)
continue

print(f”\n{‘=’*50}”)
print(f”开始处理 {grid_size} 网格: {layer_name}”)
print(f”使用字段: {count_field}”)

# 执行空间分析
success = perform_spatial_analysis(input_layer, count_field, grid_size)

if success:
print(f”\n{grid_size} 网格空间分析完成!”)
else:
print(f”\n{grid_size} 网格空间分析失败!”)

print(“\n所有网格空间分析完成!”)