ArcPy计算尖锐角

# -*-coding:gbk-*-

import arcpy.da
import json
import math


def cal_angle_a(ax, ay, bx, by, cx, cy):
    len_ab = math.sqrt(math.pow(ax - bx, 2) + math.pow(ay - by, 2))
    len_ac = math.sqrt(math.pow(ax - cx, 2) + math.pow(ay - cy, 2))
    len_bc = math.sqrt(math.pow(bx - cx, 2) + math.pow(by - cy, 2))
    cos_a = (math.pow(len_ab, 2) + math.pow(len_ac, 2) - math.pow(len_bc, 2)) / (2 * len_ab * len_ac)
    return 0 if cos_a > 0.99 else (180 if cos_a < -0.99 else (math.acos(cos_a) * 180 / math.pi))


if __name__ == '__main__':
    try:
        features = r'E:/Temp/SD_Temp.gdb/VEGA_J'
        min_angle = 80
        line_fields = ['SHAPE@JSON', 'OBJECTID']  # 图形XY坐标
        with arcpy.da.SearchCursor(features, line_fields) as line_cursor:
            for line_row in line_cursor:
                checked = False
                coord_array = json.loads(line_row[0])['rings']
                for m in range(0, len(coord_array)):
                    xys = coord_array[m]
                    if len(xys) > 2:
                        for n in range(1, len(xys) - 2):
                            angle_a = cal_angle_a(xys[n][0], xys[n][1], xys[n - 1][0], xys[n - 1][1], xys[n + 1][0],
                                                  xys[n + 1][1])
                            if angle_a < min_angle:
                                print 'Small angle at {0}'.format(line_row[1])
                                checked = True
                                break
                        if checked:
                            break
                    if checked:
                        break
                if not checked:
                    print 'Skip {0}'.format(line_row[1])
    except Exception, e:
        print e.message