ArcGIS 获得三角形内切圆

如何求三角形的内心坐标

来自：https://zhidao.baidu.com/question/158982463.html

三角形边长为a,b,c,顶点为A(x1,y1)B(x2,y2)C(x3,y3),并给出证明

内心是角平分线的交点，到三边距离相等.
设:在三角形ABC中,三顶点的坐标为：A(x1,y1),B(x2,y2),C(x3,y3) BC=a,CA=b,AB=c
内心为M （X，Y）则有aMA+bMB+cMC=0（三个向量）
MA=（X1-X，Y1-Y）
MB=(X2-X,Y2-Y)
MC=(X3-X,Y3-Y)
则：a(X1-X)+b(X2-X)+c(X3-X)=0，a(Y1-Y)+b(Y2-Y)+c(Y3-Y)=0
∴X=（aX1+bX2+cX3)/(a+b+c)，Y=（aY1+bY2+cY3)/(a+b+c)
∴M（（aX1+bX2+cX3)/(a+b+c)，（aY1+bY2+cY3)/(a+b+c)），真的很麻烦，有具体数的话，你可以求两条角分线，再求交点

三角形内切圆半径公式：r=2S/(a+b+c）

关注微信公众号，学习更多

import arcpy

def FieldExists(TableName,FieldName):
    desc = arcpy.Describe(TableName)
    FieldName=FieldName.upper()
    for field in desc.fields:
        if field.Name.upper() ==FieldName:
            return True
            break
    return False
#获得点的距离
def pointDistance(pt1,pt2):
    return math.sqrt((pt1.X-pt2.X)*(pt1.X-pt2.X)+(pt1.Y-pt2.Y)*(pt1.Y-pt2.Y))

#获得所有的点
def getallpoint(geo):
    parray = arcpy.Array()
    for part in geo:
        # Print the part number
        # print("Part {}:".format(partnum))
        # Step through each vertex in the feature
        for pnt in part:
            if pnt == None:
                # If pnt is None, this represents an interior ring
                print("Interior Ring:")
                #n = n + 1
            else:
                parray.add(pnt)
    return parray

infc = arcpy.GetParameterAsText(0)
outfc = arcpy.GetParameterAsText(1)
arcpy.AddMessage("输入："+infc)
arcpy.AddMessage("输出："+outfc)
FieldName="X"
if not FieldExists(infc,FieldName):
    arcpy.AddMessage("添加字段：" + FieldName)
    arcpy.management.AddField(infc,FieldName , "DOUBLE")
FieldName="Y"
if not FieldExists(infc,FieldName):
    arcpy.AddMessage("添加字段：" + FieldName)
    arcpy.management.AddField(infc,FieldName , "DOUBLE")
FieldName="R"
if not FieldExists(infc,FieldName):
    arcpy.AddMessage("添加字段："+FieldName)
    arcpy.management.AddField(infc,FieldName , "DOUBLE")


with arcpy.da.UpdateCursor(infc, ["OID@", "SHAPE@","X","Y","R","SHAPE@AREA"]) as cursor:
    for row in cursor:
        parray=getallpoint(row[1])
        pnum=len(parray)
        if pnum>4:
            arcpy.AddMessage("节点大于4"+str(pnum)+",FID="+str(row[0]))
            continue

        pt1=parray[0]
        pt2=parray[1]
        pt3=parray[2]
        x1=pt1.X
        y1=pt1.Y

        x2=pt2.X
        y2=pt2.Y

        x3=pt3.X
        y3=pt3.Y
        a=pointDistance(pt2,pt3)
        b=pointDistance(pt1,pt3)
        c=pointDistance(pt1,pt2)
        x=(a*x1+b*x2+c*x3)/(a+b+c)
        y=(a*y1+b*y2+c*y3)/(a+b+c)
        s=row[5] #(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2.0
        r=2.0*s/(a+b+c)
        row[2]=x
        row[3] =y
        row[4] =r
        cursor.updateRow(row)
out_layer="outyl"
sr = arcpy.Describe(infc).spatialReference
arcpy.MakeXYEventLayer_management(infc, "x", "y", out_layer, spatial_reference=sr)

arcpy.Buffer_analysis(out_layer, outfc, buffer_distance_or_field="R")