|
@@ -0,0 +1,206 @@
|
|
|
|
|
+# -*- 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()
|