Py学习  »  Python

ArcGIS Pro自相交Python代码

GISAI • 1 周前 • 62 次点击  
# -*- coding: utf-8 -*-
"""
ArcGIS自相交修复工具
功能:
1. 查找自相交点
2. 创建0.01矩形面擦除原始数据,解决自相交问题
"""
import arcpy
import os
import uuid

arcpy.env.overwriteOutput = True


def find_self_intersections(input_fc):
"""
使用检查几何工具查找自相交问题
    """
scratch_gdb = arcpy.env.scratchGDB
check_table = os.path.join(scratch_gdb"check_geom_" + uuid.uuid4().hex[:8])

检查几何
arcpy.management.CheckGeometry(
input_fc,
check_table,
"SELF_INTERSECTION"
)

统计自相交数量
count 0
self_intersect_oids = []

with arcpy.da.SearchCursor(check_table, ["FEATURE_OID""PROBLEM"]) as cursor:
for oidproblem in cursor:
if problem == "Self intersection":
count += 1
self_intersect_oids.append(oid)

清理临时表
arcpy.management.Delete(check_table)

return countself_intersect_oids


def create_small_rectangle(center_xcenter_ysizesr):
"""
在指定位置创建一个小的矩形面
    size: 矩形的边长(单位与坐标系一致)
    """
half_size size 2

创建矩形四个角点
points = [
        arcpy.Point(center_x half_sizecenter_y half_size),
        arcpy.Point(center_x half_sizecenter_y half_size),
        arcpy.Point(center_x  half_sizecenter_y half_size),
        arcpy.Point(center_x half_sizecenter_y half_size),
        arcpy.Point(center_x half_sizecenter_y half_size)  闭合
]

array = arcpy.Array(points)
polygon = arcpy.Polygon(arraysr)

return polygon


def fix_self_intersection_with_erase(input_fcoutput_fcrect_size=0.01):
"""
通过创建小矩形面擦除来修复自相交

    参数:
        input_fc: 输入面要素类
        output_fc: 输出要素类
        rect_size: 矩形边长(默认0.01
    """

desc = arcpy.Describe(input_fc)
sr desc.spatialReference

scratch_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 oidshape in cursor:
if shape is None:
continue

获取面的所有环
for part in shape:
points = [ for in part if is not None]

检查线段相交
for in range(len(points) - 1):
for in range(2, len(points) - 1):
避免相邻线段(它们共享一个端点)
if == 1:
continue

p1p2 points[i], points[1]
p3p4 points[j], points[1]

计算两条线段的交点
intersect_point = line_intersection(p1p2p3p4)

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 xyoid in self_intersect_points:
rect = create_small_rectangle(xyrect_sizesr)
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_fcoutput_fc)

清理临时数据
try:
        arcpy.management.Delete(rect_fc)
except:
pass

print(f"\n处理完成!")
    print(f"原始数据{input_fc}")
    print(f"输出数据{output_fc}")
    print(f"修复了 {len(self_intersect_points)个自相交点")


def line_intersection(p1p2p3p4):
"""
计算两条线段的交点
    返回交点坐标,如果不相交则返回None
    """
x1y1 p1.X, p1.Y
x2y2 p2.X, p2.Y
x3y3 p3.X, p3.Y
x4y4 p4.X, p4.Y

计算方向向量
dx1 x2 x1
    dy1 y2 y1
    dx2 x4 x3
    dy2 y4 y3

计算叉积
cross dx1 dy2 dy1 dx2

if abs(cross) 1e-10:
平行或共线
return None

计算参数
= ((x3 x1) * dy2 - (y3 y1) * dx2) / cross
    u = ((x3 x1) * dy1 - (y3 y1) * dx1) / cross

检查交点是否在两条线段内
if <= <= and <= <= 1:
计算交点坐标
intersect_x x1 dx1
        intersect_y y1 dy1
return arcpy.Point(intersect_xintersect_y)

return None


def 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.01

fix_self_intersection_with_erase(input_fcoutput_fcrect_size)


if __name__ == "__main__":
    main()

Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/196188