# -*- coding: utf-8 -*- import arcpy import os import datetime def calculate_spatial_relationship_optimized(): try: # 设置工作环境 arcpy.env.overwriteOutput = True # 输入参数配置 A_table = r"D:\Python\sde\KJGH.sde\SDE.ZXCQ_XXGH" # A表SDE路径 B_table = r"D:\Python\sde\KJGH.sde\SDE.ZXCQ_ZTGH" # B表SDE路径 output_folder = r"D:\test\ghcy" # 输出文件夹路径 # 创建输出文件夹 if not os.path.exists(output_folder): os.makedirs(output_folder) # 时间戳用于生成唯一文件名 timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") # 临时工作空间 temp_workspace = "in_memory" print("开始处理SDE.ZXCQ_XXGH和SDE.ZXCQ_ZTGH的空间关系...") # 验证表是否存在 if not arcpy.Exists(A_table): print("错误: A表不存在 - {0}".format(A_table)) return if not arcpy.Exists(B_table): print("错误: B表不存在 - {0}".format(B_table)) return # 检查必需字段是否存在 required_fields = ["YDYHFLDM", "YDYHFLMC"] a_fields = [f.name for f in arcpy.ListFields(A_table)] b_fields = [f.name for f in arcpy.ListFields(B_table)] for field in required_fields: if field not in a_fields: print("错误: A表中找不到 {0} 字段".format(field)) print("A表字段列表: {0}".format(a_fields)) return if field not in b_fields: print("错误: B表中找不到 {0} 字段".format(field)) print("B表字段列表: {0}".format(b_fields)) return print("字段验证通过，开始处理...") # 步骤1: 定义查询字段（只保留需要的字段） a_cursor_fields = ["OID@", "SHAPE@", "YDYHFLDM", "YDYHFLMC"] b_cursor_fields = ["OID@", "SHAPE@", "YDYHFLDM", "YDYHFLMC"] # 用于存储需要裁剪的要素 features_to_clip = [] # 创建B表的临时图层 print("正在创建临时图层...") B_layer = "B_table_layer" arcpy.MakeFeatureLayer_management(B_table, B_layer) print("开始遍历A表要素...") with arcpy.da.SearchCursor(A_table, a_cursor_fields) as a_cursor: a_count = 0 processed_count = 0 for a_row in a_cursor: a_count += 1 a_oid = a_row[0] a_geometry = a_row[1] a_ydyhfldm = a_row[2] a_ydyhflmc = a_row[3] if a_count % 100 == 0: print("已处理 {0} 个A表要素，找到 {1} 个需要裁剪的要素".format(a_count, processed_count)) # 检查几何是否有效 if a_geometry is None or a_geometry.length == 0: continue try: # 使用SelectLayerByLocation选择与A表要素相交的B表要素 arcpy.SelectLayerByLocation_management(B_layer, "INTERSECT", a_geometry) # 获取选中的要素数量 result = arcpy.GetCount_management(B_layer) selected_count = int(result.getOutput(0)) if selected_count > 0: # 遍历选中的B表要素 with arcpy.da.SearchCursor(B_layer, b_cursor_fields) as b_cursor: for b_row in b_cursor: b_oid = b_row[0] b_geometry = b_row[1] b_ydyhfldm = b_row[2] b_ydyhflmc = b_row[3] # 检查几何是否有效 if b_geometry is None or b_geometry.length == 0: continue # 处理YDYHFLDM字段比较 a_ydyhfldm_str = str(a_ydyhfldm) if a_ydyhfldm is not None else "" b_ydyhfldm_str = str(b_ydyhfldm) if b_ydyhfldm is not None else "" # 判断是否需要执行后续操作 need_process = False if b_ydyhfldm_str.endswith('0000'): # 如果b_ydyhfldm_str以0000结尾，比较前两位 if len(a_ydyhfldm_str) >= 2 and len(b_ydyhfldm_str) >= 2: if a_ydyhfldm_str[:2] != b_ydyhfldm_str[:2]: need_process = True else: # 其他情况，直接比较整个字符串 if a_ydyhfldm_str != b_ydyhfldm_str: need_process = True if need_process: processed_count += 1 if processed_count % 50 == 0: print( " 发现第 {0} 个不一致要素: A表OID={1}, B表OID={2}".format(processed_count, a_oid, b_oid)) # 获取相交部分几何 try: intersect_geom = a_geometry.intersect(b_geometry, 4) # 4表示平面相交 # 检查相交几何是否有效 if intersect_geom is not None and intersect_geom.length > 0: # 创建要素信息 a_feature = { 'geometry': intersect_geom, 'a_ydyhfldm': a_ydyhfldm_str, 'a_ydyhflmc': str(a_ydyhflmc) if a_ydyhflmc is not None else "", 'b_ydyhfldm': b_ydyhfldm_str, 'b_ydyhflmc': str(b_ydyhflmc) if b_ydyhflmc is not None else "" } features_to_clip.append(a_feature) except Exception as e: print(" 计算相交几何时出错: {0}".format(e)) continue except Exception as e: print(" 空间查询B表时出错: {0}".format(e)) continue # 清理临时图层 try: arcpy.Delete_management(B_layer) except: pass print("处理完成！总共处理了 {0} 个A表要素，找到 {1} 个需要裁剪的要素".format(a_count, len(features_to_clip))) # 如果有需要处理的要素，进行导出 if features_to_clip: print("开始导出 {0} 个裁剪要素...".format(len(features_to_clip))) # 创建输出要素类 output_fc_name = "clipped_features_{0}".format(timestamp) output_fc_path = os.path.join(temp_workspace, output_fc_name) # 获取空间参考 spatial_ref = arcpy.Describe(A_table).spatialReference # 创建输出要素类 arcpy.CreateFeatureclass_management(temp_workspace, output_fc_name, "POLYGON", None, "DISABLED", "DISABLED", spatial_ref) # 添加必需字段 arcpy.AddField_management(output_fc_path, "A_YDYHFLDM", "TEXT", "", "", 50) arcpy.AddField_management(output_fc_path, "A_YDYHFLMC", "TEXT", "", "", 100) arcpy.AddField_management(output_fc_path, "B_YDYHFLDM", "TEXT", "", "", 50) arcpy.AddField_management(output_fc_path, "B_YDYHFLMC", "TEXT", "", "", 100) # 准备插入字段列表 output_fields = [ "SHAPE@", "A_YDYHFLDM", "A_YDYHFLMC", "B_YDYHFLDM", "B_YDYHFLMC" ] print("输出字段: {0}".format(output_fields)) # 插入数据 print("正在插入裁剪后的要素...") inserted_count = 0 with arcpy.da.InsertCursor(output_fc_path, output_fields) as insert_cursor: for feature in features_to_clip: try: # 构建插入行数据 row_data = [ feature['geometry'], feature['a_ydyhfldm'], feature['a_ydyhflmc'], feature['b_ydyhfldm'], feature['b_ydyhflmc'] ] insert_cursor.insertRow(row_data) inserted_count += 1 if inserted_count % 100 == 0: print("已插入 {0} 个要素".format(inserted_count)) except Exception as e: print("插入要素时出错: {0}".format(e)) continue # 导出为Shapefile output_shp = os.path.join(output_folder, "clipped_results_{0}.shp".format(timestamp)) arcpy.CopyFeatures_management(output_fc_path, output_shp) mjField = "TBMJ" arcpy.AddField_management(output_shp, mjField, "DOUBLE") arcpy.CalculateField_management(output_shp, mjField, "!shape.area@SQUAREMETERS!", "PYTHON_9.3") print("结果已导出到: {0}".format(output_shp)) print("总共导出了 {0} 个要素".format(inserted_count)) # 清理临时数据 try: arcpy.Delete_management(output_fc_path) except: pass else: print("没有找到需要裁剪的要素") print("处理完成！") except Exception as e: print("处理过程中发生错误: {0}".format(e)) import traceback print(traceback.format_exc()) if __name__ == "__main__": calculate_spatial_relationship_optimized()