（三十三）arcpy开发&三调锐角检测工具开发

arcgis开发python视频教程

https://edu.youkuaiyun.com/course/detail/25535

最近在群里看到了小伙伴们在做三调，有需要锐角检测这样的工具。于是试着写了一个，这是利用arcpy写的一个锐角检测工具，最后做成界面如下图所示。

那么代码是如何实现呢。首先是对每个面进行遍历，分别找出面中的顶点。对顶点进行编号，然后就是利用平面向量乘积等数学相关内容，来计算角度。哎，估计时间过得太久了，一些高中的数学知识都忘记了，不禁感叹，时间过得太快，该忘的都忘了。不过这里要提示大家的是，这里使用的是投影坐标系，一般的直接坐标来计算角度的。就是使用简单向量，没有用到了大地测量复杂的运算，我想对于地球上一个非常小的图斑用简单的直角坐标计算值应该是和大地测量中地理坐标非常相近了，误差几乎可以忽略。

好了，小编这里还是将代码给公布出来了。当然，其中需要测试，也希望对大家有借鉴的地方，不足之处也可以提出来。因为平时也要忙其他的，不能及时回复，实在抱歉。好了，来看一下实现代码。

#encoding: utf-8
import arcpy
import math

class Point:
    X=0
    Y=0
def CaculateAngle(pointArray,angle):
    if len(pointArray)<3:
        return False
    for i in range(0, len(pointArray)):
        print i, pointArray[i]
        ptL=Point()
        ptM=Point()
        ptR=Point()
        if i<len(pointArray)-2:
         ptL=pointArray[i]
         ptM=pointArray[i+1]
         ptR=pointArray[i+2]
        elif i==len(pointArray)-2:
            ptL = pointArray[i]
            ptM = pointArray[i + 1]
            ptR = pointArray[0]
        elif i == len(pointArray) - 1:
            ptL = pointArray[i]
            ptM = pointArray[0]
            ptR = pointArray[1]
        AML=Point()
        AML.X=ptL.X-ptM.X
        AML.Y=ptL.Y-ptM.Y

        AMR=Point()
        AMR.X=ptR.X-ptM.X
        AMR.Y = ptR.Y - ptM.Y

        ab=(AML.X*AMR.X+AML.Y*AMR.Y)
        M=(math.sqrt(AML.X*AML.X+AML.Y*AML.Y)*math.sqrt(AMR.X*AMR.X+AMR.Y*AMR.Y));
        cosA=ab/M

        #向量乘积小于0，说明是大于90度，不用判断
        if ab<0:
            continue
        #因为0到π/2是递减,所以小于cos给定角度值，那么角度比给定角度值大,就返回了
        elif cosA<=math.cos(angle/180*math.pi):
            continue
        elif cosA>math.cos(angle/180*math.pi):
            return True
    # 遍历所有的顶点了，说明满足要求
    return False


testData=ur"D:\Data\PolygonData.shp"
inpath = 'D:\\Data\\中国国界和省界的SHP格式数据\\省界\\test.txt'
uipath = unicode(inpath, "utf8")
outFile = open(uipath, "w")
for row in arcpy.da.SearchCursor(testData, ["FID", "SHAPE@"]):
    ID = row[0]
    tempPart = []
    partnum = 0
    for part in row[1]:
        print("Part {0}:".format(partnum))
        tempPart = []
        ptnum=0
        pointArray=[]
        for pnt in part:
         if pnt:
            print("点 {0}:".format(ptnum))
            point = Point()
            point.X = pnt.X
            point.Y = pnt.Y
            pointArray.append(point)
            ptnum+=1
         else:
              print("内环:")
        pointArray.pop(len(pointArray)-1)

        #将不符合要求的写入到txt文件中
        if CaculateAngle(pointArray, 45):
            outFile.write("{0}\n".format(ID))
            print("ID为{0}的属性面含有不满足要求的锐角\n".format(ID))
    partnum += 1

---

---

                            更多内容，请微信扫二维码关注公众号，或者加入arcpy开发qq学习群：487352121