onemap-implementationSupervision-system-python/ghcy/ghcy.py

# -*- 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()