python在Arcgis Pro中 多边形锐角识别与切割脚本笔记

   

一、脚本定位与应用场景

该脚本基于 ArcGIS Pro Python 3 环境开发，是一款自动化矢量数据处理工具，核心功能是识别多边形要素中的锐角顶点并按规则切割，解决因锐角导致的 GIS 数据拓扑错误、地图可视化失真、空间分析偏差等问题，适用于建筑规划数据标准化、测绘矢量图预处理、土地利用数据优化等场景。

二、脚本结构拆解与代码解读

（一）环境配置与参数定义

1. 基础环境初始化

# -*- coding: utf-8 -*-
import arcpy  # 导入ArcGIS地理数据处理库
import math   # 导入数学计算库（用于角度、距离计算）
import os     # 导入文件路径处理库
import sys    # 导入系统控制库（用于递归限制、程序退出）
from arcpy import env  # 导入ArcGIS环境配置模块# 环境核心设置
env.overwriteOutput = True  # 允许覆盖已存在的输出文件，避免重复运行时文件冲突
sys.setrecursionlimit(1000000)  # 提升递归限制，处理顶点数量多的复杂多边形（如复杂建筑轮廓）

解读：这部分是脚本运行的基础保障，env.overwriteOutput避免手动删除旧文件的麻烦，sys.setrecursionlimit则针对复杂图形的顶点遍历，防止因递归深度不足报错。

2. 关键参数配置

# 输入输出路径（需替换为实际路径，r前缀避免转义字符问题）
input_shp = r"C:\Users\Administrator\Desktop\test\01测试数据\测试数据1\锐角测试.shp"
output_mark = r"C:\Users\Administrator\Desktop\test\内角标示.shp"  # 锐角顶点标记输出
output_cut = r"C:\Users\Administrator\Desktop\test\锐角切割.shp"    # 切割结果输出
output_mid = r"C:\Users\Administrator\Desktop\test\中间点.shp"      # 调试用中间点输出# 核心阈值参数（可根据业务需求调整）
angle_threshold = 30.0  # 锐角判断阈值（单位：度），小于该值的角判定为锐角
max_area = 1.0          # 切割后三角形最大允许面积（单位：平方米），避免切割块过大

解读：路径参数需使用绝对路径，angle_threshold和max_area是脚本的核心控制参数 —— 前者决定 “什么是锐角”，后者决定 “切割到多大合适”，需根据数据精度要求调整（如高精度测绘数据可将max_area设为 0.5）。

（二）数据校验：避免后续计算出错

# 校验1：输入文件是否存在
if not arcpy.Exists(input_shp):arcpy.AddError(f"输入文件不存在：{input_shp}")  # 向ArcGIS输出错误信息sys.exit(1)  # 终止程序# 校验2：坐标系是否为投影坐标系（确保距离、面积计算单位为米）
desc = arcpy.Describe(input_shp)  # 获取输入文件的描述信息
sr = desc.spatialReference        # 提取坐标系信息
if sr.type != "Projected":arcpy.AddError(f"请使用投影坐标系（如UTM）！当前坐标系：{sr.name}")sys.exit(1)

解读：这是脚本的 “安全检查”—— 若输入文件路径错误，或用经纬度（地理坐标系）计算面积（会出现 “度 ²” 这种无意义单位），脚本会直接报错并退出，避免生成无效结果。

（三）输出图层创建：提前搭建结果存储结构

# 清理旧数据（若之前运行过脚本，删除旧输出文件）
for shp in [output_mark, output_cut, output_mid]:if arcpy.Exists(shp):arcpy.Delete_management(shp)# 1. 创建锐角顶点标记图层（点要素）
arcpy.CreateFeatureclass_management(os.path.dirname(output_mark),  # 输出文件夹路径os.path.basename(output_mark), # 输出文件名"POINT",                       # 要素类型（点）spatial_reference=sr           # 与输入文件保持一致的坐标系
)
# 为顶点标记图层添加字段（存储原要素ID和内角值）
arcpy.AddField_management(output_mark, "ORIG_FID", "LONG")    # 关联原多边形的ID
arcpy.AddField_management(output_mark, "INT_ANGLE", "DOUBLE") # 存储锐角的具体角度# 2. 创建切割结果图层（多边形要素）
arcpy.CreateFeatureclass_management(os.path.dirname(output_cut), os.path.basename(output_cut), "POLYGON",  # 要素类型（多边形）spatial_reference=sr
)
# 为切割结果图层添加字段（记录关键信息）
arcpy.AddField_management(output_cut, "ORIG_FID", "LONG")    # 关联原多边形ID
arcpy.AddField_management(output_cut, "VERTEX_ID", "LONG")   # 关联原顶点ID
arcpy.AddField_management(output_cut, "AREA", "DOUBLE")      # 切割后三角形的面积
arcpy.AddField_management(output_cut, "MIN_EDGE", "DOUBLE")  # 切割时用的最短边长度# 3. 创建调试用中间点图层（点要素，用于查看分割点位置）
arcpy.CreateFeatureclass_management(os.path.dirname(output_mid), os.path.basename(output_mid), "POINT", spatial_reference=sr
)
arcpy.AddField_management(output_mid, "TYPE", "TEXT", 10)  # 标记点类型（顶点/分割点）

解读：脚本提前创建 3 个输出图层，并定义专属字段，确保结果能 “带属性存储”—— 比如切割后的三角形，能通过ORIG_FID找到它来自哪个原多边形，通过AREA验证是否符合面积要求，方便后续分析。

（四）核心函数：实现锐角识别与切割的 “大脑”

1. 内角计算函数（calc_angle）：判断是否为锐角

def calc_angle(p_prev, p_curr, p_next):"""计算多边形顶点p_curr的内角（核心逻辑）"""try:# 步骤1：计算当前顶点与相邻顶点的向量（v1：前一顶点→当前顶点；v2：后一顶点→当前顶点）v1 = (p_prev.X - p_curr.X, p_prev.Y - p_curr.Y)v2 = (p_next.X - p_curr.X, p_next.Y - p_curr.Y)# 步骤2：计算向量的点积（求夹角余弦值）和叉积（判断凹凸性）dot = v1[0] * v2[0] + v1[1] * v2[1]  # 点积：v1·v2 = |v1||v2|cosθcross = v1[0] * v2[1] - v1[1] * v2[0]# 叉积：v1×v2 = |v1||v2|sinθ（符号判断方向）# 步骤3：计算向量长度（避免除以0）len1, len2 = math.hypot(*v1), math.hypot(*v2)  # hypot计算直角三角形斜边（即向量长度）if len1 < 1e-6 or len2 < 1e-6:  # 若某条边长度接近0（重复顶点），返回0度return 0.0# 步骤4：计算角度（cosθ=点积/(|v1||v2|)，再转角度）cos_theta = max(-1.0, min(1.0, dot / (len1 * len2)))  # 限制cosθ在[-1,1]（避免计算误差）angle = math.degrees(math.acos(cos_theta))  # 弧度转角度# 步骤5：判断内角方向（凸角/凹角）：叉积≥0为凸角（内角=计算值），否则为凹角（内角=360-计算值）return angle if cross >= 0 else 360 - angleexcept Exception as e:arcpy.AddWarning(f"角度计算错误: {str(e)}")  # 捕获异常，输出警告（不终止程序）return 0.0

解读：这是脚本的 “核心算法”—— 通过向量运算精准计算多边形内角：

• 点积决定 “夹角大小”，叉积决定 “角的方向”（凸多边形的内角都是凸角，凹多边形会有大于 180° 的凹角）；
• 加入1e-6误差控制和异常捕获，避免因数据瑕疵（如重复顶点）导致程序崩溃。

2. 距离计算函数（get_distance）：求两点间直线距离

def get_distance(p1, p2):"""计算两个点（p1、p2）之间的直线距离"""try:return math.hypot(p1.X - p2.X, p1.Y - p2.Y)  # 本质是勾股定理：距离=√[(x2-x1)²+(y2-y1)²]except Exception as e:arcpy.AddWarning(f"距离计算错误: {str(e)}")return 0.0

解读：用于计算锐角相邻两边的长度，后续会取 “最短边长度” 作为切割时的 “基准长度”，确保切割出的两条分割边等长。

3. 分割点计算函数（get_split_point）：确定切割位置

def get_split_point(p_start, p_end, length):"""在p_start到p_end的边上，按指定长度length计算分割点（从p_start向p_end截取）"""try:dx = p_end.X - p_start.X  # 边在X轴方向的增量dy = p_end.Y - p_start.Y  # 边在Y轴方向的增量total_len = math.hypot(dx, dy)  # 边的总长度if total_len < 1e-6:  # 若边长度接近0，直接返回终点（避免错误）return p_end# 计算截取比例（若length超过总长度，取1.0，避免分割点超出边范围）ratio = min(length / total_len, 1.0)# 计算分割点坐标：p_start + 比例*增量return arcpy.Point(p_start.X + dx * ratio,p_start.Y + dy * ratio)except Exception as e:arcpy.AddWarning(f"分割点计算错误: {str(e)}")return p_end

解读：比如锐角相邻两边中，最短边长度为 2 米，就会在两条边上各从顶点截取 2 米，得到两个分割点 —— 这两个点和原顶点会构成一个等腰三角形，确保切割形状规则。

（五）主处理逻辑：串联所有功能，实现自动化处理

try:# 打开3个输出图层的插入游标（用于写入结果，Python 3上下文管理器自动释放资源）with arcpy.da.InsertCursor(output_mark, ["SHAPE@", "ORIG_FID", "INT_ANGLE"]) as mark_cursor, \arcpy.da.InsertCursor(output_cut, ["SHAPE@", "ORIG_FID", "VERTEX_ID", "AREA", "MIN_EDGE"]) as cut_cursor, \arcpy.da.InsertCursor(output_mid, ["SHAPE@", "TYPE"]) as mid_cursor:# 遍历输入多边形要素（读取每个多边形的ID和形状）with arcpy.da.SearchCursor(input_shp, ["OID@", "SHAPE@"]) as search_cursor:for oid, shape in search_cursor:arcpy.AddMessage(f"\n处理要素 ID: {oid}")  # 向ArcGIS输出处理进度# 处理多边形的每个“环”（复杂多边形可能有外环、内环，如带孔洞的地块）for part in shape:points = list(part)  # 提取环的所有顶点num_points = len(points)if num_points < 3:  # 少于3个顶点无法构成多边形，跳过arcpy.AddWarning(f"要素 {oid} 的部分顶点数不足3个，跳过")continue# 处理闭合点（ArcGIS中多边形会自动在末尾添加与第一个点相同的闭合点，需删除）if num_points > 1 and points[0].equals(points[-1]):valid_points = points[:-1]  # 去除重复的最后一个点else:valid_points = pointsnum_valid = len(valid_points)if num_valid < 3:arcpy.AddWarning(f"要素 {oid} 的有效顶点数不足3个，跳过")continue# 遍历每个有效顶点，判断是否为锐角并切割for i in range(num_valid):# 获取当前顶点的前一个、当前、后一个顶点（首尾相连，如第0个顶点的前一个是最后一个）p_prev = valid_points[i-1] if i > 0 else valid_points[-1]p_curr = valid_points[i]p_next = valid_points[(i+1) % num_valid]# 步骤1：计算当前顶点的内角angle = calc_angle(p_prev, p_curr, p_next)arcpy.AddMessage(f"顶点 {i} 内角：{angle:.2f}°")# 步骤2：若为锐角（0<角度<阈值），执行标记和切割if 0 < angle < angle_threshold:# 标记锐角顶点（写入顶点标记图层）mark_cursor.insertRow([p_curr, oid, round(angle, 2)])# 标记顶点到调试图层mid_cursor.insertRow([p_curr, "顶点"])# 步骤3：计算相邻两边的长度，取最短边作为切割基准len_prev = get_distance(p_curr, p_prev)  # 当前顶点到前一顶点的长度len_next = get_distance(p_curr, p_next)  # 当前顶点到后一顶点的长度min_len = min(len_prev, len_next)        # 最短边长度if min_len < 1e-6:  # 若最短边接近0（重复顶点），跳过切割arcpy.AddWarning(f"顶点 {i} 存在零长度边，跳过切割")continue# 步骤4：计算两条边上的分割点（从当前顶点截取最短边长度）s_prev = get_split_point(p_curr, p_prev, min_len)  # 前一条边上的分割点s_next = get_split_point(p_curr, p_next, min_len)  # 后一条边上的分割点# 标记分割点到调试图层mid_cursor.insertRow([s_prev, "分割点"])mid_cursor.insertRow([s_next, "分割点"])# 步骤5：创建切割后的三角形，并处理面积超限问题try:# 构建三角形顶点（需闭合：起点=终点）triangle_points = [p_curr, s_prev, s_next, p_curr]# 创建三角形多边形（指定坐标系）triangle = arcpy.Polygon(arcpy.Array(triangle_points), sr)area = triangle.area  # 计算三角形面积# 若面积超过max_area，按比例缩小分割点（确保面积合规）if area > max_area:scale = math.sqrt(max_area / area)  # 缩放比例（面积比开平方=长度比）# 计算缩放后的分割点（以当前顶点为原点缩放）s_prev_scaled = arcpy.Point(p_curr.X + (s_prev.X - p_curr.X) * scale,p_curr.Y + (s_prev.Y - p_curr.Y) * scale)s_next_scaled = arcpy.Point(p_curr.X + (s_next.X - p_curr.X) * scale,p_curr.Y + (s_next.Y - p_curr.Y) * scale)# 重新创建缩小后的三角形triangle_points = [p_curr, s_prev_scaled, s_next_scaled, p_curr]triangle = arcpy.Polygon(arcpy.Array(triangle_points), sr)area = triangle.areaarcpy.AddMessage(f"顶点 {i} 面积超限，已缩放至 {area:.2f}㎡")# 步骤6：将切割结果写入切割图层cut_cursor.insertRow([triangle, oid, i, round(area, 4), round(min_len, 4)])arcpy.AddMessage(f"顶点 {i} 已切割：面积={area:.2f}㎡，最短边={min_len:.2f}m")except Exception as e:arcpy.AddWarning(f"顶点 {i} 创建三角形失败：{str(e)}")# 处理完成，输出结果路径arc

（六）输出图层设计

脚本自动创建 3 个输出图层，满足不同分析需求：

图层名称	要素类型	核心字段	用途
顶点标记图层	点	ORIG_FID（原要素 ID）、INT_ANGLE（内角值）	标记锐角顶点位置及角度，用于后续核查
切割结果图层	多边形	ORIG_FID、VERTEX_ID（顶点 ID）、AREA（面积）、MIN_EDGE（最短边长度）	存储切割生成的三角形，记录关键属性
中间点图层	点	TYPE（类型：顶点 / 分割点）	调试用图层，查看分割点位置是否合理

三、完整代码

# -*- coding: utf-8 -*-
import arcpy
import os
import sys
from arcpy import env# 环境设置
env.overwriteOutput = True
env.XYResolution = "0.001 Meters"
env.XYTolerance = "0.001 Meters"# 临时文件路径
temp_workspace = r"C:\Users\Administrator\Desktop\test\临时文件"
overlap_temp = os.path.join(temp_workspace, "临时重叠区.shp")
remain_temp = os.path.join(temp_workspace, "临时剩余区.shp")
env.scratchWorkspace = temp_workspace# 自动创建临时文件夹
if not os.path.exists(temp_workspace):try:os.makedirs(temp_workspace)arcpy.AddMessage(f"自动创建临时文件夹：{temp_workspace}")except Exception as e:arcpy.AddError(f"创建临时文件夹失败：{str(e)}")sys.exit(1)# 输入输出参数
target_shp = r"C:\Users\Administrator\Desktop\test\01测试数据\测试数据1\锐角测试.shp"
tool_shp = r"C:\Users\Administrator\Desktop\test\锐角切割.shp"
output_split = r"C:\Users\Administrator\Desktop\test\最终分割结果.shp"# 基础检查
check_list = [(arcpy.Exists(target_shp), f"目标图层不存在：{target_shp}"),(arcpy.Exists(tool_shp), f"工具图层不存在：{tool_shp}"),(arcpy.Describe(target_shp).spatialReference.type == "Projected", f"目标图层不是投影坐标系：{arcpy.Describe(target_shp).spatialReference.name}"),(arcpy.Describe(tool_shp).spatialReference.name == arcpy.Describe(target_shp).spatialReference.name, f"坐标系不一致：目标={arcpy.Describe(target_shp).spatialReference.name}，工具={arcpy.Describe(tool_shp).spatialReference.name}")
]for is_valid, error_msg in check_list:if not is_valid:print(f"检查失败：{error_msg}")arcpy.AddError(error_msg)sys.exit(1)# 核心分割操作
try:# 清理旧临时文件for temp in [overlap_temp, remain_temp]:if arcpy.Exists(temp):arcpy.Delete_management(temp)# 步骤1：提取重叠区域arcpy.AddMessage("提取重叠区域...")arcpy.Intersect_analysis(in_features=[target_shp, tool_shp],out_feature_class=overlap_temp,join_attributes="ALL")arcpy.AddField_management(overlap_temp, "SPLIT_TYPE", "TEXT", 20)with arcpy.da.UpdateCursor(overlap_temp, ["SPLIT_TYPE"]) as cur:for row in cur:row[0] = "重叠切割部分"cur.updateRow(row)# 步骤2：提取剩余区域arcpy.AddMessage("提取剩余区域...")arcpy.Erase_analysis(in_features=target_shp,erase_features=tool_shp,out_feature_class=remain_temp,cluster_tolerance="0.001 Meters")arcpy.AddField_management(remain_temp, "SPLIT_TYPE", "TEXT", 20)with arcpy.da.UpdateCursor(remain_temp, ["SPLIT_TYPE"]) as cur:for row in cur:row[0] = "非重叠剩余部分"cur.updateRow(row)# 步骤3：合并结果（修正工具名称为Merge_management）arcpy.AddMessage("合并分割结果...")arcpy.Merge_management(  # 此处为修正点inputs=[overlap_temp, remain_temp],output=output_split)# 清理临时文件for temp in [overlap_temp, remain_temp]:if arcpy.Exists(temp):arcpy.Delete_management(temp)arcpy.AddMessage(f"\n分割完成！结果路径：{output_split}")arcpy.AddMessage("可通过 'SPLIT_TYPE' 字段区分两部分")except Exception as e:print(f"分割过程出错：{str(e)}")arcpy.AddError(f"分割出错：{str(e)}")for temp in [overlap_temp, remain_temp]:if arcpy.Exists(temp):arcpy.Delete_management(temp)sys.exit(1)