import arcpy import numpy as np from collections import defaultdict import traceback # 环境设置 arcpy.env.workspace rD:\data\heatlandtree\contourtree.gdb arcpy.env.overwriteOutput True try: arcpy.CheckOutExtension(Spatial) print(✓ Spatial Analyst 扩展已启用) except: raise RuntimeError(无法启用 Spatial Analyst 扩展请检查许可) def create_lst_contour_tree_safe(input_raster, output_fc, contour_interval2, min_area_km21.0): print(*60) print(LST等值线树生成修复版) print(*60) # 前置检查 print(\n[检查] 输入数据验证...) if not arcpy.Exists(input_raster): raise FileNotFoundError(f找不到输入栅格: {input_raster}) print(f✓ 输入栅格存在: {input_raster}) desc arcpy.Describe(input_raster) print(f✓ 数据类型: {desc.dataType}) sr desc.spatialReference if sr is None or sr.name Unknown: raise ValueError(栅格没有定义坐标系) print(f✓ 坐标系: {sr.name}) print(f✓ 坐标系类型: {sr.type}) print(f✓ 像元大小: {desc.meanCellWidth:.2f} x {desc.meanCellHeight:.2f}) # 【关键修复】投影转换 working_raster input_raster if sr.type Geographic: print(f\n⚠⚠⚠ 严重警告使用地理坐标系必须转换) print(f 当前单位度面积计算会错误) # 计算UTM带号 center_lon (desc.extent.XMin desc.extent.XMax) / 2 utm_zone int((center_lon 180) / 6) 1 print(f\n[自动转换] UTM Zone {utm_zone}N...) projected lst_projected # 使用EPSG代码 sr_code 32600 utm_zone # WGS84 UTM North arcpy.ProjectRaster_management(input_raster, projected, arcpy.SpatialReference(sr_code), BILINEAR, 30) working_raster projected # 验证转换 desc2 arcpy.Describe(working_raster) print(f 新像元大小: {desc2.meanCellWidth:.2f} x {desc2.meanCellHeight:.2f} 米) sr desc2.spatialReference # 【关键修复】重采样如果像元100m if desc.meanCellWidth 100: print(f\n[重采样] 像元过大重采样到100m...) resampled lst_resampled arcpy.Resample_management(working_raster, resampled, 100, BILINEAR) working_raster resampled # 步骤1: 平滑 print(\n[1] 平滑处理...) neighborhood arcpy.sa.NbrCircle(2, CELL) smoothed arcpy.sa.FocalStatistics(working_raster, neighborhood, MEAN) smoothed.save(smoothed_lst) # 步骤2: 获取温度范围 print(\n[2] 获取温度范围...) min_temp float(arcpy.GetRasterProperties_management(smoothed_lst, MINIMUM).getOutput(0)) max_temp float(arcpy.GetRasterProperties_management(smoothed_lst, MAXIMUM).getOutput(0)) print(f 温度范围: {min_temp:.2f}℃ ~ {max_temp:.2f}℃) # 步骤3: 生成等温线 print(f\n[3] 生成等温线间隔{contour_interval}℃...) arcpy.sa.Contour(smoothed_lst, temp_contours, contour_interval, min_temp) count int(arcpy.GetCount_management(temp_contours).getOutput(0)) print(f 生成: {count}条) if count 0: raise RuntimeError(未生成任何等温线) # 【关键修复】面积计算与自适应阈值 print(f\n[4] 面积计算与过滤...) arcpy.AddField_management(temp_contours, Area_KM2, DOUBLE) arcpy.CalculateField_management(temp_contours, Area_KM2, !shape.areaSQUAREKILOMETERS!, PYTHON3) # 获取面积分布 areas [row[0] for row in arcpy.da.SearchCursor(temp_contours, [Area_KM2])] areas.sort(reverseTrue) print(f 面积统计:) print(f 最大: {max(areas):.4f} km²) print(f 中位数: {areas[len(areas)//2]:.4f} km²) print(f 最小: {min(areas):.6f} km²) # 自适应阈值 actual_threshold min_area_km2 if max(areas) min_area_km2: # 最大面积都小于阈值使用75分位数 actual_threshold np.percentile(areas, 75) print(f\n ⚠ 最大面积({max(areas):.4f})小于阈值({min_area_km2})) print(f 自动调整为75分位数: {actual_threshold:.4f} km²) # 筛选 where fArea_KM2 {actual_threshold} arcpy.Select_analysis(temp_contours, valid_contours, where) valid_count int(arcpy.GetCount_management(valid_contours).getOutput(0)) print(f 有效等温线: {valid_count}条) if valid_count 0: # 最终保底使用最大的10个无论多小 print(f\n ⚠ 使用面积最大的10个等温线...) oids [str(row[0]) for row in arcpy.da.SearchCursor(temp_contours, [OID], sql_clause(None, ORDER BY Area_KM2 DESC))][:10] where fOBJECTID IN ({,.join(oids)}) arcpy.Select_analysis(temp_contours, valid_contours, where) valid_count 10 # 后续步骤修复几何、构建树略 print(\n[5-8] 构建拓扑结构...) # ... 此处省略使用您的原代码 ... # 清理 for temp in [lst_projected, lst_resampled, smoothed_lst, temp_contours, valid_contours]: if arcpy.Exists(temp): arcpy.Delete_management(temp) print(f\n✓ 完成输出: {output_fc}) return output_fc, {} # 主程序 if __name__ __main__: lst_input rD:\data\urmq_LST.tif if not arcpy.Exists(lst_input): print(f\n错误: 找不到文件: {lst_input}) else: try: output, tree_data create_lst_contour_tree_safe( input_rasterlst_input, output_fcLST_ContourTree_Result, contour_interval2, min_area_km22.0 ) except Exception as e: print(f\n✗ 处理失败: {e}) traceback.print_exc()
新的等值线树及其嵌套层次结构代码
import arcpy import numpy as np from collections import defaultdict import traceback # 环境设置 arcpy.env.workspace rD:\data\heatlandtree\contourtree.gdb arcpy.env.overwriteOutput True try: arcpy.CheckOutExtension(Spatial) print(✓ Spatial Analyst 扩展已启用) except: raise RuntimeError(无法启用 Spatial Analyst 扩展请检查许可) def create_lst_contour_tree_safe(input_raster, output_fc, contour_interval2, min_area_km21.0): print(*60) print(LST等值线树生成修复版) print(*60) # 前置检查 print(\n[检查] 输入数据验证...) if not arcpy.Exists(input_raster): raise FileNotFoundError(f找不到输入栅格: {input_raster}) print(f✓ 输入栅格存在: {input_raster}) desc arcpy.Describe(input_raster) print(f✓ 数据类型: {desc.dataType}) sr desc.spatialReference if sr is None or sr.name Unknown: raise ValueError(栅格没有定义坐标系) print(f✓ 坐标系: {sr.name}) print(f✓ 坐标系类型: {sr.type}) print(f✓ 像元大小: {desc.meanCellWidth:.2f} x {desc.meanCellHeight:.2f}) # 【关键修复】投影转换 working_raster input_raster if sr.type Geographic: print(f\n⚠⚠⚠ 严重警告使用地理坐标系必须转换) print(f 当前单位度面积计算会错误) # 计算UTM带号 center_lon (desc.extent.XMin desc.extent.XMax) / 2 utm_zone int((center_lon 180) / 6) 1 print(f\n[自动转换] UTM Zone {utm_zone}N...) projected lst_projected # 使用EPSG代码 sr_code 32600 utm_zone # WGS84 UTM North arcpy.ProjectRaster_management(input_raster, projected, arcpy.SpatialReference(sr_code), BILINEAR, 30) working_raster projected # 验证转换 desc2 arcpy.Describe(working_raster) print(f 新像元大小: {desc2.meanCellWidth:.2f} x {desc2.meanCellHeight:.2f} 米) sr desc2.spatialReference # 【关键修复】重采样如果像元100m if desc.meanCellWidth 100: print(f\n[重采样] 像元过大重采样到100m...) resampled lst_resampled arcpy.Resample_management(working_raster, resampled, 100, BILINEAR) working_raster resampled # 步骤1: 平滑 print(\n[1] 平滑处理...) neighborhood arcpy.sa.NbrCircle(2, CELL) smoothed arcpy.sa.FocalStatistics(working_raster, neighborhood, MEAN) smoothed.save(smoothed_lst) # 步骤2: 获取温度范围 print(\n[2] 获取温度范围...) min_temp float(arcpy.GetRasterProperties_management(smoothed_lst, MINIMUM).getOutput(0)) max_temp float(arcpy.GetRasterProperties_management(smoothed_lst, MAXIMUM).getOutput(0)) print(f 温度范围: {min_temp:.2f}℃ ~ {max_temp:.2f}℃) # 步骤3: 生成等温线 print(f\n[3] 生成等温线间隔{contour_interval}℃...) arcpy.sa.Contour(smoothed_lst, temp_contours, contour_interval, min_temp) count int(arcpy.GetCount_management(temp_contours).getOutput(0)) print(f 生成: {count}条) if count 0: raise RuntimeError(未生成任何等温线) # 【关键修复】面积计算与自适应阈值 print(f\n[4] 面积计算与过滤...) arcpy.AddField_management(temp_contours, Area_KM2, DOUBLE) arcpy.CalculateField_management(temp_contours, Area_KM2, !shape.areaSQUAREKILOMETERS!, PYTHON3) # 获取面积分布 areas [row[0] for row in arcpy.da.SearchCursor(temp_contours, [Area_KM2])] areas.sort(reverseTrue) print(f 面积统计:) print(f 最大: {max(areas):.4f} km²) print(f 中位数: {areas[len(areas)//2]:.4f} km²) print(f 最小: {min(areas):.6f} km²) # 自适应阈值 actual_threshold min_area_km2 if max(areas) min_area_km2: # 最大面积都小于阈值使用75分位数 actual_threshold np.percentile(areas, 75) print(f\n ⚠ 最大面积({max(areas):.4f})小于阈值({min_area_km2})) print(f 自动调整为75分位数: {actual_threshold:.4f} km²) # 筛选 where fArea_KM2 {actual_threshold} arcpy.Select_analysis(temp_contours, valid_contours, where) valid_count int(arcpy.GetCount_management(valid_contours).getOutput(0)) print(f 有效等温线: {valid_count}条) if valid_count 0: # 最终保底使用最大的10个无论多小 print(f\n ⚠ 使用面积最大的10个等温线...) oids [str(row[0]) for row in arcpy.da.SearchCursor(temp_contours, [OID], sql_clause(None, ORDER BY Area_KM2 DESC))][:10] where fOBJECTID IN ({,.join(oids)}) arcpy.Select_analysis(temp_contours, valid_contours, where) valid_count 10 # 后续步骤修复几何、构建树略 print(\n[5-8] 构建拓扑结构...) # ... 此处省略使用您的原代码 ... # 清理 for temp in [lst_projected, lst_resampled, smoothed_lst, temp_contours, valid_contours]: if arcpy.Exists(temp): arcpy.Delete_management(temp) print(f\n✓ 完成输出: {output_fc}) return output_fc, {} # 主程序 if __name__ __main__: lst_input rD:\data\urmq_LST.tif if not arcpy.Exists(lst_input): print(f\n错误: 找不到文件: {lst_input}) else: try: output, tree_data create_lst_contour_tree_safe( input_rasterlst_input, output_fcLST_ContourTree_Result, contour_interval2, min_area_km22.0 ) except Exception as e: print(f\n✗ 处理失败: {e}) traceback.print_exc()
