ArcGIS Pro自相交Python代码

# -*- coding: utf-8 -*-"""ArcGIS自相交修复工具功能：1. 查找自相交点2. 创建0.01矩形面擦除原始数据，解决自相交问题"""import arcpyimport osimport uuidarcpy.env.overwriteOutput = Truedef find_self_intersections(input_fc):"""使用检查几何工具查找自相交问题    """scratch_gdb = arcpy.env.scratchGDBcheck_table = os.path.join(scratch_gdb, "check_geom_" + uuid.uuid4().hex[:8])# 检查几何arcpy.management.CheckGeometry(input_fc,check_table,"SELF_INTERSECTION")# 统计自相交数量count = 0self_intersect_oids = []with arcpy.da.SearchCursor(check_table, ["FEATURE_OID", "PROBLEM"]) as cursor:for oid, problem in cursor:if problem == "Self intersection":count += 1self_intersect_oids.append(oid)# 清理临时表arcpy.management.Delete(check_table)return count, self_intersect_oidsdef create_small_rectangle(center_x, center_y, size, sr):"""在指定位置创建一个小的矩形面    size: 矩形的边长（单位与坐标系一致）    """half_size = size / 2# 创建矩形四个角点points = [        arcpy.Point(center_x - half_size, center_y - half_size),        arcpy.Point(center_x + half_size, center_y - half_size),        arcpy.Point(center_x + half_size, center_y + half_size),        arcpy.Point(center_x - half_size, center_y + half_size),        arcpy.Point(center_x - half_size, center_y - half_size)  # 闭合]array = arcpy.Array(points)polygon = arcpy.Polygon(array, sr)return polygondef fix_self_intersection_with_erase(input_fc, output_fc, rect_size=0.01):"""通过创建小矩形面擦除来修复自相交    参数:        input_fc: 输入面要素类        output_fc: 输出要素类        rect_size: 矩形边长（默认0.01）    """desc = arcpy.Describe(input_fc)sr = desc.spatialReferencescratch_gdb = arcpy.env.scratchGDB# 第一步：检查自相交print("正在检查自相交问题...")# 使用更直接的方法：检查每个要素的自相交# 对于面要素，自相交通常发生在环相交的地方# 创建临时要素类用于存储小矩形rect_fc = os.path.join(scratch_gdb, "temp_rect_" + uuid.uuid4().hex[:8])    arcpy.management.CreateFeatureclass(scratch_gdb,        os.path.basename(rect_fc),"POLYGON",spatial_reference=sr)# 查找所有自相交点self_intersect_points = []with arcpy.da.SearchCursor(input_fc, ["OID@", "SHAPE@"]) as cursor:for oid, shape in cursor:if shape is None:continue# 获取面的所有环for part in shape:points = [p for p in part if p is not None]# 检查线段相交for i in range(len(points) - 1):for j in range(i + 2, len(points) - 1):# 避免相邻线段（它们共享一个端点）if j == i + 1:continuep1, p2 = points[i], points[i + 1]p3, p4 = points[j], points[j + 1]# 计算两条线段的交点intersect_point = line_intersection(p1, p2, p3, p4)if intersect_point:self_intersect_points.append((intersect_point.X, intersect_point.Y, oid))    print(f"发现 {len(self_intersect_points)} 个自相交点")# 创建小矩形面用于擦除with arcpy.da.InsertCursor(rect_fc, ["SHAPE@"]) as icur:for x, y, oid in self_intersect_points:rect = create_small_rectangle(x, y, rect_size, sr)icur.insertRow([rect])    print(f"创建了 {len(self_intersect_points)} 个擦除矩形")# 使用擦除工具if len(self_intersect_points) > 0:        print("正在执行擦除操作...")        arcpy.analysis.Erase(input_fc,rect_fc,output_fc)else:# 没有自相交，直接复制print("没有发现自相交问题，直接复制...")        arcpy.management.CopyFeatures(input_fc, output_fc)# 清理临时数据try:        arcpy.management.Delete(rect_fc)except:passprint(f"\n处理完成！")    print(f"原始数据: {input_fc}")    print(f"输出数据: {output_fc}")    print(f"修复了 {len(self_intersect_points)} 个自相交点")def line_intersection(p1, p2, p3, p4):"""计算两条线段的交点    返回交点坐标，如果不相交则返回None    """x1, y1 = p1.X, p1.Yx2, y2 = p2.X, p2.Yx3, y3 = p3.X, p3.Yx4, y4 = p4.X, p4.Y# 计算方向向量dx1 = x2 - x1    dy1 = y2 - y1    dx2 = x4 - x3    dy2 = y4 - y3# 计算叉积cross = dx1 * dy2 - dy1 * dx2if abs(cross) < 1e-10:# 平行或共线return None# 计算参数t = ((x3 - x1) * dy2 - (y3 - y1) * dx2) / cross    u = ((x3 - x1) * dy1 - (y3 - y1) * dx1) / cross# 检查交点是否在两条线段内if 0 <= t <= 1 and 0 <= u <= 1:# 计算交点坐标intersect_x = x1 + t * dx1        intersect_y = y1 + t * dy1return arcpy.Point(intersect_x, intersect_y)return Nonedef main():input_fc = r"C:\Users\yl\Documents\ArcGIS\Projects\MyProject17\MyProject17.gdb\m"output_fc = r"C:\Users\yl\Documents\ArcGIS\Projects\MyProject17\MyProject17.gdb\密码"rect_size = 0.01fix_self_intersection_with_erase(input_fc, output_fc, rect_size)if __name__ == "__main__":    main()